Merge pull request #1165 from aprokop/add_segment

Add experimental Segment geometry

src/geometry/ArborX_DetailsAlgorithms.hpp

@@ -30,6 +30,7 @@ using GeometryTraits::BoxTag;
 using GeometryTraits::KDOPTag;
 using GeometryTraits::PointTag;
 using GeometryTraits::RayTag;
+using GeometryTraits::SegmentTag;
 using GeometryTraits::SphereTag;
 using GeometryTraits::TetrahedronTag;
 using GeometryTraits::TriangleTag;

src/geometry/ArborX_GeometryTraits.hpp

@@ -64,13 +64,15 @@ DEFINE_GEOMETRY(triangle, TriangleTag);
 DEFINE_GEOMETRY(kdop, KDOPTag);
 DEFINE_GEOMETRY(tetrahedron, TetrahedronTag);
 DEFINE_GEOMETRY(ray, RayTag);
+DEFINE_GEOMETRY(segment, SegmentTag);
 #undef DEFINE_GEOMETRY
 
 template <typename Geometry>
 inline constexpr bool is_valid_geometry =
     (is_point_v<Geometry> || is_box_v<Geometry> || is_sphere_v<Geometry> ||
      is_kdop_v<Geometry> || is_triangle_v<Geometry> ||
-     is_tetrahedron_v<Geometry> || is_ray_v<Geometry>);
+     is_tetrahedron_v<Geometry> || is_ray_v<Geometry> ||
+     is_segment_v<Geometry>);
 
 template <typename Geometry>
 using DimensionNotSpecializedArchetypeAlias =

src/geometry/ArborX_Segment.hpp

@@ -0,0 +1,130 @@
+/****************************************************************************
+ * Copyright (c) 2024 by the ArborX authors                                 *
+ * All rights reserved.                                                     *
+ *                                                                          *
+ * This file is part of the ArborX library. ArborX is                       *
+ * distributed under a BSD 3-clause license. For the licensing terms see    *
+ * the LICENSE file in the top-level directory.                             *
+ *                                                                          *
+ * SPDX-License-Identifier: BSD-3-Clause                                    *
+ ****************************************************************************/
+#ifndef ARBORX_SEGMENT_HPP
+#define ARBORX_SEGMENT_HPP
+
+#include <ArborX_DetailsAlgorithms.hpp>
+#include <ArborX_DetailsVector.hpp>
+#include <ArborX_GeometryTraits.hpp>
+
+#include <Kokkos_Clamp.hpp>
+#include <Kokkos_Macros.hpp>
+
+namespace ArborX::Experimental
+{
+template <int DIM, typename Coordinate = float>
+struct Segment
+{
+  ArborX::Point<DIM, Coordinate> _start;
+  ArborX::Point<DIM, Coordinate> _end;
+};
+
+template <int DIM, typename Coordinate>
+#if KOKKOS_VERSION >= 40400
+KOKKOS_DEDUCTION_GUIDE
+#else
+KOKKOS_FUNCTION
+#endif
+    Segment(Point<DIM, Coordinate>, Point<DIM, Coordinate>)
+        -> Segment<DIM, Coordinate>;
+
+template <int N, typename T>
+#if KOKKOS_VERSION >= 40400
+KOKKOS_DEDUCTION_GUIDE
+#else
+KOKKOS_FUNCTION
+#endif
+Segment(T const (&)[N], T const (&)[N]) -> Segment<N, T>;
+
+} // namespace ArborX::Experimental
+
+template <int DIM, class Coordinate>
+struct ArborX::GeometryTraits::dimension<
+    ArborX::Experimental::Segment<DIM, Coordinate>>
+{
+  static constexpr int value = DIM;
+};
+template <int DIM, class Coordinate>
+struct ArborX::GeometryTraits::tag<
+    ArborX::Experimental::Segment<DIM, Coordinate>>
+{
+  using type = SegmentTag;
+};
+template <int DIM, class Coordinate>
+struct ArborX::GeometryTraits::coordinate_type<
+    ArborX::Experimental::Segment<DIM, Coordinate>>
+{
+  using type = Coordinate;
+};
+
+namespace ArborX::Details::Dispatch
+{
+
+using GeometryTraits::BoxTag;
+using GeometryTraits::PointTag;
+using GeometryTraits::SegmentTag;
+
+template <typename Box, typename Segment>
+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);
+  }
+};
+
+template <typename Point, typename Segment>
+struct distance<PointTag, SegmentTag, Point, Segment>
+{
+  KOKKOS_FUNCTION static auto apply(Point const &point, Segment const &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);
+
+    auto const dir = segment._end - segment._start;
+
+    // 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),
+                      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];
+
+    return Details::distance(point, projection);
+  }
+};
+
+template <typename Segment>
+struct centroid<SegmentTag, Segment>
+{
+  KOKKOS_FUNCTION static auto apply(Segment const &segment)
+  {
+    constexpr int DIM = GeometryTraits::dimension_v<Segment>;
+    using Coordinate = GeometryTraits::coordinate_type_t<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;
+    return point;
+  }
+};
+
+} // namespace ArborX::Details::Dispatch
+
+#endif

test/tstCompileOnlyGeometry.cpp

@@ -12,6 +12,7 @@
 #include <ArborX_Box.hpp>
 #include <ArborX_GeometryTraits.hpp>
 #include <ArborX_Point.hpp>
+#include <ArborX_Segment.hpp>
 #include <ArborX_Sphere.hpp>
 
 namespace ArborX::GeometryTraits
@@ -248,6 +249,18 @@ void test_point_ctad()
                                Point<3, double>>);
 }
 
+void test_segment_ctad()
+{
+  using ArborX::Point;
+  using ArborX::Experimental::Segment;
+
+  static_assert(std::is_same_v<decltype(Segment{{0.f, 2.f}, {2.f, 3.f}}),
+                               Segment<2, float>>);
+  static_assert(std::is_same_v<decltype(Segment{Point<2, double>{0.f, 2.f},
+                                                Point<2, double>{2.f, 3.f}}),
+                               Segment<2, double>>);
+}
+
 void test_box_ctad()
 {
   using ArborX::Box;

test/tstDetailsAlgorithms.cpp

@@ -13,6 +13,7 @@
 #include <ArborX_DetailsAlgorithms.hpp>
 #include <ArborX_DetailsKokkosExtArithmeticTraits.hpp>
 #include <ArborX_Point.hpp>
+#include <ArborX_Segment.hpp>
 #include <ArborX_Sphere.hpp>
 #include <ArborX_Tetrahedron.hpp>
 #include <ArborX_Triangle.hpp>
@@ -65,6 +66,27 @@ BOOST_AUTO_TEST_CASE(distance_point_sphere)
   BOOST_TEST(distance(Point{{1., 1., 1.}}, sphere) == std::sqrt(3.f) - 1.f);
 }
 
+BOOST_AUTO_TEST_CASE(distance_point_segment)
+{
+  using ArborX::Details::distance;
+
+  using ArborX::Point;
+  using ArborX::Experimental::Segment;
+
+  constexpr Segment segment0{{0.f, 0.f}, {0.f, 0.f}};
+  BOOST_TEST(distance(Point{0.f, 0.f}, segment0) == 0.f);
+  BOOST_TEST(distance(Point{1.f, 0.f}, segment0) == 1.f);
+  BOOST_TEST(distance(Point{-1.f, 0.f}, segment0) == 1.f);
+
+  constexpr Segment segment1{{0.f, 0.f}, {1.f, 1.f}};
+  BOOST_TEST(distance(Point{0.f, 0.f}, segment1) == 0.f);
+  BOOST_TEST(distance(Point{1.f, 1.f}, segment1) == 0.f);
+  BOOST_TEST(distance(Point{0.5f, 0.5f}, segment1) == 0.f);
+  BOOST_TEST(distance(Point{1.0f, 0.f}, segment1) == std::sqrt(0.5f));
+  BOOST_TEST(distance(Point{0.f, -1.f}, segment1) == 1.f);
+  BOOST_TEST(distance(Point{1.5f, 1.f}, segment1) == 0.5f);
+}
+
 BOOST_AUTO_TEST_CASE(distance_point_triangle)
 {
   using ArborX::Details::distance;
@@ -346,6 +368,14 @@ BOOST_AUTO_TEST_CASE(expand)
   BOOST_TEST(equals(box, Box{{-5, -2, 2}, {1, 4, 7}}));
   expand(box, Tetrahedron{{-3, -5, 2}, {2, 6, -1}, {3, 2, 3}, {5, 8, -3}});
   BOOST_TEST(equals(box, Box{{-5, -5, -3}, {5, 8, 7}}));
+
+  // expand box with segments
+  using ArborX::Experimental::Segment;
+  box = Box{};
+  expand(box, Segment{{0.f, 0.f, 0.f}, {1.f, 1.f, 1.f}});
+  BOOST_TEST(equals(box, Box{{0, 0, 0}, {1, 1, 1}}));
+  expand(box, Segment{{-1.f, 3.0f, 0.0f}, {2.f, 1.f, 0.f}});
+  BOOST_TEST(equals(box, Box{{-1, 0, 0}, {2, 3, 1}}));
 }
 
 BOOST_AUTO_TEST_CASE(convert)
@@ -362,7 +392,9 @@ BOOST_AUTO_TEST_CASE(convert)
 
 BOOST_AUTO_TEST_CASE(centroid)
 {
+  using ArborX::Details::equals;
   using ArborX::Details::returnCentroid;
+
   Box box{{{-10.0, 0.0, 10.0}}, {{0.0, 10.0, 20.0}}};
   auto center = returnCentroid(box);
   BOOST_TEST(center[0] == -5.0);
@@ -385,6 +417,11 @@ BOOST_AUTO_TEST_CASE(centroid)
   BOOST_TEST(tet_center[0] == -1);
   BOOST_TEST(tet_center[1] == 3);
   BOOST_TEST(tet_center[2] == 2);
+
+  using ArborX::Experimental::Segment;
+  Segment segment{{-1.f, -1.f}, {3.f, 3.f}};
+  auto seg_center = returnCentroid(segment);
+  BOOST_TEST(equals(seg_center, ArborX::Point{1.f, 1.f}));
 }
 
 BOOST_AUTO_TEST_CASE(is_valid)