VecGeom/base/AABB.h

0001 /// \file AABB.h
0002 /// \author Guilherme Amadio
0003
0004 #ifndef VECGEOM_BASE_AABB_H_
0005 #define VECGEOM_BASE_AABB_H_
0006
0007 #include "VecGeom/base/Config.h"
0008 #include "VecGeom/base/Vector3D.h"
0009 #include "VecGeom/volumes/kernel/GenericKernels.h"
0010
0011 #ifdef VECGEOM_ENABLE_CUDA
0012 #include "VecGeom/backend/cuda/Interface.h"
0013 #endif
0014
0015 #include <algorithm>
0016
0017 namespace vecgeom {
0018 namespace cuda {
0019 template <typename Real_t>
0020 class AABB;
0021 }
0022 VECGEOM_DEVICE_DECLARE_CONV_TEMPLATE(class, AABB, typename);
0023 inline namespace VECGEOM_IMPL_NAMESPACE {
0024
0025 /**
0026  * @brief Simple class to represent Axis-Aligned Bounding Boxes (AABB).
0027  * @details The AABB is represented internally using the minimum and maximum corners.
0028  */
0029
0030 template <typename Real_t>
0031 class AABB {
0032 private:
0033   Vector3D<Real_t> fMin; ///< Minimum coordinates of the AABB.
0034   Vector3D<Real_t> fMax; ///< Maximum coordinates of the AABB.
0035 public:
0036   /** Default constructor. Required to use AABBs as elements in standard containers. */
0037   AABB() = default;
0038   /** Constructor. */
0039   VECCORE_ATT_HOST_DEVICE
0040   AABB(Vector3D<Real_t> Min, Vector3D<Real_t> Max) : fMin(Min), fMax(Max) {}
0041
0042   /** Returns the minimum coordinates of the AABB. */
0043   VECCORE_ATT_HOST_DEVICE
0044   Vector3D<Real_t> Min() const { return fMin; }
0045
0046   /** Returns the maximum coordinates of the AABB. */
0047   VECCORE_ATT_HOST_DEVICE
0048   Vector3D<Real_t> Max() const { return fMax; }
0049
0050   /** Returns the center of the AABB. */
0051   VECCORE_ATT_HOST_DEVICE
0052   Vector3D<Real_t> Center() const { return static_cast<Real_t>(0.5) * (fMax + fMin); }
0053
0054   /** Returns the extents of the AABB along each axis. */
0055   VECCORE_ATT_HOST_DEVICE
0056   Vector3D<Real_t> Size() const { return fMax - fMin; }
0057
0058   /** Returns the surface area of the box. */
0059   VECCORE_ATT_HOST_DEVICE
0060   Real_t SurfaceArea() const
0061   {
0062     const auto extent = Size();
0063     return static_cast<Real_t>(2.) * (extent[0] * extent[1] + extent[1] * extent[2] + extent[2] * extent[0]);
0064   }
0065
0066   /** Expand AABB. @param s Amount by which to expand in each direction. */
0067   VECCORE_ATT_HOST_DEVICE
0068   void Expand(Real_t s)
0069   {
0070     s *= static_cast<Real_t>(0.5);
0071     fMin -= s;
0072     fMax += s;
0073   }
0074
0075   /** Check whether a point is contained by the AABB. */
0076   VECCORE_ATT_HOST_DEVICE
0077   bool Contains(Vector3D<Real_t> p) const
0078   {
0079     return p[0] >= fMin[0] && p[0] <= fMax[0] && p[1] >= fMin[1] && p[1] <= fMax[1] && p[2] >= fMin[2] &&
0080            p[2] <= fMax[2];
0081   }
0082
0083   /**
0084    * Compute a safety margin from a point to AABB's surface.
0085    * The AABB is guaranteed to be further than the safety.
0086    * @param[in] point Input point.
0087    * @remark Returns a negative value if point is inside AABB.
0088    */
0089   VECCORE_ATT_HOST_DEVICE
0090   Real_t Safety(Vector3D<Real_t> point) const
0091   {
0092     return ((point - Center()).Abs() - static_cast<Real_t>(0.5) * Size()).Max();
0093   }
0094
0095   /**
0096    * Compute distance from a point to AABB's surface along the given direction.
0097    * @param[in] point Starting point for input ray.
0098    * @param[in] direction Direction of the input ray.
0099    * @param[in] step Maximum distance for which an intersection should be reported.
0100    * @remark Returns a negative value if starting point is already inside AABB.
0101    */
0102   VECCORE_ATT_HOST_DEVICE
0103   Real_t Distance(Vector3D<Real_t> point, Vector3D<Real_t> direction) const
0104   {
0105     Real_t tmin, tmax;
0106     ComputeIntersection(point, direction, tmin, tmax);
0107     return (tmin < tmax && tmax > static_cast<Real_t>(0.0)) ? tmin : vecgeom::InfinityLength<Real_t>();
0108   }
0109
0110   /**
0111    * Compute distance from a point to AABB's surface along the given direction.
0112    * @param[in] point Starting point for input ray.
0113    * @param[in] invdir Inverse of direction vector of the input ray.
0114    * @param[in] step Maximum distance for which an intersection should be reported.
0115    * @remark Returns a negative value if starting point is already inside AABB.
0116    */
0117   VECCORE_ATT_HOST_DEVICE
0118   Real_t DistanceInvDir(Vector3D<Real_t> point, Vector3D<Real_t> invdir) const
0119   {
0120     Real_t tmin, tmax;
0121     ComputeIntersectionInvDir(point, invdir, tmin, tmax);
0122     return (tmin < tmax && tmax > static_cast<Real_t>(0.0)) ? tmin : vecgeom::InfinityLength<Real_t>();
0123   }
0124
0125   /**
0126    * Compute intersection interval with a line, given a point and a direction defining it.
0127    * @param[in] point Starting point on the line.
0128    * @param[in] direction Direction of the line.
0129    * @param tmin[out] Minimum `t` such that `point + t * direction` intersects the AABB.
0130    * @param tmax[out] Maximum `t` such that `point + t * direction` intersects the AABB.
0131    */
0132   VECCORE_ATT_HOST_DEVICE
0133   void ComputeIntersection(Vector3D<Real_t> point, Vector3D<Real_t> direction, Real_t &tmin, Real_t &tmax) const
0134   {
0135     Vector3D<Real_t> invdir(1.0 / NonZero(direction[0]), 1.0 / NonZero(direction[1]), 1.0 / NonZero(direction[2]));
0136     ComputeIntersectionInvDir(point, invdir, tmin, tmax);
0137   }
0138
0139   /**
0140    * Compute intersection interval with a line, given a point and the inverse of the direction vector defining it.
0141    * @param[in] point Starting point on the line.
0142    * @param[in] invdir Inverse of direction vector of the input ray.
0143    * @param tmin[out] Minimum `t` such that `point + t * direction` intersects the AABB.
0144    * @param tmax[out] Maximum `t` such that `point + t * direction` intersects the AABB.
0145    */
0146   VECCORE_ATT_HOST_DEVICE
0147   void ComputeIntersectionInvDir(Vector3D<Real_t> point, Vector3D<Real_t> invdir, Real_t &tmin, Real_t &tmax) const
0148   {
0149     auto swap = [](Real_t &a, Real_t &b) {
0150       Real_t tmp = a;
0151       a          = b;
0152       b          = tmp;
0153     };
0154
0155     Vector3D<Real_t> t0 = (fMin - point) * invdir;
0156     Vector3D<Real_t> t1 = (fMax - point) * invdir;
0157
0158     if (t0[0] > t1[0]) swap(t0[0], t1[0]);
0159     if (t0[1] > t1[1]) swap(t0[1], t1[1]);
0160     if (t0[2] > t1[2]) swap(t0[2], t1[2]);
0161
0162     tmin = t0.Max();
0163     tmax = t1.Min() * (static_cast<Real_t>(1.) + vecgeom::kToleranceDist<Real_t>); // The 2 epsilon prevent false misses
0164   }
0165
0166   /**
0167    * Check whether the line intersects AABB.
0168    * @param[in] point Starting point on the line.
0169    * @param[in] direction Direction of the line.
0170    */
0171   VECCORE_ATT_HOST_DEVICE
0172   bool Intersect(Vector3D<Real_t> point, Vector3D<Real_t> direction) const
0173   {
0174     Real_t tmin, tmax;
0175     ComputeIntersection(point, direction, tmin, tmax);
0176     return tmin <= tmax && tmax >= static_cast<Real_t>(0.0);
0177   }
0178
0179   /**
0180    * Check whether the line intersects AABB.
0181    * @param[in] point Starting point on the line.
0182    * @param[in] invdir Inverse of direction vector of the input ray.
0183    */
0184   VECCORE_ATT_HOST_DEVICE
0185   bool IntersectInvDir(Vector3D<Real_t> point, Vector3D<Real_t> invdir) const
0186   {
0187     Real_t tmin, tmax;
0188     ComputeIntersectionInvDir(point, invdir, tmin, tmax);
0189     return tmin <= tmax && tmax >= static_cast<Real_t>(0.0);
0190   }
0191
0192   /**
0193    * Check whether the ray intersects AABB within given step length.
0194    * @param[in] point Starting point for input ray.
0195    * @param[in] direction Direction of the input ray.
0196    * @param[in] step Maximum distance for which an intersection should be reported.
0197    * @remark Does not report an intersection if the AABB lies fully behind the ray.
0198    */
0199   VECCORE_ATT_HOST_DEVICE
0200   bool Intersect(Vector3D<Real_t> point, Vector3D<Real_t> direction, Real_t step) const
0201   {
0202     Real_t tmin, tmax;
0203     ComputeIntersection(point, direction, tmin, tmax);
0204     return tmin <= tmax && tmax >= static_cast<Real_t>(0.0) && tmin < step;
0205   }
0206
0207   /**
0208    * Check whether the ray intersects AABB within given step length.
0209    * @param[in] point Starting point for input ray.
0210    * @param[in] invdir Inverse of direction vector of the input ray.
0211    * @param[in] step Maximum distance for which an intersection should be reported.
0212    * @remark Does not report an intersection if the AABB lies fully behind the ray.
0213    */
0214   VECCORE_ATT_HOST_DEVICE
0215   bool IntersectInvDir(Vector3D<Real_t> point, Vector3D<Real_t> invdir, Real_t step) const
0216   {
0217     Real_t tmin, tmax;
0218     ComputeIntersectionInvDir(point, invdir, tmin, tmax);
0219     bool hit = tmin <= tmax && tmax >= Real_t(0.0);
0220     if (hit && tmin > step) {
0221       // Estimate maximum error of the result to correct the step limit
0222       auto err = vecgeom::kEpsilonT<Real_t> * point.Abs().Max() / tmin;
0223       hit      = tmin < step + err;
0224     }
0225     return hit;
0226   }
0227
0228   /**
0229    * Check whether the ray intersects AABB within given step length. Gives also the approach distance with epsilon
0230    * before the box hit
0231    * @param[in] point Starting point for input ray.
0232    * @param[in] invdir Inverse of direction vector of the input ray.
0233    * @param[in] step Maximum distance for which an intersection should be reported.
0234    * @remark Does not report an intersection if the AABB lies fully behind the ray.
0235    */
0236   VECCORE_ATT_HOST_DEVICE
0237   bool IntersectInvDirApproach(Vector3D<Real_t> point, Vector3D<Real_t> invdir, Real_t step, Real_t &approach) const
0238   {
0239     Real_t tmin, tmax;
0240     approach = Real_t(0.0);
0241     ComputeIntersectionInvDir(point, invdir, tmin, tmax);
0242     if (tmax < Real_t(0.0) || tmin > tmax) return false;
0243
0244     Real_t ulp = ULP<Real_t>(vecCore::math::Max(point.Abs().Max(), tmin));
0245     // Do not approach if distance less than unit
0246     // Overestimate error to 10 ULP (corresponding to 20 roundings in the bad direction)
0247     if (tmin < (step + (Real_t(10.) * ulp)) && tmin > Real_t(1.)) {
0248       // For use with the surface model, it was found that a larger error overestimation is needed in order to avoid
0249       // approaches that overshoot the surface
0250       approach = vecCore::math::Max(tmin - (Real_t(500.) * ulp), Real_t(0.));
0251     }
0252     return true;
0253   }
0254
0255   /**
0256    * Compute minimum AABB that encloses the two input AABBs, A and B.
0257    */
0258   VECCORE_ATT_HOST_DEVICE
0259   static AABB Union(AABB const &A, AABB const &B)
0260   {
0261     using vecCore::math::Max;
0262     using vecCore::math::Min;
0263     Vector3D<Real_t> MinC(Min(A.fMin[0], B.fMin[0]), Min(A.fMin[1], B.fMin[1]), Min(A.fMin[2], B.fMin[2]));
0264     Vector3D<Real_t> MaxC(Max(A.fMax[0], B.fMax[0]), Max(A.fMax[1], B.fMax[1]), Max(A.fMax[2], B.fMax[2]));
0265     return {MinC, MaxC};
0266   }
0267 };
0268
0269 } // namespace VECGEOM_IMPL_NAMESPACE
0270 } // namespace vecgeom
0271
0272 #endif