VecGeom/base/BVH.h

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
0015 #include <vector>
0016
0017 namespace vecgeom {
0018 VECGEOM_DEVICE_FORWARD_DECLARE(class BVH;);
0019 VECGEOM_DEVICE_DECLARE_CONV(class, BVH);
0020 inline namespace VECGEOM_IMPL_NAMESPACE {
0021
0022 class LogicalVolume;
0023 class VPlacedVolume;
0024
0025 /**
0026  * @brief Bounding Volume Hierarchy class to represent an axis-aligned bounding volume hierarchy.
0027  * @details BVH instances can be associated with logical volumes to accelerate queries to their child volumes.
0028  */
0029
0030 class BVH {
0031 public:
0032   /** Maximum depth. */
0033   static constexpr int BVH_MAX_DEPTH = 32;
0034
0035   /**
0036    * Constructor.
0037    * @param volume Pointer to logical volume for which the BVH will be created.
0038    * @param depth Depth of the BVH binary tree. Defaults to zero, in which case
0039    * the actual depth will be chosen dynamically based on the number of child volumes.
0040    * When a fixed depth is chosen, it cannot be larger than @p BVH_MAX_DEPTH.
0041    */
0042   BVH(LogicalVolume const &volume, int depth = 0);
0043
0044   /** Destructor. */
0045   ~BVH();
0046
0047 #ifdef VECGEOM_ENABLE_CUDA
0048   /**
0049    * Constructor for GPU. Takes as input pre-constructed BVH buffers.
0050    * @param volume  Reference to logical volume on the device
0051    * @param depth Depth of the BVH binary tree stored in the device buffers.
0052    * @param dPrimId Device buffer with child volume ids
0053    * @param dAABBs  Device buffer with AABBs of child volumes
0054    * @param dOffset Device buffer with offsets in @c dPrimId for first child of each BVH node
0055    * @param dNChild Device buffer with number of children for each BVH node
0056    * @param dNodes AABBs of BVH nodes
0057    */
0058   VECCORE_ATT_DEVICE
0059   BVH(LogicalVolume const *volume, int depth, int *dPrimId, AABB *dAABBs, int *dOffset, int *NChild, AABB *dNodes);
0060 #endif
0061
0062 #ifdef VECGEOM_CUDA_INTERFACE
0063   /** Copy and construct an instance of this BVH on the device, at the device address @p addr. */
0064   DevicePtr<cuda::BVH> CopyToGpu(void *addr) const;
0065 #endif
0066
0067   /** Print a summary of BVH contents */
0068   VECCORE_ATT_HOST_DEVICE
0069   void Print(bool verbose = false) const;
0070
0071   /**
0072    * Check ray defined by <tt>localpoint + t * localdir</tt> for intersections with children
0073    * of the logical volume associated with the BVH, and within a maximum distance of @p step
0074    * along the ray, while ignoring the @p last volume.
0075    * @param[in] localpoint Point in the local coordinates of the logical volume.
0076    * @param[in] localdir Direction in the local coordinates of the logical volume.
0077    * @param[in,out] step Maximum step distance for which intersections should be considered.
0078    * @param[in] last Last volume. This volume is ignored when reporting intersections.
0079    * @param[out] hitcandidate Pointer to volume for which closest intersection was found.
0080    */
0081   /*
0082    * BVH::ComputeDaughterIntersections() computes the intersection of a ray against all children of
0083    * the logical volume. A stack is kept of the node ids that need to be checked. It needs to be at
0084    * most as deep as the binary tree itself because we always first pop the current node, and then
0085    * add at most the two children. For example, for depth two, we pop the root node, then at most we
0086    * add both of its leaves onto the stack to be checked. We initialize ptr with &stack[1] such that
0087    * when we pop the first time as we enter the loop, the position we read from is the first position
0088    * of the stack, which contains the id 0 for the root node. When we pop the stack such that ptr
0089    * points before &stack[0], it means we've checked all we needed and the loop can be terminated.
0090    * In order to determine if a node of the tree is internal or not, we check if the node id of its
0091    * left child is past the end of the array (in which case we know we are at the maximum depth), or
0092    * if the sum of children in both leaves is the same as in the current node, as for leaf nodes, the
0093    * sum of children in the left+right child nodes will be less than for the current node.
0094    */
0095   VECCORE_ATT_HOST_DEVICE
0096   void CheckDaughterIntersections(Vector3D<Precision> localpoint, Vector3D<Precision> localdir, Precision &step,
0097                                   VPlacedVolume const *last, VPlacedVolume const *&hitcandidate) const
0098   {
0099     unsigned int stack[BVH_MAX_DEPTH], *ptr = &stack[1];
0100     stack[0] = 0;
0101
0102     /* Calculate and reuse inverse direction to save on divisions */
0103     Vector3D<Precision> invdir(1.0 / NonZero(localdir[0]), 1.0 / NonZero(localdir[1]), 1.0 / NonZero(localdir[2]));
0104
0105     do {
0106       const unsigned int id = *--ptr; /* pop next node id to be checked from the stack */
0107
0108       if (fNChild[id] >= 0) {
0109         /* For leaf nodes, loop over children */
0110         for (int i = 0; i < fNChild[id]; ++i) {
0111           const int prim = fPrimId[fOffset[id] + i];
0112           /* Check AABB first, then the volume itself if needed */
0113           if (fAABBs[prim].IntersectInvDir(localpoint, invdir, step)) {
0114             const auto vol  = fLV.GetDaughters()[prim];
0115             const auto dist = vol->DistanceToIn(localpoint, localdir, step);
0116             /* If distance to current child is smaller than current step, update step and hitcandidate */
0117             if (dist < step && !(dist <= 0.0 && vol == last)) step = dist, hitcandidate = vol;
0118           }
0119         }
0120       } else {
0121         const unsigned int childL = 2 * id + 1;
0122         const unsigned int childR = 2 * id + 2;
0123
0124         /* For internal nodes, check AABBs to know if we need to traverse left and right children */
0125         Precision tminL = kInfLength, tmaxL = -kInfLength, tminR = kInfLength, tmaxR = -kInfLength;
0126
0127         fNodes[childL].ComputeIntersectionInvDir(localpoint, invdir, tminL, tmaxL);
0128         fNodes[childR].ComputeIntersectionInvDir(localpoint, invdir, tminR, tmaxR);
0129
0130         const bool traverseL = tminL <= tmaxL && tmaxL >= 0.0 && tminL < step;
0131         const bool traverseR = tminR <= tmaxR && tmaxR >= 0.0 && tminR < step;
0132
0133         /*
0134          * If both left and right nodes need to be checked, check closest one first.
0135          * This ensures step gets short as fast as possible so we can skip more nodes without checking.
0136          */
0137         if (tminR < tminL) {
0138           if (traverseR) *ptr++ = childR;
0139           if (traverseL) *ptr++ = childL;
0140         } else {
0141           if (traverseL) *ptr++ = childL;
0142           if (traverseR) *ptr++ = childR;
0143         }
0144       }
0145     } while (ptr > stack);
0146   }
0147
0148   /**
0149    * Check ray defined by <tt>localpoint + t * localdir</tt> for intersections with bounding
0150    * boxes of children of the logical volume associated with the BVH, and within a maximum
0151    * distance of @p step along the ray. Returns the distance to the first crossed box.
0152    * @param[in] localpoint Point in the local coordinates of the logical volume.
0153    * @param[in] localdir Direction in the local coordinates of the logical volume.
0154    * @param[in,out] step Maximum step distance for which intersections should be considered.
0155    * @param[in] last Last volume. This volume is ignored when reporting intersections.
0156    */
0157   VECCORE_ATT_HOST_DEVICE
0158   void ApproachNextDaughter(Vector3D<Precision> localpoint, Vector3D<Precision> localdir, Precision &step,
0159                             VPlacedVolume const *last) const;
0160
0161   /**
0162    * Compute safety against children of the logical volume associated with the BVH.
0163    * @param[in] localpoint Point in the local coordinates of the logical volume.
0164    * @param[in] safety Maximum safety. Volumes further than this are not checked.
0165    * @returns Minimum between safety to the closest child of logical volume and input @p safety.
0166    */
0167   /*
0168    * BVH::ComputeSafety is very similar to the method above regarding traversal of the tree, but it
0169    * computes only the safety instead of the intersection using a ray, so the logic is a bit simpler.
0170    */
0171   VECCORE_ATT_HOST_DEVICE
0172   Precision ComputeSafety(Vector3D<Precision> localpoint, Precision safety) const
0173   {
0174     unsigned int stack[BVH_MAX_DEPTH], *ptr = &stack[1];
0175     stack[0] = 0;
0176
0177     do {
0178       const unsigned int id = *--ptr;
0179
0180       if (fNChild[id] >= 0) {
0181         for (int i = 0; i < fNChild[id]; ++i) {
0182           const int prim = fPrimId[fOffset[id] + i];
0183           if (fAABBs[prim].Safety(localpoint) < safety) {
0184             const Precision dist = fLV.GetDaughters()[prim]->SafetyToIn(localpoint);
0185             if (dist < safety) safety = dist;
0186           }
0187         }
0188       } else {
0189         const unsigned int childL = 2 * id + 1;
0190         const unsigned int childR = 2 * id + 2;
0191
0192         const Precision safetyL = fNodes[childL].Safety(localpoint);
0193         const Precision safetyR = fNodes[childR].Safety(localpoint);
0194
0195         const bool traverseL = safetyL < safety;
0196         const bool traverseR = safetyR < safety;
0197
0198         if (safetyR < safetyL) {
0199           if (traverseR) *ptr++ = childR;
0200           if (traverseL) *ptr++ = childL;
0201         } else {
0202           if (traverseL) *ptr++ = childL;
0203           if (traverseR) *ptr++ = childR;
0204         }
0205       }
0206     } while (ptr > stack);
0207
0208     return safety;
0209   }
0210
0211   /**
0212    * Find child volume inside which the given point @p localpoint is located.
0213    * @param[in] localpoint Point in the local coordinates of the logical volume.
0214    * @param[out] daughterpvol Placed volume in which @p localpoint is contained
0215    * @param[out] daughterlocalpoint Point in the local coordinates of @p daughterpvol
0216    * @returns Whether @p localpoint falls within a child volume of @p lvol.
0217    */
0218   VECCORE_ATT_HOST_DEVICE
0219   bool LevelLocate(Vector3D<Precision> const &localpoint, VPlacedVolume const *&pvol,
0220                    Vector3D<Precision> &daughterlocalpoint) const
0221   {
0222     VPlacedVolume const *exclvol = nullptr;
0223     return LevelLocate(exclvol, localpoint, pvol, daughterlocalpoint);
0224   }
0225
0226   /**
0227    * Find child volume inside which the given point @p localpoint is located.
0228    * @param[in] localpoint Point in the local coordinates of the logical volume.
0229    * @param[out] outstate Navigation state. Gets updated if point is relocated to another volume.
0230    * @param[out] daughterlocalpoint Point in the local coordinates of newly located volume.
0231    * @returns Whether @p localpoint falls within a child volume of @p lvol.
0232    */
0233   VECCORE_ATT_HOST_DEVICE
0234   bool LevelLocate(Vector3D<Precision> const &localpoint, NavigationState &state,
0235                    Vector3D<Precision> &daughterlocalpoint) const
0236   {
0237     VPlacedVolume const *exclvol = nullptr;
0238     VPlacedVolume const *pvol    = nullptr;
0239     bool Result                  = LevelLocate(exclvol, localpoint, pvol, daughterlocalpoint);
0240     if (Result) {
0241       state.Push(pvol);
0242     }
0243     return Result;
0244   }
0245
0246   /**
0247    * Find child volume inside which the given point @p localpoint is located.
0248    * @param[in] exclvol Placed volume that should be ignored.
0249    * @param[in] localpoint Point in the local coordinates of the logical volume.
0250    * @param[out] pvol Placed volume in which @p localpoint is contained
0251    * @param[out] daughterlocalpoint Point in the local coordinates of @p daughterpvol
0252    * @returns Whether @p localpoint falls within a child volume of @p lvol.
0253    */
0254   VECCORE_ATT_HOST_DEVICE
0255   bool LevelLocate(VPlacedVolume const *exclvol, Vector3D<Precision> const &localpoint, VPlacedVolume const *&pvol,
0256                    Vector3D<Precision> &daughterlocalpoint) const
0257   {
0258     unsigned int stack[BVH_MAX_DEPTH], *ptr = &stack[1];
0259     stack[0] = 0;
0260
0261     do {
0262       const unsigned int id = *--ptr;
0263
0264       if (fNChild[id] >= 0) {
0265         for (int i = 0; i < fNChild[id]; ++i) {
0266           const int prim = fPrimId[fOffset[id] + i];
0267           if (fAABBs[prim].Contains(localpoint)) {
0268             const auto vol = fLV.GetDaughters()[prim];
0269             if (vol != exclvol && vol->Contains(localpoint, daughterlocalpoint)) {
0270               pvol = vol;
0271               return true;
0272             }
0273           }
0274         }
0275       } else {
0276         const unsigned int childL = 2 * id + 1;
0277         if (fNodes[childL].Contains(localpoint)) *ptr++ = childL;
0278
0279         const unsigned int childR = 2 * id + 2;
0280         if (fNodes[childR].Contains(localpoint)) *ptr++ = childR;
0281       }
0282     } while (ptr > stack);
0283
0284     return false;
0285   }
0286
0287   /**
0288    * Find child volume inside which the given point @p localpoint is located.
0289    * @param[in] exclvol Placed volume that should be ignored.
0290    * @param[in] localpoint Point in the local coordinates of the logical volume.
0291    * @param[in] localdir Direction in the local coordinates of the logical volume.
0292    * @param[out] pvol Placed volume in which @p localpoint is contained
0293    * @param[out] daughterlocalpoint Point in the local coordinates of @p daughterpvol
0294    * @returns Whether @p localpoint falls within a child volume of @p lvol.
0295    */
0296   VECCORE_ATT_HOST_DEVICE
0297   bool LevelLocate(VPlacedVolume const *exclvol, Vector3D<Precision> const &localpoint,
0298                    Vector3D<Precision> const &localdirection, VPlacedVolume const *&pvol,
0299                    Vector3D<Precision> &daughterlocalpoint) const
0300   {
0301     unsigned int stack[BVH_MAX_DEPTH], *ptr = &stack[1];
0302     stack[0] = 0;
0303
0304     do {
0305       const unsigned int id = *--ptr;
0306
0307       if (fNChild[id] >= 0) {
0308         for (int i = 0; i < fNChild[id]; ++i) {
0309           const int prim = fPrimId[fOffset[id] + i];
0310           if (fAABBs[prim].Contains(localpoint)) {
0311             const auto v = fLV.GetDaughters()[prim];
0312
0313             if (v == exclvol) continue;
0314
0315             const auto T = v->GetTransformation();
0316             const auto u = v->GetUnplacedVolume();
0317             const auto p = T->Transform(localpoint);
0318
0319             auto Entering = [&]() {
0320               const Vector3D<Precision> dir = T->TransformDirection(localdirection);
0321               Vector3D<Precision> normal;
0322               u->Normal(p, normal);
0323               return Vector3D<Precision>::Dot(normal, dir) < 0.0;
0324             };
0325
0326             const auto inside = u->Inside(p);
0327
0328             if (inside == kInside || (inside == kSurface && Entering())) {
0329               pvol = v, daughterlocalpoint = p;
0330               return true;
0331             }
0332           }
0333         }
0334       } else {
0335         const unsigned int childL = 2 * id + 1;
0336         if (fNodes[childL].Contains(localpoint)) *ptr++ = childL;
0337
0338         const unsigned int childR = 2 * id + 2;
0339         if (fNodes[childR].Contains(localpoint)) *ptr++ = childR;
0340       }
0341     } while (ptr > stack);
0342
0343     return false;
0344   }
0345
0346 private:
0347   enum class ConstructionAlgorithm : unsigned int;
0348   /**
0349    * Compute internal nodes of the BVH recursively.
0350    * @param[in] id Node id of node to be computed.
0351    * @param[in] first Iterator pointing to the position of this node's first volume in @c fPrimId.
0352    * @param[in] last Iterator pointing to the position of this node's last volume in @c fPrimId.
0353    * @param[in] nodes Number of nodes for this BVH.
0354    * @param[in] constructionAlgorithm Index of the splitting function to use.
0355    *
0356    * @remark This function computes the bounding box of a node, then chooses a split plane and reorders
0357    * the elements of @c fPrimId within the first,last range such that volumes for its left child all come
0358    * before volumes for its right child, then launches itself to compute bounding boxes of each child.
0359    * Recursion stops when all children lie on one side of the splitting plane, or when the current node
0360    * contains only a single child volume.
0361    */
0362   void ComputeNodes(unsigned int id, int *first, int *last, unsigned int nodes, ConstructionAlgorithm);
0363
0364   LogicalVolume const &fLV; ///< Logical volume this BVH was constructed for
0365   int *fPrimId;             ///< Child volume ids for each BVH node
0366   int *fOffset;             ///< Offset in @c fPrimId for first child of each BVH node
0367   int *fNChild;             ///< Number of children for each BVH node
0368   AABB *fNodes;             ///< AABBs of BVH nodes
0369   AABB *fAABBs;             ///< AABBs of children of logical volume @c fLV
0370   int fDepth;               ///< Depth of the BVH
0371 };
0372
0373 } // namespace VECGEOM_IMPL_NAMESPACE
0374 } // namespace vecgeom
0375
0376 #endif