Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-19 09:40:14 

0001 /// \file BVH.h
0002 /// \author Guilherme Amadio
0003 
0004 #ifndef VECGEOM_BASE_BVH_H_
0005 #define VECGEOM_BASE_BVH_H_
0006 
0007 #include "VecGeom/base/AABB.h"
0008 #include "VecGeom/base/Config.h"
0009 #include "VecGeom/base/Cuda.h"
0010 #include "VecGeom/navigation/NavStateIndex.h"
0011 #include "VecGeom/navigation/NavigationState.h"
0012 #include "VecGeom/volumes/LogicalVolume.h"
0013 #include "VecGeom/volumes/PlacedVolume.h"
0014 // #include "VecGeom/surfaces/Model.h"
0015 
0016 #include <vector>
0017 
0018 // Forward-declare CPUsurfData
0019 namespace vgbrep {
0020 template <typename Real_t>
0021 struct CPUsurfData;
0022 }
0023 
0024 namespace vecgeom {
0025 namespace cuda {
0026 template <typename Real_t>
0027 class BVH;
0028 }
0029 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE(class, BVH, typename);
0030 inline namespace VECGEOM_IMPL_NAMESPACE {
0031 
0032 class LogicalVolume;
0033 class VPlacedVolume;
0034 
0035 /**
0036  * @brief Bounding Volume Hierarchy class to represent an axis-aligned bounding volume hierarchy.
0037  * @details BVH instances can be associated with logical volumes to accelerate queries to their child volumes.
0038  */
0039 template <typename Real_t>
0040 class BVH {
0041 private:
0042   uint fRootId    = 0;           ///< Id of the root element this BVH was constructed for
0043   int fRootNChild = 0;           ///< Number of children of the root element
0044   int fDepth      = 0;           ///< Depth of the BVH
0045   int *fPrimId{nullptr};         ///< Child volume ids for each BVH node
0046   int *fOffset{nullptr};         ///< Offset in @c fPrimId for first child of each BVH node
0047   int *fNChild{nullptr};         ///< Number of children for each BVH node
0048   AABB<Real_t> *fNodes{nullptr}; ///< AABBs of BVH nodes
0049   AABB<Real_t> *fAABBs{nullptr}; ///< AABBs of children of the BVH root element
0050 
0051 public:
0052   // Default constructor
0053   BVH()
0054       : fRootId(0), fRootNChild(0), fDepth(0), fPrimId(nullptr), fOffset(nullptr), fNChild(nullptr), fNodes(nullptr),
0055         fAABBs(nullptr)
0056   {
0057   }
0058 
0059   uint GetRootId() const { return fRootId; }
0060   int GetRootNChild() const { return fRootNChild; };
0061   int GetDepth() const { return fDepth; };
0062   const int *GetPrimId() const { return fPrimId; };
0063   const int *GetOffset() const { return fOffset; };
0064   const int *GetNChild() const { return fNChild; };
0065   const AABB<Real_t> *GetAABBs() const { return fAABBs; };
0066   const AABB<Real_t> *GetNodes() const { return fNodes; };
0067 
0068   /** Maximum depth. */
0069   static constexpr int BVH_MAX_DEPTH = 32;
0070   /**
0071    * Constructor.
0072    * @param volume Pointer to logical volume for which the BVH will be created.
0073    * @param ptrAABB Container of AABBs for this volume
0074    * @param nChild Number of children of the volume
0075    * @param depth Depth of the BVH binary tree. Defaults to zero, in which case
0076    * the actual depth will be chosen dynamically based on the number of child volumes.
0077    * When a fixed depth is chosen, it cannot be larger than @p BVH_MAX_DEPTH.
0078    */
0079   BVH(LogicalVolume const &volume, Vector3D<Precision> *ptrAABB, int nChild, int depth = 0);
0080 
0081   /** Destructor. */
0082   ~BVH() { Clear(); }
0083   void Clear();
0084 
0085   /**
0086    * Initializer used by BVHcreator. Takes as input pre-constructed BVH buffers.
0087    * @param id  Id of the logical volume
0088    * @param nchild Number of children of the volume
0089    * @param depth Depth of the BVH binary tree stored in the device buffers.
0090    * @param dPrimId Device buffer with child volume ids
0091    * @param dAABBs  Device buffer with AABBs of child volumes
0092    * @param dOffset Device buffer with offsets in @c dPrimId for first child of each BVH node
0093    * @param dNChild Device buffer with number of children for each BVH node
0094    * @param dNodes AABBs of BVH nodes
0095    */
0096   VECCORE_ATT_HOST_DEVICE
0097   void Set(int id, int nchild, int depth, int *dPrimId, vecgeom::AABB<Real_t> *dAABBs, int *dOffset, int *dNChild,
0098            vecgeom::AABB<Real_t> *dNodes)
0099   {
0100     fRootId     = id;
0101     fRootNChild = nchild;
0102     fDepth      = depth;
0103     SetPointers(dPrimId, dOffset, dNChild, dAABBs, dNodes);
0104   }
0105 
0106   /**
0107    * Setter for all member arrays. Used to update the pointers after a copy from host to device
0108    * @param dPrimId Device buffer with child volume ids
0109    * @param dAABBs  Device buffer with AABBs of child volumes
0110    * @param dOffset Device buffer with offsets in @c dPrimId for first child of each BVH node
0111    * @param dNChild Device buffer with number of children for each BVH node
0112    * @param dNodes AABBs of BVH nodes
0113    */
0114   VECCORE_ATT_HOST_DEVICE
0115   void SetPointers(int *dPrimId, int *dOffset, int *dNChild, AABB<Real_t> *dAABBs, AABB<Real_t> *dNodes)
0116   {
0117     fPrimId = dPrimId;
0118     fOffset = dOffset;
0119     fNChild = dNChild;
0120     fAABBs  = dAABBs;
0121     fNodes  = dNodes;
0122   }
0123 
0124 #ifdef VECGEOM_ENABLE_CUDA
0125   /**
0126    * Constructor for GPU. Takes as input pre-constructed BVH buffers.
0127    * @param volume  Reference to logical volume on the device
0128    * @param depth Depth of the BVH binary tree stored in the device buffers.
0129    * @param dPrimId Device buffer with child volume ids
0130    * @param dAABBs  Device buffer with AABBs of child volumes
0131    * @param dOffset Device buffer with offsets in @c dPrimId for first child of each BVH node
0132    * @param dNChild Device buffer with number of children for each BVH node
0133    * @param dNodes AABBs of BVH nodes
0134    */
0135   VECCORE_ATT_DEVICE
0136   BVH(LogicalVolume const *volume, int depth, int *dPrimId, AABB<Real_t> *dAABBs, int *dOffset, int *NChild,
0137       AABB<Real_t> *dNodes);
0138 
0139   /**
0140    * Constructor for GPU. Takes as input pre-constructed BVH buffers.
0141    * @param id  Id of the logical volume
0142    * @param nchild Number of children of the volume
0143    * @param depth Depth of the BVH binary tree stored in the device buffers.
0144    * @param dPrimId Device buffer with child volume ids
0145    * @param dAABBs  Device buffer with AABBs of child volumes
0146    * @param dOffset Device buffer with offsets in @c dPrimId for first child of each BVH node
0147    * @param dNChild Device buffer with number of children for each BVH node
0148    * @param dNodes AABBs of BVH nodes
0149    */
0150   VECCORE_ATT_HOST_DEVICE
0151   BVH(int id, int nchild, int depth, int *dPrimId, vecgeom::AABB<Real_t> *dAABBs, int *dOffset, int *dNChild,
0152       vecgeom::AABB<Real_t> *dNodes)
0153       : fRootId(id), fRootNChild(nchild), fDepth(depth), fPrimId(dPrimId), fOffset(dOffset), fNChild(dNChild),
0154         fNodes(dNodes), fAABBs(dAABBs)
0155   {
0156   }
0157 
0158 #endif
0159 
0160 #ifdef VECGEOM_CUDA_INTERFACE
0161   /** Copy and construct an instance of this BVH on the device, at the device address @p addr. */
0162   DevicePtr<cuda::BVH<Real_t>> CopyToGpu(void *addr) const;
0163 #endif
0164 
0165   // void CopyToGpu(BVH *dBVH) const;
0166 
0167   /** Print a summary of BVH contents */
0168   VECCORE_ATT_HOST_DEVICE
0169   void Print(bool verbose = false) const;
0170 
0171   uint GetAllocatedSize() const
0172   {
0173     uint nodes = (2 << fDepth) - 1;
0174     uint size{0};
0175     // fAABBs
0176     size += fRootNChild * sizeof(vecgeom::AABB<Real_t>);
0177     // fNodes
0178     size += nodes * sizeof(vecgeom::AABB<Real_t>);
0179     // fNChild
0180     size += nodes * sizeof(int);
0181     // fOffset
0182     size += nodes * sizeof(int);
0183     // fPrimId
0184     size += fRootNChild * sizeof(int);
0185     return size;
0186   }
0187 
0188   /**
0189    * Check ray defined by <tt>localpoint + t * localdir</tt> for intersections with children
0190    * of the root element of the BVH, and within a maximum distance of @p step
0191    * along the ray, while ignoring the @p last_exited_id volume.
0192    * @param[in] localpoint Point in the local coordinates of the BVH root.
0193    * @param[in] localdir Direction in the local coordinates of the BVH root.
0194    * @param[in,out] step Maximum step distance for which intersections should be considered.
0195    * @param[in] last_exited_id Last exited element. This element is ignored when reporting intersections.
0196    * @param[out] hitcandidate_index Index of element for which closest intersection was found. -1 if no intersection
0197    * is found within the current step distance.
0198    */
0199   /*
0200    * BVH::ComputeDaughterIntersections() computes the intersection of a ray against all children of
0201    * the logical volume. A stack is kept of the node ids that need to be checked. It needs to be at
0202    * most as deep as the binary tree itself because we always first pop the current node, and then
0203    * add at most the two children. For example, for depth two, we pop the root node, then at most we
0204    * add both of its leaves onto the stack to be checked. We initialize ptr with &stack[1] such that
0205    * when we pop the first time as we enter the loop, the position we read from is the first position
0206    * of the stack, which contains the id 0 for the root node. When we pop the stack such that ptr
0207    * points before &stack[0], it means we've checked all we needed and the loop can be terminated.
0208    * In order to determine if a node of the tree is internal or not, we check if the node id of its
0209    * left child is past the end of the array (in which case we know we are at the maximum depth), or
0210    * if the sum of children in both leaves is the same as in the current node, as for leaf nodes, the
0211    * sum of children in the left+right child nodes will be less than for the current node.
0212    */
0213   template <typename Navigator, typename Real_i>
0214   VECCORE_ATT_HOST_DEVICE void CheckDaughterIntersections(const Vector3D<Real_i> &localpoint,
0215                                                           const Vector3D<Real_i> &localdir, Real_i &step,
0216                                                           long const last_exited_id, long &hitcandidate_index) const
0217   {
0218     unsigned int stack[BVH_MAX_DEPTH], *ptr = &stack[1];
0219     stack[0] = 0;
0220 
0221     /* Calculate and reuse inverse direction to save on divisions */
0222     Vector3D<Real_t> binvdir(static_cast<Real_t>(1.0) / vecgeom::NonZero(localdir[0]),
0223                              static_cast<Real_t>(1.0) / vecgeom::NonZero(localdir[1]),
0224                              static_cast<Real_t>(1.0) / vecgeom::NonZero(localdir[2]));
0225     Vector3D<Real_t> blocalpoint(static_cast<Real_t>(localpoint[0]), static_cast<Real_t>(localpoint[1]),
0226                                  static_cast<Real_t>(localpoint[2]));
0227     Vector3D<Real_t> blocaldir(static_cast<Real_t>(localdir[0]), static_cast<Real_t>(localdir[1]),
0228                                static_cast<Real_t>(localdir[2]));
0229     Real_t bstep = static_cast<Real_t>(step);
0230 
0231     do {
0232       const unsigned int id = *--ptr; /* pop next node id to be checked from the stack */
0233 
0234       // If the current distance is shorter than the distance to the node we can safely ignore it
0235       Real_t min{vecgeom::InfinityLength<Real_t>()}, max{-vecgeom::InfinityLength<Real_t>()};
0236       fNodes[id].ComputeIntersectionInvDir(blocalpoint, binvdir, min, max);
0237       if (min > max || max < Real_t{0.} || min >= step) {
0238         continue;
0239       }
0240 
0241       if (fNChild[id] >= 0) {
0242 
0243         /* For leaf nodes, loop over children */
0244         for (int i = 0; i < fNChild[id]; ++i) {
0245           const int prim = fPrimId[fOffset[id] + i];
0246           if (last_exited_id >= 0 && Navigator::SkipItem(fRootId, prim, last_exited_id)) continue;
0247           /* Check AABB first, then the element itself if needed */
0248           Real_t approach;
0249           if (fAABBs[prim].IntersectInvDirApproach(blocalpoint, binvdir, bstep, approach)) {
0250             auto dist = Navigator::CandidateDistanceToIn(
0251                 fRootId, prim, localpoint + static_cast<Real_i>(approach) * localdir, localdir, step);
0252             // Only compensate with the approach distance if the distance is positive (i.e. not a wrong-side error)
0253             dist += (dist > 0.) * static_cast<Real_i>(approach);
0254             /* If distance to current child is smaller than current step, update step and hitcandidate */
0255             if (dist < step) {
0256               dist               = vecCore::Max(dist, 0.);
0257               step               = dist;
0258               bstep              = static_cast<Real_t>(dist);
0259               hitcandidate_index = prim;
0260             }
0261           }
0262         }
0263       } else {
0264         const unsigned int childL = 2 * id + 1;
0265         const unsigned int childR = 2 * id + 2;
0266 
0267         /* For internal nodes, check AABBs to know if we need to traverse left and right children */
0268         Real_t tminL = vecgeom::InfinityLength<Real_t>(), tmaxL = -vecgeom::InfinityLength<Real_t>(),
0269                tminR = vecgeom::InfinityLength<Real_t>(), tmaxR = -vecgeom::InfinityLength<Real_t>();
0270 
0271         fNodes[childL].ComputeIntersectionInvDir(blocalpoint, binvdir, tminL, tmaxL);
0272         fNodes[childR].ComputeIntersectionInvDir(blocalpoint, binvdir, tminR, tmaxR);
0273 
0274         const bool traverseL = tminL <= tmaxL && tmaxL >= static_cast<Real_t>(0.0) && tminL < bstep;
0275         const bool traverseR = tminR <= tmaxR && tmaxR >= static_cast<Real_t>(0.0) && tminR < bstep;
0276 
0277         /*
0278          * If both left and right nodes need to be checked, check closest one first.
0279          * This ensures step gets short as fast as possible so we can skip more nodes without checking.
0280          */
0281         if (tminR < tminL) {
0282           if (traverseL) *ptr++ = childL;
0283           if (traverseR) *ptr++ = childR;
0284         } else {
0285           if (traverseR) *ptr++ = childR;
0286           if (traverseL) *ptr++ = childL;
0287         }
0288       }
0289     } while (ptr > stack);
0290   }
0291 
0292   /**
0293    * Compute safety against children of the root element associated with the BVH.
0294    * @param[in] localpoint Point in the local coordinates of the root element.
0295    * @param[in] safety Maximum safety. Elements further than this are not checked.
0296    * @param[in] limit Do not call the primitive safety if farther than this value
0297    * @returns Minimum between safety to the closest child of root element and input @p safety.
0298    */
0299   /*
0300    * BVH::ComputeSafety is very similar to CheckDaughterIntersections regarding traversal of the tree, but
0301    * it computes only the safety instead of the intersection using a ray, so the logic is a bit simpler.
0302    */
0303   template <typename Navigator>
0304   VECCORE_ATT_HOST_DEVICE Precision ComputeSafety(Vector3D<Precision> localpoint, Precision safety,
0305                                                   Precision limit = InfinityLength<Precision>()) const
0306   {
0307     unsigned int stack[BVH_MAX_DEPTH], *ptr = &stack[1];
0308     stack[0] = 0;
0309 
0310     do {
0311       const unsigned int id = *--ptr;
0312 
0313       // We can safely ignore nodes that are farther than the current safety
0314       if (fNodes[id].Safety(localpoint) > safety) continue;
0315 
0316       if (fNChild[id] >= 0) {
0317         for (int i = 0; i < fNChild[id]; ++i) {
0318           const int prim = fPrimId[fOffset[id] + i];
0319           if (fAABBs[prim].Safety(localpoint) < safety) {
0320             auto safety_leaf = fAABBs[prim].Safety(localpoint);
0321             // If the distance to the current node is larger than the safety we can ignore it
0322             if (safety_leaf >= safety) continue;
0323             // Don't check daughters if the safety is larger than the accuracy limit
0324             if (safety_leaf > limit) {
0325               safety = safety_leaf;
0326               continue;
0327             }
0328             const Precision dist = Navigator::CandidateSafetyToIn(fRootId, prim, localpoint);
0329             // FIXME: A check for negative distances is needed for using the BVH with surfaces,
0330             // however, it causes unexpected navigation issues with solids
0331             // if (dist > -vecgeom::kToleranceDist<Precision>) {
0332             if (dist < safety) safety = dist;
0333             // }
0334           }
0335         }
0336       } else {
0337         const unsigned int childL = 2 * id + 1;
0338         const unsigned int childR = 2 * id + 2;
0339 
0340         const Real_t safetyL = fNodes[childL].Safety(localpoint);
0341         const Real_t safetyR = fNodes[childR].Safety(localpoint);
0342 
0343         const bool traverseL = safetyL < safety;
0344         const bool traverseR = safetyR < safety;
0345 
0346         if (safetyR < safetyL) {
0347           if (traverseR) *ptr++ = childR;
0348           if (traverseL) *ptr++ = childL;
0349         } else {
0350           if (traverseL) *ptr++ = childL;
0351           if (traverseR) *ptr++ = childR;
0352         }
0353       }
0354     } while (ptr > stack);
0355 
0356     return safety;
0357   }
0358 
0359   /**
0360    * Find child element inside which the given point @p localpoint is located.
0361    * @param[in] exclude_item_id Element that should be ignored.
0362    * @param[in] localpoint Point in the local coordinates of the BVH root element.
0363    * @param[out] container_id Id of the element in which @p localpoint is contained
0364    * @param[out] path Navigation state of the container element
0365    * @returns Whether @p localpoint falls within a child element of this BVH.
0366    */
0367   template <typename Navigator>
0368   VECCORE_ATT_HOST_DEVICE bool LevelLocate(int const exclude_item_id, Vector3D<Real_t> const &localpoint,
0369                                            int &container_id, vecgeom::NavigationState &path) const
0370   {
0371     unsigned int stack[BVH_MAX_DEPTH], *ptr = &stack[1];
0372     stack[0] = 0;
0373 
0374     do {
0375       const unsigned int id = *--ptr;
0376 
0377       if (fNChild[id] >= 0) {
0378         for (int i = 0; i < fNChild[id]; ++i) {
0379           const int prim = fPrimId[fOffset[id] + i];
0380           if (fAABBs[prim].Contains(localpoint)) {
0381             if (!Navigator::SkipItem(fRootId, prim, exclude_item_id) &&
0382                 Navigator::CandidateContains(fRootId, prim, localpoint, path)) {
0383               container_id = Navigator::ItemId(fRootId, prim);
0384               return true;
0385             }
0386           }
0387         }
0388       } else {
0389         const unsigned int childL = 2 * id + 1;
0390         if (fNodes[childL].Contains(localpoint)) *ptr++ = childL;
0391 
0392         const unsigned int childR = 2 * id + 2;
0393         if (fNodes[childR].Contains(localpoint)) *ptr++ = childR;
0394       }
0395     } while (ptr > stack);
0396 
0397     return false;
0398   }
0399 
0400   /**
0401    * Find child element inside which the given point @p localpoint is located.
0402    * @param[in] exclude_item_id Element that should be ignored.
0403    * @param[in] localpoint Point in the local coordinates of the BVH root element.
0404    * @param[out] container_id Id of the element in which @p localpoint is contained
0405    * @param[out] daughterlocalpoint Point in the local coordinates of the container element
0406    * @returns Whether @p localpoint falls within a child element of this BVH.
0407    */
0408   template <typename Navigator>
0409   VECCORE_ATT_HOST_DEVICE vecgeom::Inside_t LevelInside(long const exclude_item_id,
0410                                                         Vector3D<Precision> const &localpoint, long &container_id,
0411                                                         Vector3D<Precision> &daughterlocalpoint) const
0412   {
0413     unsigned int stack[BVH_MAX_DEPTH], *ptr = &stack[1];
0414     stack[0] = 0;
0415 
0416     do {
0417       const unsigned int id = *--ptr;
0418 
0419       if (fNChild[id] >= 0) {
0420         for (int i = 0; i < fNChild[id]; ++i) {
0421           const int prim = fPrimId[fOffset[id] + i];
0422           if (fAABBs[prim].Contains(localpoint)) {
0423             if (Navigator::SkipItem(fRootId, prim, exclude_item_id)) continue;
0424             auto inside = Navigator::CandidateInside(fRootId, prim, localpoint, daughterlocalpoint);
0425             if (inside != kOutside) {
0426               container_id = Navigator::ItemId(fRootId, prim);
0427               return inside;
0428             }
0429           }
0430         }
0431       } else {
0432         const unsigned int childL = 2 * id + 1;
0433         if (fNodes[childL].Contains(localpoint)) *ptr++ = childL;
0434 
0435         const unsigned int childR = 2 * id + 2;
0436         if (fNodes[childR].Contains(localpoint)) *ptr++ = childR;
0437       }
0438     } while (ptr > stack);
0439 
0440     return kOutside;
0441   }
0442 
0443   /**
0444    * Check ray defined by <tt>localpoint + t * localdir</tt> for intersections with bounding
0445    * boxes of children of the root element of the BVH, and within a maximum
0446    * distance of @p step along the ray. Returns the distance to the first crossed box.
0447    * @param[in] localpoint Point in the local coordinates of the root element.
0448    * @param[in] localdir Direction in the local coordinates of the root element.
0449    * @param[in,out] step Maximum step distance for which intersections should be considered.
0450    * @param[in] last_exited_id Last exited element. This element is ignored when reporting intersections.
0451    */
0452   /*
0453    * BVH::ApproachNextDaughter is very similar to CheckDaughterIntersections but computes the first
0454    * hit daughter bounding box instead of the next hit shape. This lighter computation is used to
0455    * first approach the next hit solid before computing the actual distance, in the attempt to
0456    * reduce the numerical rounding error due to propagation to boundary.
0457    */
0458   template <typename Navigator>
0459   VECCORE_ATT_HOST_DEVICE void ApproachNextDaughter(Vector3D<Precision> localpoint, Vector3D<Precision> localdir,
0460                                                     Precision &step, long const last_exited_id) const
0461   {
0462     unsigned int stack[BVH_MAX_DEPTH] = {0}, *ptr = &stack[1];
0463 
0464     /* Calculate and reuse inverse direction to save on divisions */
0465     Vector3D<Real_t> invlocaldir(static_cast<Real_t>(1.0 / NonZero(localdir[0])),
0466                                  static_cast<Real_t>(1.0 / NonZero(localdir[1])),
0467                                  static_cast<Real_t>(1.0 / NonZero(localdir[2])));
0468 
0469     do {
0470       unsigned int id = *--ptr; /* pop next node id to be checked from the stack */
0471 
0472       if (fNChild[id] >= 0) {
0473         /* For leaf nodes, loop over children */
0474         for (int i = 0; i < fNChild[id]; ++i) {
0475           int prim = fPrimId[fOffset[id] + i];
0476           /* Check vecgeom::AABB first, then the element itself if needed */
0477           if (fAABBs[prim].IntersectInvDir(localpoint, invlocaldir, step)) {
0478             const auto dist = Navigator::CandidateApproachSolid(fRootId, prim, localpoint, localdir);
0479             /* If distance to current child is smaller than current step, update step and hitcandidate */
0480             if (dist < step && !(dist <= 0.0 && Navigator::SkipItem(fRootId, prim, last_exited_id))) step = dist;
0481           }
0482         }
0483       } else {
0484         unsigned int childL = 2 * id + 1;
0485         unsigned int childR = 2 * id + 2;
0486 
0487         /* For internal nodes, check AABBs to know if we need to traverse left and right children */
0488         Real_t tminL = vecgeom::InfinityLength<Real_t>(), tmaxL = -vecgeom::InfinityLength<Real_t>(),
0489                tminR = vecgeom::InfinityLength<Real_t>(), tmaxR = -vecgeom::InfinityLength<Real_t>();
0490 
0491         fNodes[childL].ComputeIntersectionInvDir(localpoint, invlocaldir, tminL, tmaxL);
0492         fNodes[childR].ComputeIntersectionInvDir(localpoint, invlocaldir, tminR, tmaxR);
0493 
0494         bool traverseL = tminL <= tmaxL && tmaxL >= 0.0 && tminL < step;
0495         bool traverseR = tminR <= tmaxR && tmaxR >= 0.0 && tminR < step;
0496 
0497         /*
0498          * If both left and right nodes need to be checked, check closest one first.
0499          * This ensures step gets short as fast as possible so we can skip more nodes without checking.
0500          */
0501         if (tminR < tminL) {
0502           if (traverseR) *ptr++ = childR;
0503           if (traverseL) *ptr++ = childL;
0504         } else {
0505           if (traverseL) *ptr++ = childL;
0506           if (traverseR) *ptr++ = childR;
0507         }
0508       }
0509     } while (ptr > stack);
0510   }
0511 
0512   /**
0513    * Check ray defined by <tt>localpoint + t * localdir</tt> for intersections with children
0514    * of the root element of the BVH, and within a maximum distance of @p step
0515    * along the ray, while ignoring the @p last_exited_id volume.
0516    * @param[in] localpoint Point in the local coordinates of the BVH root.
0517    * @param[in] localdir Direction in the local coordinates of the BVH root.
0518    * @param[in,out] step Maximum step distance for which intersections should be considered.
0519    * @param[in] last_exited_id Last exited element. This element is ignored when reporting intersections.
0520    * @param[out] hitcandidate_index Index of element for which closest intersection was found. -1 if no intersection
0521    * is found within the current step distance.
0522    */
0523   /*
0524    * This function is meant to be used for benchmarking of the BVH, and not for actual navigation. It gathers stats
0525    * on the traversal of the BVH tree
0526    */
0527   template <typename Navigator, typename Real_i>
0528   VECCORE_ATT_HOST_DEVICE void CheckDaughterIntersectionsBenchmark(const Vector3D<Real_i> &localpoint,
0529                                                                    const Vector3D<Real_i> &localdir, Real_i &step,
0530                                                                    long const last_exited_id, long &hitcandidate_index,
0531                                                                    long &total_visited_children,
0532                                                                    long &total_visited_leaves, long &total_cut_nodes,
0533                                                                    long &total_stacked_nodes) const
0534   {
0535     total_visited_children = 0;
0536     total_visited_leaves   = 0;
0537     total_cut_nodes        = 0; // Number of nodes put on the stack but skipped when visited due to being too far
0538     total_stacked_nodes    = 0; // Number of nodes put on the stack to visit
0539 
0540     unsigned int stack[BVH_MAX_DEPTH], *ptr = &stack[1];
0541     stack[0] = 0;
0542 
0543     /* Calculate and reuse inverse direction to save on divisions */
0544     Vector3D<Real_t> binvdir(static_cast<Real_t>(1.0) / vecgeom::NonZero(localdir[0]),
0545                              static_cast<Real_t>(1.0) / vecgeom::NonZero(localdir[1]),
0546                              static_cast<Real_t>(1.0) / vecgeom::NonZero(localdir[2]));
0547     Vector3D<Real_t> blocalpoint(static_cast<Real_t>(localpoint[0]), static_cast<Real_t>(localpoint[1]),
0548                                  static_cast<Real_t>(localpoint[2]));
0549     Vector3D<Real_t> blocaldir(static_cast<Real_t>(localdir[0]), static_cast<Real_t>(localdir[1]),
0550                                static_cast<Real_t>(localdir[2]));
0551     Real_t bstep = static_cast<Real_t>(step);
0552 
0553     do {
0554       const unsigned int id = *--ptr; /* pop next node id to be checked from the stack */
0555 
0556       // If the current distance is shorter than the distance to the node we can safely ignore it
0557       Real_t min{vecgeom::InfinityLength<Real_t>()}, max{-vecgeom::InfinityLength<Real_t>()};
0558       fNodes[id].ComputeIntersectionInvDir(blocalpoint, binvdir, min, max);
0559       if (!(min <= max && max >= 0.0 && min < step)) {
0560         total_cut_nodes++;
0561         continue;
0562       }
0563 
0564       if (fNChild[id] >= 0) {
0565         total_visited_leaves++;
0566         /* For leaf nodes, loop over children */
0567         for (int i = 0; i < fNChild[id]; ++i) {
0568           const int prim = fPrimId[fOffset[id] + i];
0569           /* Check AABB first, then the element itself if needed */
0570           Real_t approach;
0571           if (fAABBs[prim].IntersectInvDirApproach(blocalpoint, binvdir, bstep, approach)) {
0572             auto dist = Navigator::CandidateDistanceToIn(
0573                 fRootId, prim, localpoint + static_cast<Real_i>(approach) * localdir, localdir, step);
0574             dist += static_cast<Real_i>(approach);
0575             /* If distance to current child is smaller than current step, update step and hitcandidate */
0576             if (dist < step &&
0577                 !(dist <= vecgeom::kToleranceDist<Real_i> && Navigator::SkipItem(fRootId, prim, last_exited_id))) {
0578               step               = dist;
0579               bstep              = static_cast<Real_t>(dist);
0580               hitcandidate_index = prim;
0581             }
0582           }
0583         }
0584       } else {
0585         const unsigned int childL = 2 * id + 1;
0586         const unsigned int childR = 2 * id + 2;
0587 
0588         /* For internal nodes, check AABBs to know if we need to traverse left and right children */
0589         Real_t tminL = vecgeom::InfinityLength<Real_t>(), tmaxL = -vecgeom::InfinityLength<Real_t>(),
0590                tminR = vecgeom::InfinityLength<Real_t>(), tmaxR = -vecgeom::InfinityLength<Real_t>();
0591 
0592         fNodes[childL].ComputeIntersectionInvDir(blocalpoint, binvdir, tminL, tmaxL);
0593         fNodes[childR].ComputeIntersectionInvDir(blocalpoint, binvdir, tminR, tmaxR);
0594 
0595         const bool traverseL = tminL <= tmaxL && tmaxL >= 0.0 && tminL < step;
0596         const bool traverseR = tminR <= tmaxR && tmaxR >= 0.0 && tminR < step;
0597 
0598         /*
0599          * If both left and right nodes need to be checked, check closest one first.
0600          * This ensures step gets short as fast as possible so we can skip more nodes without checking.
0601          */
0602         if (tminR < tminL) {
0603           if (traverseR) *ptr++ = childR;
0604           if (traverseL) *ptr++ = childL;
0605         } else {
0606           if (traverseL) *ptr++ = childL;
0607           if (traverseR) *ptr++ = childR;
0608         }
0609 
0610         if (traverseR) total_stacked_nodes++;
0611         if (traverseL) total_stacked_nodes++;
0612       }
0613     } while (ptr > stack);
0614   }
0615 
0616 private:
0617   enum class ConstructionAlgorithm : unsigned int;
0618   /**
0619    * Compute internal nodes of the BVH recursively.
0620    * @param[in] id Node id of node to be computed.
0621    * @param[in] first Iterator pointing to the position of this node's first volume in @c fPrimId.
0622    * @param[in] last Iterator pointing to the position of this node's last volume in @c fPrimId.
0623    * @param[in] nodes Number of nodes for this BVH.
0624    * @param[in] constructionAlgorithm Index of the splitting function to use.
0625    *
0626    * @remark This function computes the bounding box of a node, then chooses a split plane and reorders
0627    * the elements of @c fPrimId within the first,last range such that volumes for its left child all come
0628    * before volumes for its right child, then launches itself to compute bounding boxes of each child.
0629    * Recursion stops when all children lie on one side of the splitting plane, or when the current node
0630    * contains only a single child volume.
0631    */
0632   void ComputeNodes(unsigned int id, int *first, int *last, unsigned int nodes, ConstructionAlgorithm);
0633 };
0634 
0635 } // namespace VECGEOM_IMPL_NAMESPACE
0636 } // namespace vecgeom
0637 
0638 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!