Merge pull request #1207 from aprokop/segment_intersection

Add segment-segment intersection (2D only)

src/geometry/ArborX_Segment.hpp

@@ -22,8 +22,8 @@ namespace ArborX::Experimental
 template <int DIM, typename Coordinate = float>
 struct Segment
 {
-  ArborX::Point<DIM, Coordinate> _start;
-  ArborX::Point<DIM, Coordinate> _end;
+  ArborX::Point<DIM, Coordinate> a;
+  ArborX::Point<DIM, Coordinate> b;
 };
 
 template <int DIM, typename Coordinate>

src/geometry/algorithms/ArborX_Centroid.hpp

@@ -119,7 +119,7 @@ struct centroid<SegmentTag, Segment>
     // WARNING implicit requirement on KDOP first DIM directions
     Point<DIM, Coordinate> point;
     for (int d = 0; d < DIM; ++d)
-      point[d] = (segment._start[d] + segment._end[d]) / 2;
+      point[d] = (segment.a[d] + segment.b[d]) / 2;
     return point;
   }
 };

src/geometry/algorithms/ArborX_Distance.hpp

@@ -255,20 +255,20 @@ struct distance<PointTag, SegmentTag, Point, Segment>
     constexpr int DIM = GeometryTraits::dimension_v<Point>;
     using Coordinate = GeometryTraits::coordinate_type_t<Point>;
 
-    if (Details::equals(segment._start, segment._end))
-      return Details::distance(point, segment._start);
+    if (Details::equals(segment.a, segment.b))
+      return Details::distance(point, segment.a);
 
-    auto const dir = segment._end - segment._start;
+    auto const dir = segment.b - segment.a;
 
     // The line of the segment [a,b] is parametrized as a + t * (b - a).
     // Find the projection of the point to that line, and clamp it.
     auto t =
-        Kokkos::clamp(dir.dot(point - segment._start) / dir.dot(dir),
+        Kokkos::clamp(dir.dot(point - segment.a) / dir.dot(dir),
                       static_cast<Coordinate>(0), static_cast<Coordinate>(1));
 
     Point projection;
     for (int d = 0; d < DIM; ++d)
-      projection[d] = segment._start[d] + t * dir[d];
+      projection[d] = segment.a[d] + t * dir[d];
 
     return Details::distance(point, projection);
   }

src/geometry/algorithms/ArborX_Expand.hpp

@@ -230,8 +230,8 @@ struct expand<BoxTag, SegmentTag, Box, Segment>
 {
   KOKKOS_FUNCTION static void apply(Box &box, Segment const &segment)
   {
-    Details::expand(box, segment._start);
-    Details::expand(box, segment._end);
+    Details::expand(box, segment.a);
+    Details::expand(box, segment.b);
   }
 };
 

src/geometry/algorithms/ArborX_Intersects.hpp

@@ -412,6 +412,70 @@ struct intersects<PointTag, TetrahedronTag, Point, Tetrahedron>
   }
 };
 
+template <typename Segment>
+struct intersects<SegmentTag, SegmentTag, Segment, Segment>
+{
+  // The algorithm is described in
+  // https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/
+
+  // Given a segment and a collinear point, check if the point lies on the
+  // segment
+  template <typename Point>
+  KOKKOS_FUNCTION static bool on_segment(Segment const &segment, Point const &p)
+  {
+    using Kokkos::max;
+    using Kokkos::min;
+    auto const &a = segment.a;
+    auto const &b = segment.b;
+    return (p[0] >= min(a[0], b[0]) && p[0] <= max(a[0], b[0])) &&
+           (p[1] >= min(a[1], b[1]) && p[1] <= max(a[1], b[1]));
+  }
+
+  // Find orientation of an ordered tuple (a, b, c)
+  // -  0: points are collinear
+  // -  1: clockwise
+  // - -1: counter-clockwise
+  template <typename Point>
+  KOKKOS_FUNCTION static int orientation(Point const &a, Point const &b,
+                                         Point const &c)
+  {
+    auto x = (b[1] - a[1]) * (c[0] - b[0]) - (b[0] - a[0]) * (c[1] - b[1]);
+    return (0 < x) - (x < 0); // sgn
+  }
+
+  KOKKOS_FUNCTION static constexpr bool apply(Segment const &segment1,
+                                              Segment const &segment2)
+  {
+    static_assert(GeometryTraits::dimension_v<Segment> == 2);
+
+    int o1 = orientation(segment1.a, segment1.b, segment2.a);
+    int o2 = orientation(segment1.a, segment1.b, segment2.b);
+    int o3 = orientation(segment2.a, segment2.b, segment1.a);
+    int o4 = orientation(segment2.a, segment2.b, segment1.b);
+
+    // General case (no collinearity)
+    if (o1 != o2 && o3 != o4)
+      return true;
+
+    // Special cases
+
+    // segment2.a is collinear to segment1 and is within
+    if (o1 == 0 && on_segment(segment1, segment2.a))
+      return true;
+    // segment2.b is collinear to segment1 and is within
+    if (o2 == 0 && on_segment(segment1, segment2.b))
+      return true;
+    // segment1.a is collinear to segment2 and is within
+    if (o3 == 0 && on_segment(segment2, segment1.a))
+      return true;
+    // segment1.b is collinear to segment2 and is within
+    if (o4 == 0 && on_segment(segment2, segment1.b))
+      return true;
+
+    return false;
+  }
+};
+
 } // namespace Dispatch
 
 } // namespace ArborX::Details

test/tstDetailsAlgorithms.cpp

@@ -313,6 +313,20 @@ BOOST_AUTO_TEST_CASE(intersects)
   BOOST_TEST(!intersects(Point{0, 0, 1.1}, tet));
   BOOST_TEST(!intersects(Point{0.5, 0.5, 0.5}, tet));
   BOOST_TEST(!intersects(Point{-0.5, 0.5, 0.5}, tet));
+
+  // segment
+  using Segment2 = ArborX::Experimental::Segment<2>;
+  constexpr Segment2 seg{{1, 1}, {2, 2}};
+  BOOST_TEST(intersects(Segment2{{2, 2}, {3, 3}}, seg));
+  BOOST_TEST(intersects(Segment2{{1.5, 1.5}, {1.7, 1.7}}, seg));
+  BOOST_TEST(intersects(Segment2{{0, 0}, {1, 1}}, seg));
+  BOOST_TEST(intersects(Segment2{{1, 2}, {2, 1}}, seg));
+  BOOST_TEST(intersects(Segment2{{2, 0}, {0, 2}}, seg));
+  BOOST_TEST(intersects(Segment2{{1, 3}, {3, 1}}, seg));
+  BOOST_TEST(!intersects(Segment2{{0, 0}, {0.9, 0.9}}, seg));
+  BOOST_TEST(!intersects(Segment2{{1.1, 1}, {2, 1}}, seg));
+  BOOST_TEST(!intersects(Segment2{{1, 0}, {2, 1}}, seg));
+  BOOST_TEST(!intersects(Segment2{{1, 3}, {3, 1.1}}, seg));
 }
 
 BOOST_AUTO_TEST_CASE(equals)