Fix ORANGE bounding box bumping to be consistent with tracking tolera…

…nces (#940)

* Implement bounding box bumper

* Use relative points instead of width due to infinities

* Refactor bumper to ensure encloses post-bump point

1e-8 is too small to change the single-precision representation. The
new implementation combines the relative + absolute bump to the
full-precision point, then uses `std::nextafter` to ensure that the
truncated point is outside that value.

src/orange/BoundingBoxUtils.hh

@@ -11,6 +11,7 @@
 
 #include "corecel/cont/Range.hh"
 #include "corecel/math/Algorithms.hh"
+#include "corecel/math/SoftEqual.hh"
 #include "orange/BoundingBox.hh"
 #include "orange/OrangeTypes.hh"
 
@@ -181,30 +182,74 @@ calc_intersection(BoundingBox<T> const& a, BoundingBox<T> const& b)
 
 //---------------------------------------------------------------------------//
 /*!
- * Convert bbox with U type values to a bumped bbox with T type values.
+ * Bump a bounding box outward and possibly convert to another type.
+ * \tparam T destination type
+ * \tparam U source type
  *
- * Each U lower value is bumped to the greatest T value less than the U value.
- * Each U upper value is bumped to the lowest T value greater than the U value.
- * Infinite values are unchanged.
+ * The upper and lower coordinates are bumped outward independently using the
+ * relative and absolute tolerances. To ensure that the outward bump is
+ * not truncated in the destination type, the "std::nextafter" function
+ * advances to the next floating point representable number.
  */
-template<class T, class U>
-inline BoundingBox<T> calc_bumped(BoundingBox<U> const& bbox)
+template<class T, class U = real_type>
+class BoundingBoxBumper
 {
-    CELER_EXPECT(bbox);
+  public:
+    //!@{
+    //! \name Type aliases
+    using result_type = BoundingBox<T>;
+    using argument_type = BoundingBox<U>;
+    //!@}
 
-    Array<T, 3> lower;
-    Array<T, 3> upper;
+  public:
+    //! Construct with default "soft equal" tolerances
+    BoundingBoxBumper()
+        : rel_{SoftEqual<U>{}.rel()}, abs_{SoftEqual<U>{}.abs()}
+    {
+    }
 
-    for (auto ax : range(to_int(Axis::size_)))
+    //! Construct with a single bump tolerance used for both relative and abs
+    explicit BoundingBoxBumper(U tol) : rel_{tol}, abs_{tol}
     {
-        lower[ax] = std::nextafter(static_cast<T>(bbox.lower()[ax]),
-                                   -numeric_limits<T>::infinity());
-        upper[ax] = std::nextafter(static_cast<T>(bbox.upper()[ax]),
-                                   numeric_limits<T>::infinity());
+        CELER_EXPECT(rel_ > 0 && abs_ > 0);
     }
 
-    return BoundingBox<T>::from_unchecked(lower, upper);
-}
+    //! Construct with relative and absolute bump tolerances
+    BoundingBoxBumper(U rel, U abs) : rel_{rel}, abs_{abs}
+    {
+        CELER_EXPECT(rel_ > 0 && abs_ > 0);
+    }
+
+    //! Return the expanded and converted bounding box
+    result_type operator()(argument_type const& bbox)
+    {
+        CELER_EXPECT(bbox);
+
+        Array<T, 3> lower;
+        Array<T, 3> upper;
+
+        for (auto ax : range(to_int(Axis::size_)))
+        {
+            lower[ax] = this->bumped<-1>(bbox.lower()[ax]);
+            upper[ax] = this->bumped<+1>(bbox.upper()[ax]);
+        }
+
+        return result_type::from_unchecked(lower, upper);
+    }
+
+  private:
+    U rel_;
+    U abs_;
+
+    //! Calculate the bump distance given a point: see detail::BumpCalculator
+    template<int S>
+    T bumped(U value) const
+    {
+        U bumped = value + S * celeritas::max(abs_, rel_ * std::fabs(value));
+        return std::nextafter(static_cast<T>(bumped),
+                              S * numeric_limits<T>::infinity());
+    }
+};
 
 //---------------------------------------------------------------------------//
 // Calculate the bounding box of a transformed box

src/orange/detail/UnitInserter.cc

@@ -131,6 +131,8 @@ UnitInserter::UnitInserter(Data* orange_data)
     , insert_transform_{&orange_data_->transforms, &orange_data_->reals}
 {
     CELER_EXPECT(orange_data);
+    CELER_EXPECT(orange_data->scalars.bump_rel > 0);
+    CELER_EXPECT(orange_data->scalars.bump_abs > 0);
 
     // Initialize scalars
     orange_data_->scalars.max_faces = 1;
@@ -148,6 +150,12 @@ SimpleUnitId UnitInserter::operator()(UnitInput const& inp)
     // Insert surfaces
     unit.surfaces = this->insert_surfaces(inp.surfaces);
 
+    // Bounding box bumper and converter: expand to twice the potential bump
+    // distance from a boundary so that the bbox will enclose the point even
+    // after a potential bump
+    BoundingBoxBumper<float> calc_bumped{2 * orange_data_->scalars.bump_rel,
+                                         2 * orange_data_->scalars.bump_abs};
+
     // Define volumes
     std::vector<VolumeRecord> vol_records(inp.volumes.size());
     std::vector<std::set<LocalVolumeId>> connectivity(inp.surfaces.size());
@@ -160,7 +168,7 @@ SimpleUnitId UnitInserter::operator()(UnitInput const& inp)
         // Store the bbox or an infinite bbox placeholder
         if (inp.volumes[i].bbox)
         {
-            bboxes.push_back(calc_bumped<fast_real_type>(inp.volumes[i].bbox));
+            bboxes.push_back(calc_bumped(inp.volumes[i].bbox));
         }
         else
         {

test/Test.hh

@@ -48,6 +48,7 @@ class Test : public ::testing::Test
 
     // Define "inf" value for subclass testing
     static constexpr double inf = HUGE_VAL;
+    static constexpr float inff = HUGE_VALF;
 
   private:
     int filename_counter_ = 0;

test/orange/BoundingBoxUtils.test.cc

@@ -176,38 +176,57 @@ TEST_F(BoundingBoxUtilsTest, bbox_intersection)
 
 TEST_F(BoundingBoxUtilsTest, bumped)
 {
-    auto inf_double = std::numeric_limits<double>::infinity();
-    auto inf_float = std::numeric_limits<float>::infinity();
+    BoundingBox<double> const ref{{-inf, 0, -100}, {0, 0.11223344556677, inf}};
 
-    double long_number = 0.11223344556677;
-
-    BBox unbumped = {{-inf_double, 0, -100}, {0, long_number, inf_double}};
-
-    auto bumped = calc_bumped<float>(unbumped);
-
-    // Test lower corner
-    EXPECT_EQ(-inf_float, bumped.lower()[0]);
-
-    EXPECT_SOFT_EQ(0, bumped.lower()[1]);
-    EXPECT_TRUE(bumped.lower()[1] < 0);
-
-    EXPECT_SOFT_EQ(-100, bumped.lower()[2]);
-    EXPECT_TRUE(bumped.lower()[2] < -100);
-
-    // Test upper corner
-
-    EXPECT_SOFT_EQ(0, bumped.upper()[0]);
-    EXPECT_TRUE(bumped.upper()[0] > 0);
-
-    EXPECT_SOFT_EQ(long_number, bumped.upper()[1]);
-    EXPECT_TRUE(bumped.upper()[1] > long_number);
-
-    EXPECT_EQ(inf_float, bumped.upper()[2]);
-
-    // Test the bounds are inside
-    EXPECT_TRUE(is_inside(bumped, Array<double, 3>{-inf_double, 0, -100}));
-    EXPECT_TRUE(
-        is_inside(bumped, Array<double, 3>{0, long_number, inf_double}));
+    {
+        SCOPED_TRACE("default precision");
+        BoundingBoxBumper<float, double> calc_bumped{};
+        auto bumped = calc_bumped(ref);
+        static float const expected_lower[] = {-inff, -1e-14f, -100.0f};
+        static float const expected_upper[] = {1e-14f, 0.1122335f, inff};
+        EXPECT_VEC_SOFT_EQ(expected_lower, bumped.lower());
+        EXPECT_VEC_SOFT_EQ(expected_upper, bumped.upper());
+
+        EXPECT_TRUE(is_inside(bumped, ref.lower()));
+        EXPECT_TRUE(is_inside(bumped, ref.upper()));
+    }
+    {
+        SCOPED_TRACE("double precise");
+        BoundingBoxBumper<double> calc_bumped{1e-10};
+        auto bumped = calc_bumped(ref);
+        static double const expected_lower[] = {-inf, -1e-10, -100.00000001};
+        static double const expected_upper[] = {1e-10, 0.11223344566677, inf};
+        EXPECT_VEC_SOFT_EQ(expected_lower, bumped.lower());
+        EXPECT_VEC_SOFT_EQ(expected_upper, bumped.upper());
+
+        EXPECT_TRUE(is_inside(bumped, ref.lower() - 1e-11));
+        EXPECT_TRUE(is_inside(bumped, ref.upper() + 1e-11));
+    }
+    {
+        SCOPED_TRACE("float loose");
+        BoundingBoxBumper<float, double> calc_bumped{1e-3, 1e-5};
+        auto bumped = calc_bumped(ref);
+        static float const expected_lower[] = {-inff, -1e-05f, -100.1f};
+        static float const expected_upper[] = {1e-05f, 0.1123457f, inff};
+        EXPECT_VEC_SOFT_EQ(expected_lower, bumped.lower());
+        EXPECT_VEC_SOFT_EQ(expected_upper, bumped.upper());
+
+        EXPECT_TRUE(is_inside(bumped, ref.lower() - 1e-6));
+        EXPECT_TRUE(is_inside(bumped, ref.upper() + 1e-6));
+    }
+    {
+        SCOPED_TRACE("float orange");
+        BBox const ref{{-2, -6, -1}, {8, 4, 2}};
+        BoundingBoxBumper<float, double> calc_bumped{2e-8, 2e-8};
+        auto bumped = calc_bumped(ref);
+        static float const expected_lower[] = {-2.f, -6.f, -1.f};
+        static float const expected_upper[] = {8.000001f, 4.f, 2.f};
+        EXPECT_VEC_SOFT_EQ(expected_lower, bumped.lower());
+        EXPECT_VEC_SOFT_EQ(expected_upper, bumped.upper());
+
+        EXPECT_TRUE(is_inside(bumped, ref.lower() - 1e-8));
+        EXPECT_TRUE(is_inside(bumped, ref.upper() + 1e-8));
+    }
 }
 
 TEST_F(BoundingBoxUtilsTest, bbox_translate)