From 92caa1aa6ef40a31512253504f547ae6633c2ce0 Mon Sep 17 00:00:00 2001 From: "Seth R. Johnson" Date: Sat, 16 Sep 2023 07:31:41 -0400 Subject: [PATCH] Fix ORANGE bounding box bumping to be consistent with tracking tolerances (#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 | 77 ++++++++++++++++++++------ src/orange/detail/UnitInserter.cc | 10 +++- test/Test.hh | 1 + test/orange/BoundingBoxUtils.test.cc | 81 +++++++++++++++++----------- 4 files changed, 121 insertions(+), 48 deletions(-) diff --git a/src/orange/BoundingBoxUtils.hh b/src/orange/BoundingBoxUtils.hh index ba8c1ad809..41d3b96844 100644 --- a/src/orange/BoundingBoxUtils.hh +++ b/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 const& a, BoundingBox 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 -inline BoundingBox calc_bumped(BoundingBox const& bbox) +template +class BoundingBoxBumper { - CELER_EXPECT(bbox); + public: + //!@{ + //! \name Type aliases + using result_type = BoundingBox; + using argument_type = BoundingBox; + //!@} - Array lower; - Array upper; + public: + //! Construct with default "soft equal" tolerances + BoundingBoxBumper() + : rel_{SoftEqual{}.rel()}, abs_{SoftEqual{}.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(bbox.lower()[ax]), - -numeric_limits::infinity()); - upper[ax] = std::nextafter(static_cast(bbox.upper()[ax]), - numeric_limits::infinity()); + CELER_EXPECT(rel_ > 0 && abs_ > 0); } - return BoundingBox::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 lower; + Array 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 + T bumped(U value) const + { + U bumped = value + S * celeritas::max(abs_, rel_ * std::fabs(value)); + return std::nextafter(static_cast(bumped), + S * numeric_limits::infinity()); + } +}; //---------------------------------------------------------------------------// // Calculate the bounding box of a transformed box diff --git a/src/orange/detail/UnitInserter.cc b/src/orange/detail/UnitInserter.cc index b6e2758df3..77b2cf8298 100644 --- a/src/orange/detail/UnitInserter.cc +++ b/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 calc_bumped{2 * orange_data_->scalars.bump_rel, + 2 * orange_data_->scalars.bump_abs}; + // Define volumes std::vector vol_records(inp.volumes.size()); std::vector> 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(inp.volumes[i].bbox)); + bboxes.push_back(calc_bumped(inp.volumes[i].bbox)); } else { diff --git a/test/Test.hh b/test/Test.hh index 3b80fa12d9..89dfb5757e 100644 --- a/test/Test.hh +++ b/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; diff --git a/test/orange/BoundingBoxUtils.test.cc b/test/orange/BoundingBoxUtils.test.cc index c5962d72b9..fa55c0359e 100644 --- a/test/orange/BoundingBoxUtils.test.cc +++ b/test/orange/BoundingBoxUtils.test.cc @@ -176,38 +176,57 @@ TEST_F(BoundingBoxUtilsTest, bbox_intersection) TEST_F(BoundingBoxUtilsTest, bumped) { - auto inf_double = std::numeric_limits::infinity(); - auto inf_float = std::numeric_limits::infinity(); + BoundingBox 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(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{-inf_double, 0, -100})); - EXPECT_TRUE( - is_inside(bumped, Array{0, long_number, inf_double})); + { + SCOPED_TRACE("default precision"); + BoundingBoxBumper 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 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 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 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)