volumes/utilities/ClippedExtent.h

0001 // This file is part of VecGeom and is distributed under the
0002 // conditions in the file LICENSE.txt in the top directory.
0003 // For the full list of authors see CONTRIBUTORS.txt and `git log`.
0004
0005 /// \brief Utilities for host-side clipped extent queries.
0006 /// \file volumes/utilities/ClippedExtent.h
0007
0008 #ifndef VECGEOM_VOLUMES_UTILITIES_CLIPPEDEXTENT_H_
0009 #define VECGEOM_VOLUMES_UTILITIES_CLIPPEDEXTENT_H_
0010
0011 #include "VecGeom/base/Global.h"
0012 #include "VecGeom/base/Transformation3D.h"
0013 #include "VecGeom/base/Vector3D.h"
0014
0015 namespace vecgeom {
0016
0017 /// @brief Coordinate axis selector for clipped extent queries.
0018 enum class Axis { kX = 0, kY = 1, kZ = 2 };
0019
0020 /// @brief Convert an axis selector to a `Vector3D` component index.
0021 /// @param axis Axis selector.
0022 /// @return Component index in `[0, 2]`.
0023 inline int AxisIndex(Axis axis) { return static_cast<int>(axis); }
0024
0025 /// @brief Convert a component index to an axis selector.
0026 /// @param axis Component index in `[0, 2]`.
0027 /// @return Axis selector with the same numeric value.
0028 inline Axis AxisFromIndex(int axis) { return static_cast<Axis>(axis); }
0029
0030 /// @brief Return the lower-case coordinate name of an axis.
0031 /// @param axis Axis selector.
0032 /// @return `"x"`, `"y"`, or `"z"`.
0033 inline const char *AxisName(Axis axis) { return axis == Axis::kX ? "x" : (axis == Axis::kY ? "y" : "z"); }
0034
0035 /// @brief One-dimensional conservative extent interval.
0036 struct ExtentInterval {
0037   /// True when `min` and `max` describe a non-empty interval.
0038   bool exists = false;
0039   /// Lower interval bound when `exists` is true.
0040   Precision min = kInfLength;
0041   /// Upper interval bound when `exists` is true.
0042   Precision max = -kInfLength;
0043
0044   ExtentInterval() = default;
0045
0046   /// @brief Construct an interval from explicit bounds.
0047   /// @param min_value Lower interval bound.
0048   /// @param max_value Upper interval bound.
0049   ExtentInterval(Precision min_value, Precision max_value)
0050       : exists(min_value <= max_value), min(min_value), max(max_value)
0051   {
0052   }
0053 };
0054
0055 /// @brief Conservative axis-aligned box returned by clipped extent queries.
0056 struct ClippedBox {
0057   /// True when `min` and `max` describe a non-empty box.
0058   bool exists = false;
0059   /// Lower box corner when `exists` is true.
0060   Vector3D<Precision> min{kInfLength, kInfLength, kInfLength};
0061   /// Upper box corner when `exists` is true.
0062   Vector3D<Precision> max{-kInfLength, -kInfLength, -kInfLength};
0063
0064   ClippedBox() = default;
0065
0066   /// @brief Construct a box from explicit lower and upper corners.
0067   /// @param min_value Lower box corner.
0068   /// @param max_value Upper box corner.
0069   ClippedBox(Vector3D<Precision> const &min_value, Vector3D<Precision> const &max_value)
0070       : exists(min_value.x() <= max_value.x() && min_value.y() <= max_value.y() && min_value.z() <= max_value.z()),
0071         min(min_value), max(max_value)
0072   {
0073   }
0074
0075   /// @brief Return this box projected onto one axis.
0076   /// @param axis Projection axis.
0077   /// @return Existing interval when this box exists, otherwise an empty interval.
0078   ExtentInterval Interval(Axis axis) const
0079   {
0080     const int index = AxisIndex(axis);
0081     return exists ? ExtentInterval(min[index], max[index]) : ExtentInterval();
0082   }
0083
0084   /// @brief Expand this box to include a point.
0085   /// @param point Point to include.
0086   void Include(Vector3D<Precision> const &point)
0087   {
0088     if (!exists) {
0089       min    = point;
0090       max    = point;
0091       exists = true;
0092       return;
0093     }
0094     for (int axis = 0; axis < 3; ++axis) {
0095       if (point[axis] < min[axis]) min[axis] = point[axis];
0096       if (point[axis] > max[axis]) max[axis] = point[axis];
0097     }
0098   }
0099
0100   /// @brief Expand this box to include another box.
0101   /// @param other Box to include; empty boxes are ignored.
0102   void Include(ClippedBox const &other)
0103   {
0104     if (!other.exists) return;
0105     Include(other.min);
0106     Include(other.max);
0107   }
0108
0109   /// @brief Enlarge all existing box bounds by a symmetric margin.
0110   /// @param margin Non-negative margin to subtract from `min` and add to `max`.
0111   void Inflate(Precision margin)
0112   {
0113     if (!exists || margin <= Precision(0.)) return;
0114     min -= Vector3D<Precision>(margin);
0115     max += Vector3D<Precision>(margin);
0116   }
0117
0118   /// @brief Test whether a point is inside the box within tolerance.
0119   /// @param point Point in the same frame as the box.
0120   /// @param tolerance Symmetric tolerance applied to all faces.
0121   /// @return True when the box exists and contains the point within tolerance.
0122   bool Contains(Vector3D<Precision> const &point, Precision tolerance = Precision(0.)) const
0123   {
0124     if (!exists) return false;
0125     return point.x() >= min.x() - tolerance && point.x() <= max.x() + tolerance && point.y() >= min.y() - tolerance &&
0126            point.y() <= max.y() + tolerance && point.z() >= min.z() - tolerance && point.z() <= max.z() + tolerance;
0127   }
0128 };
0129
0130 /// @brief Axis-aligned mother-frame clipping limits for extent queries.
0131 class VoxelLimits {
0132 public:
0133   /// @brief Construct unlimited, non-empty clipping limits.
0134   VoxelLimits()
0135   {
0136     for (int axis = 0; axis < 3; ++axis) {
0137       fMin[axis]     = -kInfLength;
0138       fMax[axis]     = kInfLength;
0139       fLimited[axis] = false;
0140     }
0141   }
0142
0143   /// @brief Intersect the current limits with one axis interval.
0144   /// @param axis Axis to restrict.
0145   /// @param min_value Lower accepted coordinate on `axis`.
0146   /// @param max_value Upper accepted coordinate on `axis`.
0147   void AddLimit(Axis axis, Precision min_value, Precision max_value)
0148   {
0149     const int index = AxisIndex(axis);
0150     if (min_value > max_value) {
0151       fEmpty = true;
0152       return;
0153     }
0154     if (fLimited[index]) {
0155       if (min_value > fMin[index]) fMin[index] = min_value;
0156       if (max_value < fMax[index]) fMax[index] = max_value;
0157     } else {
0158       fMin[index]     = min_value;
0159       fMax[index]     = max_value;
0160       fLimited[index] = true;
0161     }
0162     if (fMin[index] > fMax[index]) fEmpty = true;
0163   }
0164
0165   /// @brief Query whether one axis has finite limits.
0166   /// @param axis Axis to query.
0167   /// @return True when `axis` has an active finite interval.
0168   bool IsLimited(Axis axis) const { return fLimited[AxisIndex(axis)]; }
0169
0170   /// @brief Query whether any axis has finite limits.
0171   /// @return True when at least one axis is limited.
0172   bool IsLimited() const { return fLimited[0] || fLimited[1] || fLimited[2]; }
0173
0174   /// @brief Query whether previous limit intersections made the voxel empty.
0175   /// @return True when no point can satisfy the limits.
0176   bool IsEmpty() const { return fEmpty; }
0177
0178   /// @brief Lower coordinate limit on one axis.
0179   /// @param axis Axis to query.
0180   /// @return Lower limit, or `-kInfLength` for an unlimited axis.
0181   Precision Min(Axis axis) const { return fMin[AxisIndex(axis)]; }
0182
0183   /// @brief Upper coordinate limit on one axis.
0184   /// @param axis Axis to query.
0185   /// @return Upper limit, or `kInfLength` for an unlimited axis.
0186   Precision Max(Axis axis) const { return fMax[AxisIndex(axis)]; }
0187
0188   /// @brief Test whether a point satisfies all active limits.
0189   /// @param point Mother-frame point.
0190   /// @param tolerance Symmetric tolerance applied to active limits.
0191   /// @return True when the limits are non-empty and the point is inside.
0192   bool Inside(Vector3D<Precision> const &point, Precision tolerance = Precision(0.)) const
0193   {
0194     if (fEmpty) return false;
0195     for (int axis = 0; axis < 3; ++axis) {
0196       if (!fLimited[axis]) continue;
0197       if (point[axis] < fMin[axis] - tolerance || point[axis] > fMax[axis] + tolerance) return false;
0198     }
0199     return true;
0200   }
0201
0202   /// @brief Clip a line segment to the active limits.
0203   /// @param[in,out] start Segment start; replaced by the clipped start on success.
0204   /// @param[in,out] end Segment end; replaced by the clipped end on success.
0205   /// @return True when the segment intersects the limits.
0206   bool ClipSegment(Vector3D<Precision> &start, Vector3D<Precision> &end) const
0207   {
0208     if (fEmpty) return false;
0209     const Vector3D<Precision> original_start = start;
0210     const Vector3D<Precision> delta          = end - start;
0211     Precision tmin                           = 0.;
0212     Precision tmax                           = 1.;
0213
0214     for (int axis = 0; axis < 3; ++axis) {
0215       if (!fLimited[axis]) continue;
0216       const Precision direction = delta[axis];
0217       if (Abs(direction) <= kTolerance) {
0218         if (start[axis] < fMin[axis] || start[axis] > fMax[axis]) return false;
0219         continue;
0220       }
0221       Precision lo = (fMin[axis] - start[axis]) / direction;
0222       Precision hi = (fMax[axis] - start[axis]) / direction;
0223       if (lo > hi) {
0224         const Precision tmp = lo;
0225         lo                  = hi;
0226         hi                  = tmp;
0227       }
0228       if (lo > tmin) tmin = lo;
0229       if (hi < tmax) tmax = hi;
0230       if (tmin > tmax) return false;
0231     }
0232
0233     start = original_start + tmin * delta;
0234     end   = original_start + tmax * delta;
0235     return true;
0236   }
0237
0238 private:
0239   Vector3D<Precision> fMin;
0240   Vector3D<Precision> fMax;
0241   bool fLimited[3] = {false, false, false};
0242   bool fEmpty      = false;
0243 };
0244
0245 /// @brief Intersect a conservative box with voxel limits.
0246 /// @param[in,out] box Box to clip; marked empty if there is no overlap.
0247 /// @param limits Mother-frame limits.
0248 /// @param margin Symmetric margin applied to limits before clipping.
0249 /// @return True when the clipped box remains non-empty.
0250 inline bool IntersectBoxWithLimits(ClippedBox &box, VoxelLimits const &limits, Precision margin = Precision(0.))
0251 {
0252   if (!box.exists || limits.IsEmpty()) {
0253     box.exists = false;
0254     return false;
0255   }
0256   for (int axis = 0; axis < 3; ++axis) {
0257     const Axis clipped_axis = AxisFromIndex(axis);
0258     if (!limits.IsLimited(clipped_axis)) continue;
0259     const Precision lower = limits.Min(clipped_axis) - margin;
0260     const Precision upper = limits.Max(clipped_axis) + margin;
0261     if (box.max[axis] < lower || box.min[axis] > upper) {
0262       box.exists = false;
0263       return false;
0264     }
0265     if (box.min[axis] < lower) box.min[axis] = lower;
0266     if (box.max[axis] > upper) box.max[axis] = upper;
0267   }
0268   return box.exists;
0269 }
0270
0271 /// @brief Transform a local AABB to a padded mother-frame AABB.
0272 /// @param local_min Local lower corner.
0273 /// @param local_max Local upper corner.
0274 /// @param transform Local-to-mother transform used by the existing placement path.
0275 /// @param padding Symmetric padding applied after transformation.
0276 /// @return Conservative mother-frame AABB.
0277 inline ClippedBox TransformLocalAABBToMother(Vector3D<Precision> const &local_min, Vector3D<Precision> const &local_max,
0278                                              Transformation3D const &transform, Precision padding = kAABBPadding)
0279 {
0280   const Vector3D<Precision> delta = local_max - local_min;
0281   ClippedBox box;
0282   for (int x = 0; x <= 1; ++x) {
0283     for (int y = 0; y <= 1; ++y) {
0284       for (int z = 0; z <= 1; ++z) {
0285         Vector3D<Precision> corner(local_min.x() + x * delta.x(), local_min.y() + y * delta.y(),
0286                                    local_min.z() + z * delta.z());
0287         box.Include(transform.InverseTransform(corner));
0288       }
0289     }
0290   }
0291   box.Inflate(padding);
0292   return box;
0293 }
0294
0295 /// @brief Transform a local AABB using a pointer-valued transform.
0296 /// @param local_min Local lower corner.
0297 /// @param local_max Local upper corner.
0298 /// @param transform Local-to-mother transform pointer.
0299 /// @param padding Symmetric padding applied after transformation.
0300 /// @return Conservative mother-frame AABB.
0301 inline ClippedBox TransformLocalAABBToMother(Vector3D<Precision> const &local_min, Vector3D<Precision> const &local_max,
0302                                              Transformation3D const *transform, Precision padding = kAABBPadding)
0303 {
0304   return TransformLocalAABBToMother(local_min, local_max, *transform, padding);
0305 }
0306
0307 } // namespace vecgeom
0308
0309 #endif // VECGEOM_VOLUMES_UTILITIES_CLIPPEDEXTENT_H_