include/orange/BoundingBoxUtils.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2022-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file orange/BoundingBoxUtils.hh
0007 //! \brief Host-only utilities for bounding boxes
0008 //---------------------------------------------------------------------------//
0009 #pragma once
0010
0011 #include <cmath>
0012 #include <iosfwd>
0013
0014 #include "corecel/cont/Range.hh"
0015 #include "corecel/math/Algorithms.hh"
0016 #include "corecel/math/SoftEqual.hh"
0017 #include "geocel/BoundingBox.hh"
0018
0019 #include "OrangeTypes.hh"
0020
0021 namespace celeritas
0022 {
0023 //---------------------------------------------------------------------------//
0024 class Translation;
0025 class Transformation;
0026
0027 //---------------------------------------------------------------------------//
0028 /*!
0029  * Check if a bounding box spans (-inf, inf) in every direction.
0030  */
0031 template<class T>
0032 inline bool is_infinite(BoundingBox<T> const& bbox)
0033 {
0034     auto all_equal = [](Array<T, 3> const& values, T rhs) {
0035         return all_of(
0036             values.begin(), values.end(), [rhs](T lhs) { return lhs == rhs; });
0037     };
0038     constexpr T inf = numeric_limits<T>::infinity();
0039
0040     return all_equal(bbox.lower(), -inf) && all_equal(bbox.upper(), inf);
0041 }
0042
0043 //---------------------------------------------------------------------------//
0044 /*!
0045  * Check if a bounding box has no infinities.
0046  *
0047  * \pre The bounding box cannot be null
0048  */
0049 template<class T>
0050 inline bool is_finite(BoundingBox<T> const& bbox)
0051 {
0052     CELER_EXPECT(bbox);
0053
0054     auto isinf = [](T value) { return std::isinf(value); };
0055     return !any_of(bbox.lower().begin(), bbox.lower().end(), isinf)
0056            && !any_of(bbox.upper().begin(), bbox.upper().end(), isinf);
0057 }
0058
0059 //---------------------------------------------------------------------------//
0060 /*!
0061  * Check if a bounding box has zero length in any direction.
0062  *
0063  * \pre The bounding box cannot be null
0064  */
0065 template<class T>
0066 inline bool is_degenerate(BoundingBox<T> const& bbox)
0067 {
0068     CELER_EXPECT(bbox);
0069
0070     auto axes = range(to_int(Axis::size_));
0071     return any_of(axes.begin(), axes.end(), [&bbox](int ax) {
0072         return bbox.lower()[ax] == bbox.upper()[ax];
0073     });
0074 }
0075
0076 //---------------------------------------------------------------------------//
0077 /*!
0078  * Calculate the center of a bounding box.
0079  *
0080  * \pre The bounding box cannot be null
0081  */
0082 template<class T>
0083 inline Array<T, 3> calc_center(BoundingBox<T> const& bbox)
0084 {
0085     CELER_EXPECT(bbox);
0086
0087     Array<T, 3> center;
0088     for (auto ax : range(to_int(Axis::size_)))
0089     {
0090         center[ax] = (bbox.lower()[ax] + bbox.upper()[ax]) / 2;
0091     }
0092
0093     return center;
0094 }
0095
0096 //---------------------------------------------------------------------------//
0097 /*!
0098  * Calculate the half widths of the bounding box.
0099  *
0100  * \pre The bounding box cannot be null
0101  */
0102 template<class T>
0103 inline Array<T, 3> calc_half_widths(BoundingBox<T> const& bbox)
0104 {
0105     CELER_EXPECT(bbox);
0106
0107     Array<T, 3> hw;
0108     for (auto ax : range(to_int(Axis::size_)))
0109     {
0110         hw[ax] = (bbox.upper()[ax] - bbox.lower()[ax]) / 2;
0111     }
0112
0113     return hw;
0114 }
0115
0116 //---------------------------------------------------------------------------//
0117 /*!
0118  * Calculate the surface area of a bounding box.
0119  *
0120  * \pre The bounding box cannot be null
0121  */
0122 template<class T>
0123 inline T calc_surface_area(BoundingBox<T> const& bbox)
0124 {
0125     CELER_EXPECT(bbox);
0126
0127     Array<T, 3> lengths;
0128
0129     for (auto ax : range(to_int(Axis::size_)))
0130     {
0131         lengths[ax] = bbox.upper()[ax] - bbox.lower()[ax];
0132     }
0133
0134     return 2
0135            * (lengths[to_int(Axis::x)] * lengths[to_int(Axis::y)]
0136               + lengths[to_int(Axis::x)] * lengths[to_int(Axis::z)]
0137               + lengths[to_int(Axis::y)] * lengths[to_int(Axis::z)]);
0138 }
0139
0140 //---------------------------------------------------------------------------//
0141 /*!
0142  * Calculate the volume of a bounding box.
0143  *
0144  * \pre The bounding box cannot be null
0145  */
0146 template<class T>
0147 inline T calc_volume(BoundingBox<T> const& bbox)
0148 {
0149     CELER_EXPECT(bbox);
0150
0151     T result{1};
0152
0153     for (auto ax : range(to_int(Axis::size_)))
0154     {
0155         result *= bbox.upper()[ax] - bbox.lower()[ax];
0156     }
0157
0158     return result;
0159 }
0160
0161 //---------------------------------------------------------------------------//
0162 /*!
0163  * Calculate the smallest bounding box enclosing two bounding boxes.
0164  */
0165 template<class T>
0166 inline constexpr BoundingBox<T>
0167 calc_union(BoundingBox<T> const& a, BoundingBox<T> const& b)
0168 {
0169     Array<T, 3> lower{};
0170     Array<T, 3> upper{};
0171
0172     for (auto ax : range(to_int(Axis::size_)))
0173     {
0174         lower[ax] = celeritas::min(a.lower()[ax], b.lower()[ax]);
0175         upper[ax] = celeritas::max(a.upper()[ax], b.upper()[ax]);
0176     }
0177
0178     return BoundingBox<T>::from_unchecked(lower, upper);
0179 }
0180
0181 //---------------------------------------------------------------------------//
0182 /*!
0183  * Calculate the intersection of two bounding boxes.
0184  *
0185  * If there is no intersection, the result will be a null bounding box.
0186  */
0187 template<class T>
0188 inline constexpr BoundingBox<T>
0189 calc_intersection(BoundingBox<T> const& a, BoundingBox<T> const& b)
0190 {
0191     Array<T, 3> lower{};
0192     Array<T, 3> upper{};
0193
0194     for (auto ax : range(to_int(Axis::size_)))
0195     {
0196         lower[ax] = celeritas::max(a.lower()[ax], b.lower()[ax]);
0197         upper[ax] = celeritas::min(a.upper()[ax], b.upper()[ax]);
0198     }
0199
0200     return BoundingBox<T>::from_unchecked(lower, upper);
0201 }
0202
0203 //---------------------------------------------------------------------------//
0204 /*!
0205  * Check if all points inside the small bbox are in the big bbox.
0206  *
0207  * All bounding boxes should enclose a "null" bounding box (there are no points
0208  * in the null box, so no points are outside the big box). The null bounding
0209  * box will enclose no real bounding boxes. Comparing two null bounding boxes
0210  * is unspecified (forbidden for now).
0211  */
0212 template<class T>
0213 inline bool encloses(BoundingBox<T> const& big, BoundingBox<T> const& small)
0214 {
0215     CELER_EXPECT(big || small);
0216
0217     auto axes = range(to_int(Axis::size_));
0218     return all_of(axes.begin(), axes.end(), [&big, &small](int ax) {
0219         return big.lower()[ax] <= small.lower()[ax]
0220                && big.upper()[ax] >= small.upper()[ax];
0221     });
0222 }
0223
0224 //---------------------------------------------------------------------------//
0225 /*!
0226  * Bump a bounding box outward and possibly convert to another type.
0227  * \tparam T destination type
0228  * \tparam U source type
0229  *
0230  * The upper and lower coordinates are bumped outward independently using the
0231  * relative and absolute tolerances. To ensure that the outward bump is
0232  * not truncated in the destination type, the "std::nextafter" function
0233  * advances to the next floating point representable number.
0234  */
0235 template<class T, class U = T>
0236 class BoundingBoxBumper
0237 {
0238   public:
0239     //!@{
0240     //! \name Type aliases
0241     using TolU = Tolerance<U>;
0242     using result_type = BoundingBox<T>;
0243     using argument_type = BoundingBox<U>;
0244     //!@}
0245
0246   public:
0247     //! Construct with default "soft equal" tolerances
0248     BoundingBoxBumper() : tol_{TolU::from_softequal()} {}
0249
0250     //! Construct with ORANGE tolerances
0251     explicit BoundingBoxBumper(TolU const& tol) : tol_{tol}
0252     {
0253         CELER_EXPECT(tol_);
0254     }
0255
0256     //! Return the expanded and converted bounding box
0257     result_type operator()(argument_type const& bbox)
0258     {
0259         Array<T, 3> lower;
0260         Array<T, 3> upper;
0261
0262         for (auto ax : range(to_int(Axis::size_)))
0263         {
0264             lower[ax] = this->bumped<-1>(bbox.lower()[ax]);
0265             upper[ax] = this->bumped<+1>(bbox.upper()[ax]);
0266         }
0267
0268         return result_type::from_unchecked(lower, upper);
0269     }
0270
0271   private:
0272     TolU tol_;
0273
0274     //! Calculate the bump distance given a point: see detail::BumpCalculator
0275     template<int S>
0276     T bumped(U value) const
0277     {
0278         U bumped = value
0279                    + S * celeritas::max(tol_.abs, tol_.rel * std::fabs(value));
0280         return std::nextafter(static_cast<T>(bumped),
0281                               S * numeric_limits<T>::infinity());
0282     }
0283 };
0284
0285 //---------------------------------------------------------------------------//
0286 // Calculate the bounding box of a transformed box
0287 BBox calc_transform(Translation const& tr, BBox const& a);
0288
0289 BBox calc_transform(Transformation const& tr, BBox const& a);
0290
0291 //---------------------------------------------------------------------------//
0292 template<class T>
0293 std::ostream& operator<<(std::ostream&, BoundingBox<T> const& bbox);
0294
0295 //---------------------------------------------------------------------------//
0296 }  // namespace celeritas