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
0010 #ifdef VECGEOM_ENABLE_CUDA
0011 #include "VecGeom/backend/cuda/Interface.h"
0012 #endif
0013
0014 #include <algorithm>
0015
0016 namespace vecgeom {
0017 VECGEOM_DEVICE_FORWARD_DECLARE(class AABB;);
0018 VECGEOM_DEVICE_DECLARE_CONV(class, AABB);
0019 inline namespace VECGEOM_IMPL_NAMESPACE {
0020
0021 /**
0022  * @brief Simple class to represent Axis-Aligned Bounding Boxes (AABB).
0023  * @details The AABB is represented internally using the minimum and maximum corners.
0024  */
0025
0026 class AABB {
0027 public:
0028   /** Default constructor. Required to use AABBs as elements in standard containers. */
0029   AABB() = default;
0030   /** Constructor. */
0031   VECCORE_ATT_HOST_DEVICE
0032   AABB(Vector3D<Precision> Min, Vector3D<Precision> Max) : fMin(Min), fMax(Max) {}
0033
0034   /** Returns the minimum coordinates of the AABB. */
0035   VECCORE_ATT_HOST_DEVICE
0036   Vector3D<Precision> Min() const { return fMin; }
0037
0038   /** Returns the maximum coordinates of the AABB. */
0039   VECCORE_ATT_HOST_DEVICE
0040   Vector3D<Precision> Max() const { return fMax; }
0041
0042   /** Returns the center of the AABB. */
0043   VECCORE_ATT_HOST_DEVICE
0044   Vector3D<Precision> Center() const { return 0.5 * (fMax + fMin); }
0045
0046   /** Returns the extents of the AABB along each axis. */
0047   VECCORE_ATT_HOST_DEVICE
0048   Vector3D<Precision> Size() const { return fMax - fMin; }
0049
0050   /** Returns the surface area of the box. */
0051   VECCORE_ATT_HOST_DEVICE
0052   Precision SurfaceArea() const
0053   {
0054     const auto extent = Size();
0055     return 2. * (extent[0]*extent[1] + extent[1]*extent[2] + extent[2]*extent[0]);
0056   }
0057
0058   /** Expand AABB. @param s Amount by which to expand in each direction. */
0059   VECCORE_ATT_HOST_DEVICE
0060   void Expand(Precision s)
0061   {
0062     s *= 0.5;
0063     fMin -= s;
0064     fMax += s;
0065   }
0066
0067   /** Check whether a point is contained by the AABB. */
0068   VECCORE_ATT_HOST_DEVICE
0069   bool Contains(Vector3D<Precision> p) const
0070   {
0071     return p[0] >= fMin[0] && p[0] <= fMax[0] && p[1] >= fMin[1] && p[1] <= fMax[1] && p[2] >= fMin[2] && p[2] <= fMax[2];
0072   }
0073
0074   /**
0075    * Compute a safety margin from a point to AABB's surface.
0076    * The AABB is guaranteed to be further than the safety.
0077    * @param[in] point Input point.
0078    * @remark Returns a negative value if point is inside AABB.
0079    */
0080   VECCORE_ATT_HOST_DEVICE
0081   Precision Safety(Vector3D<Precision> point) const { return ((point - Center()).Abs() - 0.5 * Size()).Max(); }
0082
0083   /**
0084    * Compute distance from a point to AABB's surface along the given direction.
0085    * @param[in] point Starting point for input ray.
0086    * @param[in] direction Direction of the input ray.
0087    * @param[in] step Maximum distance for which an intersection should be reported.
0088    * @remark Returns a negative value if starting point is already inside AABB.
0089    */
0090   VECCORE_ATT_HOST_DEVICE
0091   Precision Distance(Vector3D<Precision> point, Vector3D<Precision> direction) const
0092   {
0093     Precision tmin, tmax;
0094     ComputeIntersection(point, direction, tmin, tmax);
0095     return (tmin < tmax && tmax > 0.0) ? tmin : kInfLength;
0096   }
0097
0098   /**
0099    * Compute distance from a point to AABB's surface along the given direction.
0100    * @param[in] point Starting point for input ray.
0101    * @param[in] invdir Inverse of direction vector of the input ray.
0102    * @param[in] step Maximum distance for which an intersection should be reported.
0103    * @remark Returns a negative value if starting point is already inside AABB.
0104    */
0105   VECCORE_ATT_HOST_DEVICE
0106   Precision DistanceInvDir(Vector3D<Precision> point, Vector3D<Precision> invdir) const
0107   {
0108     Precision tmin, tmax;
0109     ComputeIntersectionInvDir(point, invdir, tmin, tmax);
0110     return (tmin < tmax && tmax > 0.0) ? tmin : kInfLength;
0111   }
0112
0113   /**
0114    * Compute intersection interval with a line, given a point and a direction defining it.
0115    * @param[in] point Starting point on the line.
0116    * @param[in] direction Direction of the line.
0117    * @param tmin[out] Minimum `t` such that `point + t * direction` intersects the AABB.
0118    * @param tmax[out] Maximum `t` such that `point + t * direction` intersects the AABB.
0119    */
0120   VECCORE_ATT_HOST_DEVICE
0121   void ComputeIntersection(Vector3D<Precision> point, Vector3D<Precision> direction, Precision &tmin,
0122                            Precision &tmax) const
0123   {
0124     Vector3D<Precision> invdir(1.0 / NonZero(direction[0]), 1.0 / NonZero(direction[1]), 1.0 / NonZero(direction[2]));
0125     ComputeIntersectionInvDir(point, invdir, tmin, tmax);
0126   }
0127
0128   /**
0129    * Compute intersection interval with a line, given a point and the inverse of the direction vector defining it.
0130    * @param[in] point Starting point on the line.
0131    * @param[in] invdir Inverse of direction vector of the input ray.
0132    * @param tmin[out] Minimum `t` such that `point + t * direction` intersects the AABB.
0133    * @param tmax[out] Maximum `t` such that `point + t * direction` intersects the AABB.
0134    */
0135   VECCORE_ATT_HOST_DEVICE
0136   void ComputeIntersectionInvDir(Vector3D<Precision> point, Vector3D<Precision> invdir, Precision &tmin,
0137                                  Precision &tmax) const
0138   {
0139     auto swap = [](Precision &a, Precision &b) {
0140       Precision tmp = a;
0141       a             = b;
0142       b             = tmp;
0143     };
0144
0145     Vector3D<Precision> t0 = (fMin - point) * invdir;
0146     Vector3D<Precision> t1 = (fMax - point) * invdir;
0147
0148     if (t0[0] > t1[0]) swap(t0[0], t1[0]);
0149     if (t0[1] > t1[1]) swap(t0[1], t1[1]);
0150     if (t0[2] > t1[2]) swap(t0[2], t1[2]);
0151
0152     tmin = t0.Max();
0153     tmax = t1.Min() * (Precision(1.) + 2*kEpsilon); // The 2 epsilon prevent false misses
0154   }
0155
0156   /**
0157    * Check whether the line intersects AABB.
0158    * @param[in] point Starting point on the line.
0159    * @param[in] direction Direction of the line.
0160    */
0161   VECCORE_ATT_HOST_DEVICE
0162   bool Intersect(Vector3D<Precision> point, Vector3D<Precision> direction) const
0163   {
0164     Precision tmin, tmax;
0165     ComputeIntersection(point, direction, tmin, tmax);
0166     return tmin <= tmax && tmax >= 0.0;
0167   }
0168
0169   /**
0170    * Check whether the line intersects AABB.
0171    * @param[in] point Starting point on the line.
0172    * @param[in] invdir Inverse of direction vector of the input ray.
0173    */
0174   VECCORE_ATT_HOST_DEVICE
0175   bool IntersectInvDir(Vector3D<Precision> point, Vector3D<Precision> invdir) const
0176   {
0177     Precision tmin, tmax;
0178     ComputeIntersectionInvDir(point, invdir, tmin, tmax);
0179     return tmin <= tmax && tmax >= 0.0;
0180   }
0181
0182   /**
0183    * Check whether the ray intersects AABB within given step length.
0184    * @param[in] point Starting point for input ray.
0185    * @param[in] direction Direction of the input ray.
0186    * @param[in] step Maximum distance for which an intersection should be reported.
0187    * @remark Does not report an intersection if the AABB lies fully behind the ray.
0188    */
0189   VECCORE_ATT_HOST_DEVICE
0190   bool Intersect(Vector3D<Precision> point, Vector3D<Precision> direction, Precision step) const
0191   {
0192     Precision tmin, tmax;
0193     ComputeIntersection(point, direction, tmin, tmax);
0194     return tmin <= tmax && tmax >= 0.0 && tmin < step;
0195   }
0196
0197   /**
0198    * Check whether the ray intersects AABB within given step length.
0199    * @param[in] point Starting point for input ray.
0200    * @param[in] invdir Inverse of direction vector of the input ray.
0201    * @param[in] step Maximum distance for which an intersection should be reported.
0202    * @remark Does not report an intersection if the AABB lies fully behind the ray.
0203    */
0204   VECCORE_ATT_HOST_DEVICE
0205   bool IntersectInvDir(Vector3D<Precision> point, Vector3D<Precision> invdir, Precision step) const
0206   {
0207     Precision tmin, tmax;
0208     ComputeIntersectionInvDir(point, invdir, tmin, tmax);
0209     return tmin <= tmax && tmax >= 0.0 && tmin < step;
0210   }
0211
0212   /**
0213    * Compute minimum AABB that encloses the two input AABBs, A and B.
0214    */
0215   VECCORE_ATT_HOST_DEVICE
0216   static AABB Union(AABB const &A, AABB const &B)
0217   {
0218     using vecCore::math::Max;
0219     using vecCore::math::Min;
0220     Vector3D<Precision> MinC(Min(A.fMin[0], B.fMin[0]), Min(A.fMin[1], B.fMin[1]), Min(A.fMin[2], B.fMin[2]));
0221     Vector3D<Precision> MaxC(Max(A.fMax[0], B.fMax[0]), Max(A.fMax[1], B.fMax[1]), Max(A.fMax[2], B.fMax[2]));
0222     return {MinC, MaxC};
0223   }
0224
0225 private:
0226   Vector3D<Precision> fMin; ///< Minimum coordinates of the AABB.
0227   Vector3D<Precision> fMax; ///< Maximum coordinates of the AABB.
0228 };
0229
0230 } // namespace VECGEOM_IMPL_NAMESPACE
0231 } // namespace vecgeom
0232
0233 #endif