Sort the mergeable segments before computing merged segments (#1207)

This PR is part-1 fix to \#1200. In `find_and_combine_segments` the N^2 algorithm depends on the fact that the API needs to be presorted with certain criteria. This Pr adds such sorting.

Authors:
  - Michael Wang (https://github.com/isVoid)

Approvers:
  - Mark Harris (https://github.com/harrism)
  - H. Thomson Comer (https://github.com/thomcom)

URL: #1207

cpp/include/cuspatial/detail/algorithm/is_point_in_polygon.cuh

@@ -68,8 +68,8 @@ __device__ inline bool is_point_in_polygon(vec_2d<T> const& test_point, PolygonR
       T run_to_point  = test_point.x - a.x;
 
       // point-on-edge test
-      bool is_colinear = float_equal(run * rise_to_point, run_to_point * rise);
-      if (is_colinear) {
+      bool is_collinear = float_equal(run * rise_to_point, run_to_point * rise);
+      if (is_collinear) {
         T minx = a.x;
         T maxx = b.x;
         if (minx > maxx) thrust::swap(minx, maxx);

cpp/include/cuspatial/detail/find/find_and_combine_segment.cuh

@@ -19,20 +19,33 @@
 #include <cuspatial/cuda_utils.hpp>
 #include <cuspatial/detail/utility/linestring.cuh>
 #include <cuspatial/error.hpp>
+#include <cuspatial/iterator_factory.cuh>
 
 #include <rmm/cuda_stream_view.hpp>
 #include <rmm/exec_policy.hpp>
 
 #include <ranger/ranger.hpp>
 
+#include <thrust/sort.h>
 #include <thrust/uninitialized_fill.h>
 
 namespace cuspatial {
 namespace detail {
 
 /**
  * @internal
- * @brief Kernel to merge segments, naive n^2 algorithm.
+ * @pre All segments in range @p segments , it is presorted using `segment_comparator`.
+ *
+ * @brief Kernel to merge segments. Each thread works on one segment space, within each space,
+ * `segment_comparator` guarantees that segments with same slope are grouped together. We call
+ * each of such group is a "mergeable group". Within each mergeable group, the first segment is
+ * the "leading segment". The algorithm behave as follows:
+ *
+ * 1. For each mergeable group, loop over the rest of the segments in the group and see if it is
+ * mergeable with the leading segment. If it is, overwrite the leading segment with the merged
+ * result. Then mark the segment as merged by setting the flag to 1. This makes sure the inner loop
+ * for each merged segment is not run again.
+ * 2. Repeat 1 until all mergeable group is processed.
  */
 template <typename OffsetRange, typename SegmentRange, typename OutputIt>
 void __global__ simple_find_and_combine_segments_kernel(OffsetRange offsets,
@@ -45,6 +58,9 @@ void __global__ simple_find_and_combine_segments_kernel(OffsetRange offsets,
       merged_flag[i] = 0;
     }
 
+    // For each of the segment, loop over the rest of the segment in the space and see
+    // if it is mergeable with the current segment.
+    // Note that if the current segment is already merged. Skip checking.
     for (auto i = offsets[pair_idx]; i < offsets[pair_idx + 1] && merged_flag[i] != 1; i++) {
       for (auto j = i + 1; j < offsets[pair_idx + 1]; j++) {
         auto res = maybe_merge_segments(segments[i], segments[j]);
@@ -58,20 +74,66 @@ void __global__ simple_find_and_combine_segments_kernel(OffsetRange offsets,
   }
 }
 
+/**
+ * @brief Comparator for sorting the segment range.
+ *
+ * This comparator makes sure that the segment range is sorted such that:
+ * 1. Segments with the same space id are grouped together.
+ * 2. Segments within the same space are grouped by their slope.
+ * 3. Within each slope group, segments are sorted by their lower left point.
+ */
+template <typename index_t, typename T>
+struct segment_comparator {
+  bool __device__ operator()(thrust::tuple<index_t, segment<T>> const& lhs,
+                             thrust::tuple<index_t, segment<T>> const& rhs) const
+  {
+    auto lhs_index   = thrust::get<0>(lhs);
+    auto rhs_index   = thrust::get<0>(rhs);
+    auto lhs_segment = thrust::get<1>(lhs);
+    auto rhs_segment = thrust::get<1>(rhs);
+
+    // Compare space id
+    if (lhs_index == rhs_index) {
+      // Compare slope
+      if (lhs_segment.collinear(rhs_segment)) {
+        // Sort by the lower left point of the segment.
+        return lhs_segment.lower_left() < rhs_segment.lower_left();
+      }
+      return lhs_segment.slope() < rhs_segment.slope();
+    }
+    return lhs_index < rhs_index;
+  }
+};
+
 /**
  * @internal
  * @brief For each pair of mergeable segment, overwrites the first segment with merged result,
  * sets the flag for the second segment as 1.
+ *
+ * @note This function will alter the input segment range by rearranging the order of the segments
+ * within each space so that merging kernel can take place.
  */
 template <typename OffsetRange, typename SegmentRange, typename OutputIt>
 void find_and_combine_segment(OffsetRange offsets,
                               SegmentRange segments,
                               OutputIt merged_flag,
                               rmm::cuda_stream_view stream)
 {
+  using index_t   = typename OffsetRange::value_type;
+  using T         = typename SegmentRange::value_type::value_type;
   auto num_spaces = offsets.size() - 1;
   if (num_spaces == 0) return;
 
+  // Construct a key iterator using the offsets of the segment and the segment itself.
+  auto space_id_iter         = make_geometry_id_iterator<index_t>(offsets.begin(), offsets.end());
+  auto space_id_segment_iter = thrust::make_zip_iterator(space_id_iter, segments.begin());
+
+  thrust::sort_by_key(rmm::exec_policy(stream),
+                      space_id_segment_iter,
+                      space_id_segment_iter + segments.size(),
+                      segments.begin(),
+                      segment_comparator<index_t, T>{});
+
   auto [threads_per_block, num_blocks] = grid_1d(num_spaces);
   simple_find_and_combine_segments_kernel<<<num_blocks, threads_per_block, 0, stream.value()>>>(
     offsets, segments, merged_flag);

cpp/include/cuspatial/detail/utility/linestring.cuh

@@ -93,10 +93,10 @@ __forceinline__ T __device__ point_to_segment_distance_squared(vec_2d<T> const&
  * @brief Computes shortest distance between two segments (ab and cd) that don't intersect.
  */
 template <typename T>
-__forceinline__ T __device__ segment_distance_no_intersect_or_colinear(vec_2d<T> const& a,
-                                                                       vec_2d<T> const& b,
-                                                                       vec_2d<T> const& c,
-                                                                       vec_2d<T> const& d)
+__forceinline__ T __device__ segment_distance_no_intersect_or_collinear(vec_2d<T> const& a,
+                                                                        vec_2d<T> const& b,
+                                                                        vec_2d<T> const& c,
+                                                                        vec_2d<T> const& d)
 {
   auto dist_sqr = min(
     min(point_to_segment_distance_squared(a, c, d), point_to_segment_distance_squared(b, c, d)),
@@ -123,7 +123,7 @@ __forceinline__ T __device__ squared_segment_distance(vec_2d<T> const& a,
 
   if (float_equal(denom, T{0})) {
     // Segments parallel or collinear
-    return segment_distance_no_intersect_or_colinear(a, b, c, d);
+    return segment_distance_no_intersect_or_collinear(a, b, c, d);
   }
 
   auto ac               = c - a;
@@ -132,7 +132,7 @@ __forceinline__ T __device__ squared_segment_distance(vec_2d<T> const& a,
   auto r                = r_numer * denom_reciprocal;
   auto s                = det(ac, ab) * denom_reciprocal;
   if (r >= 0 and r <= 1 and s >= 0 and s <= 1) { return 0.0; }
-  return segment_distance_no_intersect_or_colinear(a, b, c, d);
+  return segment_distance_no_intersect_or_collinear(a, b, c, d);
 }
 
 /**

cpp/include/cuspatial/geometry/segment.cuh

@@ -55,6 +55,20 @@ class alignas(sizeof(Vertex)) segment {
   /// Return the length squared of segment.
   T CUSPATIAL_HOST_DEVICE length2() const { return dot(v2 - v1, v2 - v1); }
 
+  /// Return slope of segment.
+  T CUSPATIAL_HOST_DEVICE slope() { return (v2.y - v1.y) / (v2.x - v1.x); }
+
+  /// Return the lower left vertex of segment.
+  Vertex CUSPATIAL_HOST_DEVICE lower_left() { return v1 < v2 ? v1 : v2; }
+
+  /// Returns true if two segments are on the same line
+  /// Test if the determinant of the matrix, of which the column vector is constructed from the
+  /// segments is 0.
+  bool CUSPATIAL_HOST_DEVICE collinear(segment<T> const& other)
+  {
+    return (v1.x - v2.x) * (other.v1.y - other.v2.y) == (v1.y - v2.y) * (other.v1.x - other.v2.x);
+  }
+
  private:
   friend std::ostream& operator<<(std::ostream& os, segment<T> const& seg)
   {

cpp/tests/find/find_and_combine_segments_test.cu

@@ -279,5 +279,67 @@ TYPED_TEST(FindAndCombineSegmentsTest, nooverlap3)
   CUSPATIAL_RUN_TEST(this->run_single_test,
                      segments,
                      {0, 0},
-                     {S{P{0.0, 0.0}, P{1.0, 1.0}}, S{P{0.0, 1.0}, P{1.0, 0.0}}});
+                     {S{P{0.0, 1.0}, P{1.0, 0.0}}, S{P{0.0, 0.0}, P{1.0, 1.0}}});
+}
+
+TYPED_TEST(FindAndCombineSegmentsTest, twospaces)
+{
+  using T       = TypeParam;
+  using index_t = std::size_t;
+  using P       = vec_2d<T>;
+  using S       = segment<T>;
+
+  auto segments = make_segment_array<index_t, T>({0, 2, 4},
+                                                 {S{P{0.0, 0.0}, P{1.0, 1.0}},
+                                                  S{P{1.0, 1.0}, P{2.0, 2.0}},
+                                                  S{P{1.0, 1.0}, P{0.0, 0.0}},
+                                                  S{P{2.0, 2.0}, P{1.0, 1.0}}});
+
+  CUSPATIAL_RUN_TEST(this->run_single_test,
+                     segments,
+                     {0, 1, 0, 1},
+                     {S{P{0.0, 0.0}, P{2.0, 2.0}},
+                      S{P{1.0, 1.0}, P{2.0, 2.0}},
+                      S{P{0.0, 0.0}, P{2.0, 2.0}},
+                      S{P{2.0, 2.0}, P{1.0, 1.0}}});
+}
+
+TYPED_TEST(FindAndCombineSegmentsTest, twospaces_non_contiguous_segments_with_empty)
+{
+  using T       = TypeParam;
+  using index_t = std::size_t;
+  using P       = vec_2d<T>;
+  using S       = segment<T>;
+
+  auto segments = make_segment_array<index_t, T>({0, 4, 4},
+                                                 {S{P{1.0, 1.0}, P{2.0, 2.0}},
+                                                  S{P{3.0, 3.0}, P{4.0, 4.0}},
+                                                  S{P{0.0, 0.0}, P{1.0, 1.0}},
+                                                  S{P{2.0, 2.0}, P{3.0, 3.0}}});
+
+  CUSPATIAL_RUN_TEST(this->run_single_test,
+                     segments,
+                     {0, 1, 1, 1},
+                     {S{P{0.0, 0.0}, P{4.0, 4.0}},
+                      S{P{1.0, 1.0}, P{2.0, 2.0}},
+                      S{P{2.0, 2.0}, P{3.0, 3.0}},
+                      S{P{3.0, 3.0}, P{4.0, 4.0}}});
+}
+
+TYPED_TEST(FindAndCombineSegmentsTest, onespace_non_contiguous_segments_overlaps)
+{
+  using T       = TypeParam;
+  using index_t = std::size_t;
+  using P       = vec_2d<T>;
+  using S       = segment<T>;
+
+  auto segments = make_segment_array<index_t, T>(
+    {0, 3},
+    {S{P{1.0, 1.0}, P{2.0, 2.0}}, S{P{4.0, 4.0}, P{5.0, 5.0}}, S{P{-1.0, -1.0}, P{4.0, 4.0}}});
+
+  CUSPATIAL_RUN_TEST(
+    this->run_single_test,
+    segments,
+    {0, 1, 1},
+    {S{P{-1.0, -1.0}, P{5.0, 5.0}}, S{P{1.0, 1.0}, P{2.0, 2.0}}, S{P{4.0, 4.0}, P{5.0, 5.0}}});
 }

cpp/tests/intersection/linestring_intersection_large_test.cu

@@ -27,6 +27,8 @@
 #include <cuspatial/geometry/vec_2d.hpp>
 #include <cuspatial/intersection.cuh>
 
+#include <rmm/exec_policy.hpp>
+
 #include <thrust/host_vector.h>
 
 template <typename T>

python/cuspatial/cuspatial/tests/binpreds/test_contains_properly.py

@@ -131,7 +131,7 @@ def test_float_precision_limits(point, polygon, expects):
     Unique success cases identified by @mharris. These go in a pair
     with test_float_precision_limits_failures because these are
     inconsistent results, where 0.6 fails above (as True, within the
-    polygon) and 0.66 below succeeds, though they are colinear.
+    polygon) and 0.66 below succeeds, though they are collinear.
     """
     gpdpoint = gpd.GeoSeries(point)
     gpdpolygon = gpd.GeoSeries(polygon)