Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-09-25 09:19:25

0001 // Created by: Eugeny MALTCHIKOV
0002 // Created on: 2019-04-17
0003 // Copyright (c) 2019 OPEN CASCADE SAS
0004 //
0005 // This file is part of Open CASCADE Technology software library.
0006 //
0007 // This library is free software; you can redistribute it and/or modify it under
0008 // the terms of the GNU Lesser General Public License version 2.1 as published
0009 // by the Free Software Foundation, with special exception defined in the file
0010 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0011 // distribution for complete text of the license and disclaimer of any warranty.
0012 //
0013 // Alternatively, this file may be used under the terms of Open CASCADE
0014 // commercial license or contractual agreement.
0015 
0016 #ifndef _BVH_Tools_Header
0017 #define _BVH_Tools_Header
0018 
0019 #include <BVH_Box.hxx>
0020 #include <BVH_Ray.hxx>
0021 #include <BVH_Types.hxx>
0022 
0023 //! Defines a set of static methods operating with points and bounding boxes.
0024 //! \tparam T Numeric data type
0025 //! \tparam N Vector dimension
0026 template <class T, int N>
0027 class BVH_Tools
0028 {
0029 public: //! @name public types
0030   typedef typename BVH::VectorType<T, N>::Type BVH_VecNt;
0031 
0032 public:
0033   enum BVH_PrjStateInTriangle
0034   {
0035     BVH_PrjStateInTriangle_VERTEX,
0036     BVH_PrjStateInTriangle_EDGE,
0037     BVH_PrjStateInTriangle_INNER
0038   };
0039 
0040 public: //! @name Box-Box Square distance
0041   //! Computes Square distance between Axis aligned bounding boxes
0042   static T BoxBoxSquareDistance(const BVH_Box<T, N>& theBox1, const BVH_Box<T, N>& theBox2)
0043   {
0044     if (!theBox1.IsValid() || !theBox2.IsValid())
0045     {
0046       return static_cast<T>(0);
0047     }
0048     return BoxBoxSquareDistance(theBox1.CornerMin(),
0049                                 theBox1.CornerMax(),
0050                                 theBox2.CornerMin(),
0051                                 theBox2.CornerMax());
0052   }
0053 
0054   //! Computes Square distance between Axis aligned bounding boxes
0055   static T BoxBoxSquareDistance(const BVH_VecNt& theCMin1,
0056                                 const BVH_VecNt& theCMax1,
0057                                 const BVH_VecNt& theCMin2,
0058                                 const BVH_VecNt& theCMax2)
0059   {
0060     T aDist = T(0);
0061     for (int i = 0; i < N; ++i)
0062     {
0063       if (theCMin1[i] > theCMax2[i])
0064       {
0065         T d = theCMin1[i] - theCMax2[i];
0066         aDist += d * d;
0067       }
0068       else if (theCMax1[i] < theCMin2[i])
0069       {
0070         T d = theCMin2[i] - theCMax1[i];
0071         aDist += d * d;
0072       }
0073     }
0074     return aDist;
0075   }
0076 
0077 public: //! @name Point-Box Square distance
0078   //! Computes square distance between point and bounding box
0079   static T PointBoxSquareDistance(const BVH_VecNt& thePoint, const BVH_Box<T, N>& theBox)
0080   {
0081     if (!theBox.IsValid())
0082     {
0083       return static_cast<T>(0);
0084     }
0085     return PointBoxSquareDistance(thePoint, theBox.CornerMin(), theBox.CornerMax());
0086   }
0087 
0088   //! Computes square distance between point and bounding box
0089   static T PointBoxSquareDistance(const BVH_VecNt& thePoint,
0090                                   const BVH_VecNt& theCMin,
0091                                   const BVH_VecNt& theCMax)
0092   {
0093     T aDist = T(0);
0094     for (int i = 0; i < N; ++i)
0095     {
0096       if (thePoint[i] < theCMin[i])
0097       {
0098         T d = theCMin[i] - thePoint[i];
0099         aDist += d * d;
0100       }
0101       else if (thePoint[i] > theCMax[i])
0102       {
0103         T d = thePoint[i] - theCMax[i];
0104         aDist += d * d;
0105       }
0106     }
0107     return aDist;
0108   }
0109 
0110 public: //! @name Point-Box projection
0111   //! Computes projection of point on bounding box
0112   static BVH_VecNt PointBoxProjection(const BVH_VecNt& thePoint, const BVH_Box<T, N>& theBox)
0113   {
0114     if (!theBox.IsValid())
0115     {
0116       return thePoint;
0117     }
0118     return PointBoxProjection(thePoint, theBox.CornerMin(), theBox.CornerMax());
0119   }
0120 
0121   //! Computes projection of point on bounding box
0122   static BVH_VecNt PointBoxProjection(const BVH_VecNt& thePoint,
0123                                       const BVH_VecNt& theCMin,
0124                                       const BVH_VecNt& theCMax)
0125   {
0126     return thePoint.cwiseMax(theCMin).cwiseMin(theCMax);
0127   }
0128 
0129 private: //! @name Internal helpers for point-triangle projection
0130   //! Helper to set projection state for vertex
0131   static void SetVertexState(BVH_PrjStateInTriangle* thePrjState,
0132                              int*                    theFirstNode,
0133                              int*                    theLastNode,
0134                              int                     theVertexIndex)
0135   {
0136     if (thePrjState != nullptr)
0137     {
0138       *thePrjState  = BVH_PrjStateInTriangle_VERTEX;
0139       *theFirstNode = theVertexIndex;
0140       *theLastNode  = theVertexIndex;
0141     }
0142   }
0143 
0144   //! Helper to set projection state for edge
0145   static void SetEdgeState(BVH_PrjStateInTriangle* thePrjState,
0146                            int*                    theFirstNode,
0147                            int*                    theLastNode,
0148                            int                     theStartVertex,
0149                            int                     theEndVertex)
0150   {
0151     if (thePrjState != nullptr)
0152     {
0153       *thePrjState  = BVH_PrjStateInTriangle_EDGE;
0154       *theFirstNode = theStartVertex;
0155       *theLastNode  = theEndVertex;
0156     }
0157   }
0158 
0159   //! Helper to compute projection onto edge
0160   static BVH_VecNt ProjectToEdge(const BVH_VecNt& theEdgeStart,
0161                                  const BVH_VecNt& theEdge,
0162                                  T                theDot1,
0163                                  T                theDot2)
0164   {
0165     T aT = theDot1 / (theDot1 + theDot2);
0166     return theEdgeStart + theEdge * aT;
0167   }
0168 
0169 public: //! @name Point-Triangle Square distance
0170   //! Find nearest point on a triangle for the given point.
0171   //! Uses Voronoi region testing to determine closest feature (vertex, edge, or interior).
0172   static BVH_VecNt PointTriangleProjection(const BVH_VecNt&        thePoint,
0173                                            const BVH_VecNt&        theNode0,
0174                                            const BVH_VecNt&        theNode1,
0175                                            const BVH_VecNt&        theNode2,
0176                                            BVH_PrjStateInTriangle* thePrjState          = nullptr,
0177                                            int*                    theNumberOfFirstNode = nullptr,
0178                                            int*                    theNumberOfLastNode  = nullptr)
0179   {
0180     // Compute edge vectors
0181     const BVH_VecNt aAB = theNode1 - theNode0;
0182     const BVH_VecNt aAC = theNode2 - theNode0;
0183     const BVH_VecNt aBC = theNode2 - theNode1;
0184 
0185     // Compute point-to-vertex vectors
0186     const BVH_VecNt aAP = thePoint - theNode0;
0187     const BVH_VecNt aBP = thePoint - theNode1;
0188     const BVH_VecNt aCP = thePoint - theNode2;
0189 
0190     // Compute dot products for Voronoi region tests
0191     const T aABdotAP = aAB.Dot(aAP);
0192     const T aACdotAP = aAC.Dot(aAP);
0193 
0194     // Check if P is in vertex region outside A
0195     if (aABdotAP <= T(0) && aACdotAP <= T(0))
0196     {
0197       SetVertexState(thePrjState, theNumberOfFirstNode, theNumberOfLastNode, 0);
0198       return theNode0;
0199     }
0200 
0201     const T aBAdotBP = -aAB.Dot(aBP);
0202     const T aBCdotBP = aBC.Dot(aBP);
0203 
0204     // Check if P is in vertex region outside B
0205     if (aBAdotBP <= T(0) && aBCdotBP <= T(0))
0206     {
0207       SetVertexState(thePrjState, theNumberOfFirstNode, theNumberOfLastNode, 1);
0208       return theNode1;
0209     }
0210 
0211     const T aCBdotCP = -aBC.Dot(aCP);
0212     const T aCAdotCP = -aAC.Dot(aCP);
0213 
0214     // Check if P is in vertex region outside C
0215     if (aCAdotCP <= T(0) && aCBdotCP <= T(0))
0216     {
0217       SetVertexState(thePrjState, theNumberOfFirstNode, theNumberOfLastNode, 2);
0218       return theNode2;
0219     }
0220 
0221     // Compute barycentric coordinates for edge/interior tests
0222     const T aACdotBP = aAC.Dot(aBP);
0223     const T aVC      = aABdotAP * aACdotBP + aBAdotBP * aACdotAP;
0224 
0225     // Check if P is in edge region of AB
0226     if (aVC <= T(0) && aABdotAP > T(0) && aBAdotBP > T(0))
0227     {
0228       SetEdgeState(thePrjState, theNumberOfFirstNode, theNumberOfLastNode, 0, 1);
0229       return ProjectToEdge(theNode0, aAB, aABdotAP, aBAdotBP);
0230     }
0231 
0232     const T aABdotCP = aAB.Dot(aCP);
0233     const T aVA      = aBAdotBP * aCAdotCP - aABdotCP * aACdotBP;
0234 
0235     // Check if P is in edge region of BC
0236     if (aVA <= T(0) && aBCdotBP > T(0) && aCBdotCP > T(0))
0237     {
0238       SetEdgeState(thePrjState, theNumberOfFirstNode, theNumberOfLastNode, 1, 2);
0239       return ProjectToEdge(theNode1, aBC, aBCdotBP, aCBdotCP);
0240     }
0241 
0242     const T aVB = aABdotCP * aACdotAP + aABdotAP * aCAdotCP;
0243 
0244     // Check if P is in edge region of CA
0245     if (aVB <= T(0) && aACdotAP > T(0) && aCAdotCP > T(0))
0246     {
0247       SetEdgeState(thePrjState, theNumberOfFirstNode, theNumberOfLastNode, 2, 0);
0248       return ProjectToEdge(theNode0, aAC, aACdotAP, aCAdotCP);
0249     }
0250 
0251     // P is inside triangle - compute barycentric coordinates
0252     const T aNorm = aVA + aVB + aVC;
0253 
0254     // Handle degenerate triangle (zero or near-zero area)
0255     if (aNorm
0256         <= std::numeric_limits<T>::epsilon() * (std::abs(aVA) + std::abs(aVB) + std::abs(aVC)))
0257     {
0258       SetVertexState(thePrjState, theNumberOfFirstNode, theNumberOfLastNode, 0);
0259       return (theNode0 + theNode1 + theNode2) / T(3);
0260     }
0261 
0262     if (thePrjState != nullptr)
0263     {
0264       *thePrjState = BVH_PrjStateInTriangle_INNER;
0265     }
0266 
0267     return (theNode0 * aVA + theNode1 * aVB + theNode2 * aVC) / aNorm;
0268   }
0269 
0270   //! Computes square distance between point and triangle
0271   static T PointTriangleSquareDistance(const BVH_VecNt& thePoint,
0272                                        const BVH_VecNt& theNode0,
0273                                        const BVH_VecNt& theNode1,
0274                                        const BVH_VecNt& theNode2)
0275   {
0276     const BVH_VecNt aProj = PointTriangleProjection(thePoint, theNode0, theNode1, theNode2);
0277     const BVH_VecNt aPP   = aProj - thePoint;
0278     return aPP.Dot(aPP);
0279   }
0280 
0281 public: //! @name Ray-Box Intersection
0282   //! Computes hit time of ray-box intersection.
0283   //! Uses precomputed reciprocal direction from BVH_Ray for optimal performance.
0284   static bool RayBoxIntersection(const BVH_Ray<T, N>& theRay,
0285                                  const BVH_Box<T, N>& theBox,
0286                                  T&                   theTimeEnter,
0287                                  T&                   theTimeLeave)
0288   {
0289     if (!theBox.IsValid())
0290     {
0291       return false;
0292     }
0293     return RayBoxIntersection(theRay,
0294                               theBox.CornerMin(),
0295                               theBox.CornerMax(),
0296                               theTimeEnter,
0297                               theTimeLeave);
0298   }
0299 
0300   //! Computes hit time of ray-box intersection.
0301   //! Uses precomputed reciprocal direction from BVH_Ray for optimal performance.
0302   //! Handles parallel rays (infinite inverse direction) explicitly to avoid NaN from 0*inf
0303   //! when ray origin is exactly on a slab boundary.
0304   static bool RayBoxIntersection(const BVH_Ray<T, N>& theRay,
0305                                  const BVH_VecNt&     theBoxCMin,
0306                                  const BVH_VecNt&     theBoxCMax,
0307                                  T&                   theTimeEnter,
0308                                  T&                   theTimeLeave)
0309   {
0310     T aTimeEnter = (std::numeric_limits<T>::lowest)();
0311     T aTimeLeave = (std::numeric_limits<T>::max)();
0312 
0313     for (int i = 0; i < N; ++i)
0314     {
0315       // Handle parallel rays (infinite inverse direction) to avoid NaN from 0*inf
0316       if (std::isinf(theRay.InvDirect[i]))
0317       {
0318         if (theRay.Origin[i] < theBoxCMin[i] || theRay.Origin[i] > theBoxCMax[i])
0319         {
0320           return false;
0321         }
0322         continue;
0323       }
0324       T aT1      = (theBoxCMin[i] - theRay.Origin[i]) * theRay.InvDirect[i];
0325       T aT2      = (theBoxCMax[i] - theRay.Origin[i]) * theRay.InvDirect[i];
0326       aTimeEnter = (std::max)(aTimeEnter, (std::min)(aT1, aT2));
0327       aTimeLeave = (std::min)(aTimeLeave, (std::max)(aT1, aT2));
0328       if (aTimeEnter > aTimeLeave)
0329       {
0330         return false;
0331       }
0332     }
0333 
0334     // Check if intersection is behind the ray origin
0335     if (aTimeLeave < T(0))
0336     {
0337       return false;
0338     }
0339 
0340     theTimeEnter = aTimeEnter;
0341     theTimeLeave = aTimeLeave;
0342     return true;
0343   }
0344 
0345   //! Computes hit time of ray-box intersection
0346   static bool RayBoxIntersection(const BVH_VecNt&     theRayOrigin,
0347                                  const BVH_VecNt&     theRayDirection,
0348                                  const BVH_Box<T, N>& theBox,
0349                                  T&                   theTimeEnter,
0350                                  T&                   theTimeLeave)
0351   {
0352     if (!theBox.IsValid())
0353     {
0354       return false;
0355     }
0356     return RayBoxIntersection(theRayOrigin,
0357                               theRayDirection,
0358                               theBox.CornerMin(),
0359                               theBox.CornerMax(),
0360                               theTimeEnter,
0361                               theTimeLeave);
0362   }
0363 
0364   //! Computes hit time of ray-box intersection.
0365   //! Uses optimized single-pass algorithm with early exit.
0366   //! @param theRayOrigin ray origin point
0367   //! @param theRayDirection ray direction vector
0368   //! @param theBoxCMin minimum corner of the box
0369   //! @param theBoxCMax maximum corner of the box
0370   //! @param theTimeEnter time of ray entering the box
0371   //! @param theTimeLeave time of ray leaving the box
0372   //! @return true if ray intersects the box
0373   static bool RayBoxIntersection(const BVH_VecNt& theRayOrigin,
0374                                  const BVH_VecNt& theRayDirection,
0375                                  const BVH_VecNt& theBoxCMin,
0376                                  const BVH_VecNt& theBoxCMax,
0377                                  T&               theTimeEnter,
0378                                  T&               theTimeLeave)
0379   {
0380     T aTimeEnter = (std::numeric_limits<T>::lowest)();
0381     T aTimeLeave = (std::numeric_limits<T>::max)();
0382 
0383     for (int i = 0; i < N; ++i)
0384     {
0385       if (theRayDirection[i] == T(0))
0386       {
0387         // Ray is parallel to this axis slab - check if origin is within bounds
0388         if (theRayOrigin[i] < theBoxCMin[i] || theRayOrigin[i] > theBoxCMax[i])
0389         {
0390           return false; // Ray misses the slab entirely
0391         }
0392         // Ray is within the slab, doesn't constrain the intersection interval
0393         continue;
0394       }
0395 
0396       // Compute intersection distances for this axis
0397       T aT1 = (theBoxCMin[i] - theRayOrigin[i]) / theRayDirection[i];
0398       T aT2 = (theBoxCMax[i] - theRayOrigin[i]) / theRayDirection[i];
0399 
0400       // Ensure aT1 <= aT2 (handle negative direction)
0401       T aTMin = (std::min)(aT1, aT2);
0402       T aTMax = (std::max)(aT1, aT2);
0403 
0404       // Update intersection interval
0405       aTimeEnter = (std::max)(aTimeEnter, aTMin);
0406       aTimeLeave = (std::min)(aTimeLeave, aTMax);
0407 
0408       // Early exit if no intersection
0409       if (aTimeEnter > aTimeLeave)
0410       {
0411         return false;
0412       }
0413     }
0414 
0415     // Check if intersection is behind the ray origin
0416     if (aTimeLeave < T(0))
0417     {
0418       return false;
0419     }
0420 
0421     theTimeEnter = aTimeEnter;
0422     theTimeLeave = aTimeLeave;
0423     return true;
0424   }
0425 };
0426 
0427 #endif // _BVH_Tools_Header