Add distance-to-bbox function (#1446)

src/geocel/Types.hh

@@ -72,10 +72,11 @@ enum class Axis
  * Here, lo/hi correspond to left/right, back/front, bottom/top. It's used for
  * the two points in a bounding box.
  */
-enum class Bound : bool
+enum class Bound : unsigned char
 {
     lo,
-    hi
+    hi,
+    size_
 };
 
 //---------------------------------------------------------------------------//

src/orange/BoundingBoxUtils.hh

@@ -221,6 +221,72 @@ inline bool encloses(BoundingBox<T> const& big, BoundingBox<T> const& small)
     });
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate the distance to the inside of the bbox from a pos and dir.
+ *
+ * The supplied position is expected to be outside of the bbox. If there is no
+ * intersection, the result will be inf.
+ */
+template<class T, class U>
+inline U calc_dist_to_inside(BoundingBox<T> const& bbox,
+                             Array<U, 3> const& pos,
+                             Array<U, 3> const& dir)
+{
+    CELER_EXPECT(!is_inside(bbox, pos));
+
+    // Test if an intersection is outside the bbox for a given axis
+    auto out_of_bounds = [&bbox](U intersect, int ax) {
+        return !(intersect >= bbox.lower()[ax]
+                 && intersect <= bbox.upper()[ax]);
+    };
+
+    // Check that the intersection point occurs within the region
+    // bounded by the planes of the other two axes
+    auto in_bounds = [&](int ax, U dist) {
+        for (auto other_ax : range(to_int(Axis::size_)))
+        {
+            if (other_ax == ax)
+                continue;
+
+            auto intersect = pos[other_ax] + dist * dir[other_ax];
+            if (out_of_bounds(intersect, other_ax))
+                return false;
+        }
+        return true;
+    };
+
+    // Loop over all 6 planes to find the minimum intersection
+    U min_dist = numeric_limits<U>::infinity();
+    for (auto bound : range(to_int(Bound::size_)))
+    {
+        for (auto ax : range(to_int(Axis::size_)))
+        {
+            if (dir[ax] == 0)
+            {
+                // Short circut if there is not movement in this dir
+                continue;
+            }
+
+            U dist = (bbox.point(static_cast<Bound>(bound))[ax] - pos[ax])
+                     / dir[ax];
+
+            if (dist < 0)
+            {
+                // Short circut if the plane is behind us
+                continue;
+            }
+
+            if (in_bounds(ax, dist))
+            {
+                min_dist = celeritas::min(min_dist, dist);
+            }
+        }
+    }
+
+    return min_dist;
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Bump a bounding box outward and possibly convert to another type.

test/orange/BoundingBoxUtils.test.cc

@@ -175,6 +175,39 @@ TEST_F(BoundingBoxUtilsTest, bbox_intersection)
     }
 }
 
+TEST_F(BoundingBoxUtilsTest, bbox_dist_to_inside)
+{
+    using Real3 = Array<double, 3>;
+
+    auto bbox = BBox{{0., 0., 0.}, {1, 1, 1}};
+
+    // Basic case
+    Real3 pos{1.1, 0.5, 0.5};
+    Real3 dir{-1, 0, 0};
+    EXPECT_SOFT_EQ(0.1, calc_dist_to_inside(bbox, pos, dir));
+
+    // Coming in from an angle
+    dir = Real3{-std::sqrt(2) / 2, -std::sqrt(2) / 2, 0};
+    EXPECT_SOFT_EQ(0.1 * std::sqrt(2), calc_dist_to_inside(bbox, pos, dir));
+
+    // First intersection point occurs outside box, but second intersection
+    // point is valid
+    pos = Real3{3, 2.5, 0.5};
+    EXPECT_SOFT_EQ(2 * std::sqrt(2), calc_dist_to_inside(bbox, pos, dir));
+
+    // No intersection
+    dir = Real3{0, -1, 0};
+    EXPECT_EQ(numeric_limits<double>::infinity(),
+              calc_dist_to_inside(bbox, pos, dir));
+
+    // Already inside
+    if (CELERITAS_DEBUG)
+    {
+        pos = Real3{0.5, 0.6, 0.7};
+        EXPECT_THROW(calc_dist_to_inside(bbox, pos, dir), DebugError);
+    }
+}
+
 TEST_F(BoundingBoxUtilsTest, bbox_encloses)
 {
     EXPECT_TRUE(encloses(BBox::from_infinite(), BBox{{-1, -1, -1}, {1, 1, 1}}));