From a425c6882a723b8ce080c8f8cd96a3edd644c920 Mon Sep 17 00:00:00 2001 From: Pasukhin Dmitry Date: Mon, 22 Dec 2025 16:51:46 +0000 Subject: [PATCH] Foundation Classes - Improve BVH Box and Rays (#882) - Refactored Add method in BVH_Box to utilize in-place component-wise operations, reducing temporary vector creation. - Introduced precomputed reciprocal direction in BVH_Ray for faster ray-box intersection tests. - Updated RayBoxIntersection methods to leverage the new reciprocal direction for optimal performance. - Enhanced documentation for clarity on new functionalities and optimizations. --- src/FoundationClasses/TKMath/BVH/BVH_Box.hxx | 331 ++++++------------ src/FoundationClasses/TKMath/BVH/BVH_Ray.hxx | 47 ++- .../TKMath/BVH/BVH_Tools.hxx | 54 ++- 3 files changed, 199 insertions(+), 233 deletions(-) diff --git a/src/FoundationClasses/TKMath/BVH/BVH_Box.hxx b/src/FoundationClasses/TKMath/BVH/BVH_Box.hxx index 8f10f27e28..1aad593689 100644 --- a/src/FoundationClasses/TKMath/BVH/BVH_Box.hxx +++ b/src/FoundationClasses/TKMath/BVH/BVH_Box.hxx @@ -120,51 +120,48 @@ class BVH_Box : public BVH_BaseBox public: typedef typename BVH::VectorType::Type BVH_VecNt; +private: + //! Returns the minimum point sentinel value for invalid box. + static constexpr T minSentinel() noexcept { return (std::numeric_limits::max)(); } + + //! Returns the maximum point sentinel value for invalid box. + static constexpr T maxSentinel() noexcept { return (std::numeric_limits::lowest)(); } + public: //! Creates uninitialized bounding box. constexpr BVH_Box() noexcept - : myIsInited(Standard_False) + : myMinPoint(BVH_VecNt(minSentinel())), + myMaxPoint(BVH_VecNt(maxSentinel())) { } //! Creates bounding box of given point. constexpr BVH_Box(const BVH_VecNt& thePoint) noexcept : myMinPoint(thePoint), - myMaxPoint(thePoint), - myIsInited(Standard_True) + myMaxPoint(thePoint) { } //! Creates bounding box from corner points. constexpr BVH_Box(const BVH_VecNt& theMinPoint, const BVH_VecNt& theMaxPoint) noexcept : myMinPoint(theMinPoint), - myMaxPoint(theMaxPoint), - myIsInited(Standard_True) + myMaxPoint(theMaxPoint) { } public: //! Clears bounding box. - constexpr void Clear() noexcept { myIsInited = Standard_False; } + constexpr void Clear() noexcept + { + myMinPoint = BVH_VecNt(minSentinel()); + myMaxPoint = BVH_VecNt(maxSentinel()); + } //! Is bounding box valid? - constexpr Standard_Boolean IsValid() const noexcept { return myIsInited; } + constexpr Standard_Boolean IsValid() const noexcept { return myMinPoint[0] <= myMaxPoint[0]; } //! Appends new point to the bounding box. - void Add(const BVH_VecNt& thePoint) - { - if (!myIsInited) - { - myMinPoint = thePoint; - myMaxPoint = thePoint; - myIsInited = Standard_True; - } - else - { - myMinPoint = myMinPoint.cwiseMin(thePoint); - myMaxPoint = myMaxPoint.cwiseMax(thePoint); - } - } + void Add(const BVH_VecNt& thePoint); //! Combines bounding box with another one. void Combine(const BVH_Box& theBox); @@ -198,13 +195,14 @@ public: void DumpJson(Standard_OStream& theOStream, Standard_Integer theDepth = -1) const { (void)theDepth; - OCCT_DUMP_FIELD_VALUE_NUMERICAL(theOStream, myIsInited) + const Standard_Integer anIsValid = IsValid() ? 1 : 0; + OCCT_DUMP_FIELD_VALUE_NUMERICAL(theOStream, anIsValid) constexpr int n = (N < 3) ? N : 3; if constexpr (n == 1) { OCCT_DUMP_FIELD_VALUE_NUMERICAL(theOStream, myMinPoint[0]) - OCCT_DUMP_FIELD_VALUE_NUMERICAL(theOStream, myMinPoint[0]) + OCCT_DUMP_FIELD_VALUE_NUMERICAL(theOStream, myMaxPoint[0]) } else if constexpr (n == 2) { @@ -233,18 +231,26 @@ public: { Standard_Integer aPos = theStreamPos; - Standard_Integer anIsInited = 0; + Standard_Integer anIsValid = 0; TCollection_AsciiString aStreamStr = Standard_Dump::Text(theSStream); - OCCT_INIT_FIELD_VALUE_INTEGER(aStreamStr, aPos, anIsInited); - myIsInited = anIsInited != 0; + OCCT_INIT_FIELD_VALUE_INTEGER(aStreamStr, aPos, anIsValid); + + if (anIsValid == 0) + { + Clear(); // Set to invalid state using sentinel values + theStreamPos = aPos; + return Standard_True; + } constexpr int n = (N < 3) ? N : 3; if constexpr (n == 1) { - Standard_Real aValue; - OCCT_INIT_FIELD_VALUE_REAL(aStreamStr, aPos, aValue); - myMinPoint[0] = (T)aValue; + Standard_Real aMinValue, aMaxValue; + OCCT_INIT_FIELD_VALUE_REAL(aStreamStr, aPos, aMinValue); + OCCT_INIT_FIELD_VALUE_REAL(aStreamStr, aPos, aMaxValue); + myMinPoint[0] = (T)aMinValue; + myMaxPoint[0] = (T)aMaxValue; } else if constexpr (n == 2) { @@ -271,6 +277,14 @@ public: myMaxPoint[2] = (T)aValue3; } + // For N > 3, initialize remaining dimensions to unbounded range + // so they don't affect intersection checks + for (int i = n; i < N; ++i) + { + myMinPoint[i] = (std::numeric_limits::lowest)(); + myMaxPoint[i] = (std::numeric_limits::max)(); + } + theStreamPos = aPos; return Standard_True; } @@ -291,19 +305,9 @@ public: if (!IsValid()) return Standard_True; - if constexpr (N >= 1) + for (int i = 0; i < N; ++i) { - if (myMinPoint[0] > theMaxPoint[0] || myMaxPoint[0] < theMinPoint[0]) - return Standard_True; - } - if constexpr (N >= 2) - { - if (myMinPoint[1] > theMaxPoint[1] || myMaxPoint[1] < theMinPoint[1]) - return Standard_True; - } - if constexpr (N >= 3) - { - if (myMinPoint[2] > theMaxPoint[2] || myMaxPoint[2] < theMinPoint[2]) + if (myMinPoint[i] > theMaxPoint[i] || myMaxPoint[i] < theMinPoint[i]) return Standard_True; } return Standard_False; @@ -330,27 +334,12 @@ public: return Standard_False; Standard_Boolean isInside = Standard_True; - - if constexpr (N >= 1) + for (int i = 0; i < N; ++i) { - hasOverlap = (myMinPoint[0] <= theMaxPoint[0] && myMaxPoint[0] >= theMinPoint[0]); + hasOverlap = (myMinPoint[i] <= theMaxPoint[i] && myMaxPoint[i] >= theMinPoint[i]); if (!hasOverlap) return Standard_False; - isInside = isInside && (myMinPoint[0] <= theMinPoint[0] && myMaxPoint[0] >= theMaxPoint[0]); - } - if constexpr (N >= 2) - { - hasOverlap = (myMinPoint[1] <= theMaxPoint[1] && myMaxPoint[1] >= theMinPoint[1]); - if (!hasOverlap) - return Standard_False; - isInside = isInside && (myMinPoint[1] <= theMinPoint[1] && myMaxPoint[1] >= theMaxPoint[1]); - } - if constexpr (N >= 3) - { - hasOverlap = (myMinPoint[2] <= theMaxPoint[2] && myMaxPoint[2] >= theMinPoint[2]); - if (!hasOverlap) - return Standard_False; - isInside = isInside && (myMinPoint[2] <= theMinPoint[2] && myMaxPoint[2] >= theMaxPoint[2]); + isInside = isInside && (myMinPoint[i] <= theMinPoint[i] && myMaxPoint[i] >= theMaxPoint[i]); } return isInside; } @@ -361,28 +350,17 @@ public: if (!IsValid()) return Standard_True; - if constexpr (N >= 1) + for (int i = 0; i < N; ++i) { - if (thePoint[0] < myMinPoint[0] || thePoint[0] > myMaxPoint[0]) - return Standard_True; - } - if constexpr (N >= 2) - { - if (thePoint[1] < myMinPoint[1] || thePoint[1] > myMaxPoint[1]) - return Standard_True; - } - if constexpr (N >= 3) - { - if (thePoint[2] < myMinPoint[2] || thePoint[2] > myMaxPoint[2]) + if (thePoint[i] < myMinPoint[i] || thePoint[i] > myMaxPoint[i]) return Standard_True; } return Standard_False; } protected: - BVH_VecNt myMinPoint; //!< Minimum point of bounding box - BVH_VecNt myMaxPoint; //!< Maximum point of bounding box - Standard_Boolean myIsInited; //!< Is bounding box initialized? + BVH_VecNt myMinPoint; //!< Minimum point of bounding box (max when invalid) + BVH_VecNt myMaxPoint; //!< Maximum point of bounding box (lowest when invalid) }; namespace BVH @@ -393,131 +371,56 @@ namespace BVH template struct CenterAxis { - // Not implemented -}; - -template -struct CenterAxis -{ - static inline T Center(const BVH_Box& theBox, const Standard_Integer theAxis) + //! Returns the center of the box along the specified axis using array access. + static inline T Center(const BVH_Box& theBox, const Standard_Integer theAxis) { - if (theAxis == 0) - { - return (theBox.CornerMin().x() + theBox.CornerMax().x()) * static_cast(0.5); - } - else if (theAxis == 1) - { - return (theBox.CornerMin().y() + theBox.CornerMax().y()) * static_cast(0.5); - } - return static_cast(0.0); - } -}; - -template -struct CenterAxis -{ - static inline T Center(const BVH_Box& theBox, const Standard_Integer theAxis) - { - if (theAxis == 0) - { - return (theBox.CornerMin().x() + theBox.CornerMax().x()) * static_cast(0.5); - } - else if (theAxis == 1) - { - return (theBox.CornerMin().y() + theBox.CornerMax().y()) * static_cast(0.5); - } - else if (theAxis == 2) - { - return (theBox.CornerMin().z() + theBox.CornerMax().z()) * static_cast(0.5); - } - return static_cast(0.0); - } -}; - -template -struct CenterAxis -{ - static inline T Center(const BVH_Box& theBox, const Standard_Integer theAxis) - { - if (theAxis == 0) - { - return (theBox.CornerMin().x() + theBox.CornerMax().x()) * static_cast(0.5); - } - else if (theAxis == 1) - { - return (theBox.CornerMin().y() + theBox.CornerMax().y()) * static_cast(0.5); - } - else if (theAxis == 2) - { - return (theBox.CornerMin().z() + theBox.CornerMax().z()) * static_cast(0.5); - } - return static_cast(0.0); + return (theBox.CornerMin()[theAxis] + theBox.CornerMax()[theAxis]) * static_cast(0.5); } }; //! Tool class for calculating surface area of the box. +//! For N=1, computes length (degenerate case). +//! For N=2, computes area (or perimeter for degenerate boxes). +//! For N>=3, computes 3D surface area using X, Y, Z components only. +//! The W component (4th dimension) is intentionally ignored as BVH surface area +//! heuristic (SAH) operates in 3D geometric space regardless of additional dimensions. //! \tparam T Numeric data type //! \tparam N Vector dimension template struct SurfaceCalculator { - // Not implemented -}; - -template -struct SurfaceCalculator -{ - static inline T Area(const typename BVH_Box::BVH_VecNt& theSize) + static inline T Area(const typename BVH_Box::BVH_VecNt& theSize) { - const T anArea = std::abs(theSize.x() * theSize.y()); - - if (anArea < std::numeric_limits::epsilon()) + if constexpr (N == 1) { - return std::abs(theSize.x()) + std::abs(theSize.y()); + // For 1D, return the length + return std::abs(theSize[0]); + } + else if constexpr (N == 2) + { + const T anArea = std::abs(theSize.x() * theSize.y()); + if (anArea < std::numeric_limits::epsilon()) + { + return std::abs(theSize.x()) + std::abs(theSize.y()); + } + return anArea; + } + else + { + // For N >= 3, compute standard 3D surface area. + const T anArea = (std::abs(theSize.x() * theSize.y()) + std::abs(theSize.x() * theSize.z()) + + std::abs(theSize.z() * theSize.y())) + * static_cast(2.0); + if (anArea < std::numeric_limits::epsilon()) + { + return std::abs(theSize.x()) + std::abs(theSize.y()) + std::abs(theSize.z()); + } + return anArea; } - - return anArea; } }; -template -struct SurfaceCalculator -{ - static inline T Area(const typename BVH_Box::BVH_VecNt& theSize) - { - const T anArea = (std::abs(theSize.x() * theSize.y()) + std::abs(theSize.x() * theSize.z()) - + std::abs(theSize.z() * theSize.y())) - * static_cast(2.0); - - if (anArea < std::numeric_limits::epsilon()) - { - return std::abs(theSize.x()) + std::abs(theSize.y()) + std::abs(theSize.z()); - } - - return anArea; - } -}; - -template -struct SurfaceCalculator -{ - static inline T Area(const typename BVH_Box::BVH_VecNt& theSize) - { - const T anArea = (std::abs(theSize.x() * theSize.y()) + std::abs(theSize.x() * theSize.z()) - + std::abs(theSize.z() * theSize.y())) - * static_cast(2.0); - - if (anArea < std::numeric_limits::epsilon()) - { - return std::abs(theSize.x()) + std::abs(theSize.y()) + std::abs(theSize.z()); - } - - return anArea; - } -}; - -//! Tool class for calculate component-wise vector minimum -//! and maximum (optimized version). +//! Tool class for computing component-wise vector minimum and maximum. //! \tparam T Numeric data type //! \tparam N Vector dimension template @@ -525,58 +428,44 @@ struct BoxMinMax { typedef typename BVH::VectorType::Type BVH_VecNt; + //! Computes component-wise minimum in-place. static inline void CwiseMin(BVH_VecNt& theVec1, const BVH_VecNt& theVec2) { - theVec1.x() = (std::min)(theVec1.x(), theVec2.x()); - theVec1.y() = (std::min)(theVec1.y(), theVec2.y()); - theVec1.z() = (std::min)(theVec1.z(), theVec2.z()); + for (int i = 0; i < N; ++i) + { + theVec1[i] = (std::min)(theVec1[i], theVec2[i]); + } } + //! Computes component-wise maximum in-place. static inline void CwiseMax(BVH_VecNt& theVec1, const BVH_VecNt& theVec2) { - theVec1.x() = (std::max)(theVec1.x(), theVec2.x()); - theVec1.y() = (std::max)(theVec1.y(), theVec2.y()); - theVec1.z() = (std::max)(theVec1.z(), theVec2.z()); - } -}; - -template -struct BoxMinMax -{ - typedef typename BVH::VectorType::Type BVH_VecNt; - - static inline void CwiseMin(BVH_VecNt& theVec1, const BVH_VecNt& theVec2) - { - theVec1.x() = (std::min)(theVec1.x(), theVec2.x()); - theVec1.y() = (std::min)(theVec1.y(), theVec2.y()); - } - - static inline void CwiseMax(BVH_VecNt& theVec1, const BVH_VecNt& theVec2) - { - theVec1.x() = (std::max)(theVec1.x(), theVec2.x()); - theVec1.y() = (std::max)(theVec1.y(), theVec2.y()); + for (int i = 0; i < N; ++i) + { + theVec1[i] = (std::max)(theVec1[i], theVec2[i]); + } } }; } // namespace BVH //================================================================================================= +template +void BVH_Box::Add(const BVH_VecNt& thePoint) +{ + BVH::BoxMinMax::CwiseMin(myMinPoint, thePoint); + BVH::BoxMinMax::CwiseMax(myMaxPoint, thePoint); +} + +//================================================================================================= + template void BVH_Box::Combine(const BVH_Box& theBox) { - if (theBox.myIsInited) + if (theBox.IsValid()) { - if (!myIsInited) - { - myMinPoint = theBox.myMinPoint; - myMaxPoint = theBox.myMaxPoint; - myIsInited = Standard_True; - } - else - { - BVH::BoxMinMax::CwiseMin(myMinPoint, theBox.myMinPoint); - BVH::BoxMinMax::CwiseMax(myMaxPoint, theBox.myMaxPoint); - } + BVH::BoxMinMax::CwiseMin(myMinPoint, theBox.myMinPoint); + BVH::BoxMinMax::CwiseMax(myMaxPoint, theBox.myMaxPoint); } } @@ -585,8 +474,8 @@ void BVH_Box::Combine(const BVH_Box& theBox) template T BVH_Box::Area() const { - return !myIsInited ? static_cast(0.0) - : BVH::SurfaceCalculator::Area(myMaxPoint - myMinPoint); + return !IsValid() ? static_cast(0.0) + : BVH::SurfaceCalculator::Area(myMaxPoint - myMinPoint); } //================================================================================================= diff --git a/src/FoundationClasses/TKMath/BVH/BVH_Ray.hxx b/src/FoundationClasses/TKMath/BVH/BVH_Ray.hxx index ce418f55b7..d65181fd9a 100644 --- a/src/FoundationClasses/TKMath/BVH/BVH_Ray.hxx +++ b/src/FoundationClasses/TKMath/BVH/BVH_Ray.hxx @@ -18,6 +18,8 @@ #include +#include + //! Describes a ray based on BVH vectors. template class BVH_Ray @@ -26,23 +28,60 @@ public: typedef typename BVH::VectorType::Type BVH_VecNt; public: - BVH_VecNt Origin; //!< Ray origin point - BVH_VecNt Direct; //!< Ray direction vector + BVH_VecNt Origin; //!< Ray origin point + BVH_VecNt Direct; //!< Ray direction vector + BVH_VecNt InvDirect; //!< Reciprocal of direction (1/Direct) public: //! Creates ray with given origin and direction. constexpr BVH_Ray(const BVH_VecNt& theOrigin, const BVH_VecNt& theDirect) noexcept : Origin(theOrigin), - Direct(theDirect) + Direct(theDirect), + InvDirect(computeInvDirect(theDirect)) { } //! Default constructor (creates invalid ray at origin). constexpr BVH_Ray() noexcept : Origin(BVH_VecNt()), - Direct(BVH_VecNt()) + Direct(BVH_VecNt()), + InvDirect(BVH_VecNt()) { } + +private: + //! Computes reciprocal of direction component (returns infinity for zero). + static constexpr T invComponent(T theDir) noexcept + { + return (theDir != T(0)) ? (T(1) / theDir) : std::numeric_limits::infinity(); + } + + //! Computes reciprocal of direction vector. + static constexpr BVH_VecNt computeInvDirect(const BVH_VecNt& theDirect) noexcept + { + static_assert(N >= 1 && N <= 4, "BVH_Ray only supports dimensions 1 to 4"); + if constexpr (N == 1) + { + return BVH_VecNt(invComponent(theDirect[0])); + } + else if constexpr (N == 2) + { + return BVH_VecNt(invComponent(theDirect.x()), invComponent(theDirect.y())); + } + else if constexpr (N == 3) + { + return BVH_VecNt(invComponent(theDirect.x()), + invComponent(theDirect.y()), + invComponent(theDirect.z())); + } + else // N == 4 + { + return BVH_VecNt(invComponent(theDirect.x()), + invComponent(theDirect.y()), + invComponent(theDirect.z()), + invComponent(theDirect.w())); + } + } }; #endif // _BVH_Ray_Header diff --git a/src/FoundationClasses/TKMath/BVH/BVH_Tools.hxx b/src/FoundationClasses/TKMath/BVH/BVH_Tools.hxx index eef6e2f4d3..4eb5592299 100644 --- a/src/FoundationClasses/TKMath/BVH/BVH_Tools.hxx +++ b/src/FoundationClasses/TKMath/BVH/BVH_Tools.hxx @@ -279,7 +279,8 @@ public: //! @name Point-Triangle Square distance } public: //! @name Ray-Box Intersection - //! Computes hit time of ray-box intersection + //! Computes hit time of ray-box intersection. + //! Uses precomputed reciprocal direction from BVH_Ray for optimal performance. static Standard_Boolean RayBoxIntersection(const BVH_Ray& theRay, const BVH_Box& theBox, T& theTimeEnter, @@ -297,18 +298,48 @@ public: //! @name Ray-Box Intersection } //! Computes hit time of ray-box intersection. + //! Uses precomputed reciprocal direction from BVH_Ray for optimal performance. + //! Handles parallel rays (infinite inverse direction) explicitly to avoid NaN from 0*inf + //! when ray origin is exactly on a slab boundary. static Standard_Boolean RayBoxIntersection(const BVH_Ray& theRay, const BVH_VecNt& theBoxCMin, const BVH_VecNt& theBoxCMax, T& theTimeEnter, T& theTimeLeave) { - return RayBoxIntersection(theRay.Origin, - theRay.Direct, - theBoxCMin, - theBoxCMax, - theTimeEnter, - theTimeLeave); + T aTimeEnter = (std::numeric_limits::lowest)(); + T aTimeLeave = (std::numeric_limits::max)(); + + for (int i = 0; i < N; ++i) + { + // Handle parallel rays (infinite inverse direction) to avoid NaN from 0*inf + if (std::isinf(theRay.InvDirect[i])) + { + if (theRay.Origin[i] < theBoxCMin[i] || theRay.Origin[i] > theBoxCMax[i]) + { + return Standard_False; + } + continue; + } + T aT1 = (theBoxCMin[i] - theRay.Origin[i]) * theRay.InvDirect[i]; + T aT2 = (theBoxCMax[i] - theRay.Origin[i]) * theRay.InvDirect[i]; + aTimeEnter = (std::max)(aTimeEnter, (std::min)(aT1, aT2)); + aTimeLeave = (std::min)(aTimeLeave, (std::max)(aT1, aT2)); + if (aTimeEnter > aTimeLeave) + { + return Standard_False; + } + } + + // Check if intersection is behind the ray origin + if (aTimeLeave < T(0)) + { + return Standard_False; + } + + theTimeEnter = aTimeEnter; + theTimeLeave = aTimeLeave; + return Standard_True; } //! Computes hit time of ray-box intersection @@ -332,6 +363,13 @@ public: //! @name Ray-Box Intersection //! Computes hit time of ray-box intersection. //! Uses optimized single-pass algorithm with early exit. + //! @param theRayOrigin ray origin point + //! @param theRayDirection ray direction vector + //! @param theBoxCMin minimum corner of the box + //! @param theBoxCMax maximum corner of the box + //! @param theTimeEnter time of ray entering the box + //! @param theTimeLeave time of ray leaving the box + //! @return true if ray intersects the box static Standard_Boolean RayBoxIntersection(const BVH_VecNt& theRayOrigin, const BVH_VecNt& theRayDirection, const BVH_VecNt& theBoxCMin, @@ -386,4 +424,4 @@ public: //! @name Ray-Box Intersection } }; -#endif \ No newline at end of file +#endif // _BVH_Tools_Header