Fix a bug in segment intersection primitive where two collinear segment touch at endpoints is miscomputed as a degenerate segment

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

@@ -137,36 +137,40 @@ __forceinline__ T __device__ squared_segment_distance(vec_2d<T> const& a,
 
 /**
  * @internal
- * @brief Given two collinear or parallel segments, return their potential overlapping segment
+ * @brief Given two collinear or parallel segments, return their potential overlapping segment or
+ * point
  *
  * @p a, @p b, @p c, @p d refer to end points of segment ab and cd.
  * @p center is the geometric center of the segments, used to decondition the coordinates.
  *
- * @return optional end points of overlapping segment
+ * @return A pair of optional overlapping point or segments
  */
 template <typename T>
-__forceinline__ thrust::optional<segment<T>> __device__ collinear_or_parallel_overlapping_segments(
-  vec_2d<T> a, vec_2d<T> b, vec_2d<T> c, vec_2d<T> d, vec_2d<T> center = vec_2d<T>{})
+__forceinline__
+
+  thrust::pair<thrust::optional<vec_2d<T>>, thrust::optional<segment<T>>>
+    __device__ collinear_or_parallel_overlapping_segments(
+      vec_2d<T> a, vec_2d<T> b, vec_2d<T> c, vec_2d<T> d, vec_2d<T> center = vec_2d<T>{})
 {
   auto ab = b - a;
   auto ac = c - a;
 
   // Parallel
-  if (not float_equal(det(ab, ac), T{0})) return thrust::nullopt;
+  if (not float_equal(det(ab, ac), T{0})) return {thrust::nullopt, thrust::nullopt};
 
   // Must be on the same line, sort the endpoints
   if (b < a) thrust::swap(a, b);
   if (d < c) thrust::swap(c, d);
 
   // Test if not overlap
-  if (b < c || d < a) return thrust::nullopt;
+  if (b < c || d < a) return {thrust::nullopt, thrust::nullopt};
 
   // Compute smallest interval between the segments
   auto e0 = a > c ? a : c;
   auto e1 = b < d ? b : d;
 
-  // Decondition the coordinates
-  return segment<T>{e0 + center, e1 + center};
+  if (e0 == e1) { return {e0 + center, thrust::nullopt}; }
+  return {thrust::nullopt, segment<T>{e0 + center, e1 + center}};
 }
 
 /**
@@ -192,7 +196,7 @@ segment_intersection(segment<T> const& segment1, segment<T> const& segment2)
 
   if (float_equal(denom, T{0})) {
     // Segments parallel or collinear
-    return {thrust::nullopt, collinear_or_parallel_overlapping_segments(a, b, c, d, center)};
+    return collinear_or_parallel_overlapping_segments(a, b, c, d, center);
   }
 
   auto ac               = c - a;

cpp/include/cuspatial/geometry/vec_2d.hpp

@@ -17,6 +17,7 @@
 #pragma once
 
 #include <cuspatial/cuda_utils.hpp>
+#include <cuspatial/detail/utility/floating_point.cuh>
 
 #include <algorithm>
 #include <ostream>
@@ -58,7 +59,7 @@ class alignas(2 * sizeof(T)) vec_2d {
    */
   friend bool CUSPATIAL_HOST_DEVICE operator==(vec_2d<T> const& lhs, vec_2d<T> const& rhs)
   {
-    return (lhs.x == rhs.x) && (lhs.y == rhs.y);
+    return detail::float_equal<T>(lhs.x, rhs.x) && detail::float_equal(lhs.y, rhs.y);
   }
 
   /**

cpp/tests/operators/linestrings_test.cu

@@ -133,6 +133,32 @@ TYPED_TEST(SegmentIntersectionTest, IntersectAtEndPoint)
   run_single_intersection_test(ab, cd, points_expected, segments_expected);
 }
 
+TYPED_TEST(SegmentIntersectionTest, IntersectAtEndPoint2)
+{
+  using T = TypeParam;
+
+  segment<T> ab{{-1.0, 0.0}, {0.0, 0.0}};
+  segment<T> cd{{0.0, 0.0}, {0.0, 1.0}};
+
+  std::vector<thrust::optional<vec_2d<T>>> points_expected{vec_2d<T>{0.0, 0.0}};
+  std::vector<thrust::optional<segment<T>>> segments_expected{thrust::nullopt};
+
+  run_single_intersection_test(ab, cd, points_expected, segments_expected);
+}
+
+TYPED_TEST(SegmentIntersectionTest, IntersectAtEndPoint3)
+{
+  using T = TypeParam;
+
+  segment<T> ab{{-1.0, 0.0}, {0.0, 0.0}};
+  segment<T> cd{{1.0, 0.0}, {0.0, 0.0}};
+
+  std::vector<thrust::optional<vec_2d<T>>> points_expected{vec_2d<T>{0.0, 0.0}};
+  std::vector<thrust::optional<segment<T>>> segments_expected{thrust::nullopt};
+
+  run_single_intersection_test(ab, cd, points_expected, segments_expected);
+}
+
 TYPED_TEST(SegmentIntersectionTest, UnparallelDisjoint1)
 {
   using T = TypeParam;

cpp/tests/spatial/intersection/linestring_intersection_test.cu

@@ -278,6 +278,28 @@ TYPED_TEST(LinestringIntersectionTest, SingletoSingleOnePair)
                      expected);
 }
 
+TYPED_TEST(LinestringIntersectionTest, OnePairWithRings)
+{
+  using T = TypeParam;
+  using P = vec_2d<T>;
+
+  using index_t = typename linestring_intersection_result<T, std::size_t>::index_t;
+  using types_t = typename linestring_intersection_result<T, std::size_t>::types_t;
+
+  auto multilinestrings1 = make_multilinestring_array<T>({0, 1}, {0, 2}, {{-1, 0}, {0, 0}});
+
+  auto multilinestrings2 =
+    make_multilinestring_array<T>({0, 1}, {0, 5}, {{0, 0}, {0, 1}, {1, 1}, {1, 0}, {0, 0}});
+
+  auto expected = make_linestring_intersection_result<T, index_t, types_t>(
+    {0, 1}, {0}, {0}, {P{0.0, 0.0}}, {}, {0}, {0}, {0}, {0}, this->stream(), this->mr());
+
+  CUSPATIAL_RUN_TEST(this->template run_single_test<index_t>,
+                     multilinestrings1.range(),
+                     multilinestrings2.range(),
+                     expected);
+}
+
 TYPED_TEST(LinestringIntersectionTest, SingletoSingleOnePairWithDuplicatePoint)
 {
   using T = TypeParam;

python/cuspatial/cuspatial/tests/binops/test_intersections.py

@@ -88,9 +88,9 @@ def test_one_pair_with_overlap():
     expect_ids = pd.DataFrame(
         {
             "lhs_linestring_id": [[0]],
-            "lhs_segment_id": [[0]],
+            "lhs_segment_id": [[1]],
             "rhs_linestring_id": [[0]],
-            "rhs_segment_id": [[0]],
+            "rhs_segment_id": [[1]],
         }
     )
 
@@ -122,9 +122,9 @@ def test_two_pairs_with_intersect_and_overlap():
     expect_ids = pd.DataFrame(
         {
             "lhs_linestring_id": [[0], [0, 0]],
-            "lhs_segment_id": [[0], [1, 0]],
+            "lhs_segment_id": [[1], [1, 0]],
             "rhs_linestring_id": [[0], [0, 0]],
-            "rhs_segment_id": [[0], [0, 2]],
+            "rhs_segment_id": [[1], [0, 2]],
         }
     )