surfaces/mask/QuadrilateralMask.h

0001 #ifndef VECGEOM_SURFACE_QUADMASK_H
0002 #define VECGEOM_SURFACE_QUADMASK_H
0003
0004 #ifndef QUAD_ACCURATE_SAFETY
0005 #define QUAD_ACCURATE_SAFETY 0
0006 #endif
0007
0008 #include <VecGeom/surfaces/base/CommonTypes.h>
0009
0010 namespace vgbrep {
0011
0012 /// @brief Convex quadrilateral masks on plane surfaces.
0013 /// @tparam Real_t Precision type
0014 template <typename Real_t>
0015 struct QuadrilateralMask {
0016   using value_type = Real_t;
0017   Point2D<Real_t> p_[4] = {Real_t(0)}; ///< 2D coordinates of the vertices.
0018   Point2D<Real_t> n_[4] = {Real_t(0)}; ///< 2D coordinates of the outwards normals to segments
0019
0020   QuadrilateralMask() = default;
0021   /**
0022    * @brief Construct a new Quadrilateral Mask object.
0023    *
0024    * @tparam Real_i Precision type of inputs
0025    * @param x1 is the x coordinate of the lower left corner of the quadrilateral mask.
0026    * @param y1 is the y coordinate of the lower left corner of the quadrilateral mask.
0027    * @details The rest of the points should be entered counter-clockwise
0028    * starting from the lower left corner.
0029    */
0030   template <typename Real_i>
0031   QuadrilateralMask(Real_i x1, Real_i y1, Real_i x2, Real_i y2, Real_i x3, Real_i y3, Real_i x4, Real_i y4)
0032   {
0033     // temporary points and normals in vecgeom Precision for calculation and degeneracy check
0034     Point2D<vecgeom::Precision> p[4] = {vecgeom::Precision(0)};
0035     Point2D<vecgeom::Precision> n[4] = {vecgeom::Precision(0)};
0036     p[0].Set(x1, y1);
0037     p[1].Set(x2, y2);
0038     p[2].Set(x3, y3);
0039     p[3].Set(x4, y4);
0040
0041     // Compute outward normals
0042     for (int i = 0; i < 4; ++i) {
0043       auto j      = (i + 1) % 4;
0044       auto k      = (i + 2) % 4;
0045       auto seg_ij = p[j] - p[i];
0046       auto seg_ik = p[k] - p[i];
0047       VECGEOM_ASSERT(seg_ij.Mag2() >
0048                      vecgeom::kToleranceSquared); // use vecgeom tolerance since p is in vecgeom precision
0049       // normal in XY plane
0050       n[i].Set(seg_ij.y(), -seg_ij.x());
0051       // flip the normal so point k is 'backwards'
0052       if (n[i].Dot(seg_ik) > Real_t(0)) n[i] *= Real_t(-1);
0053       // Normalize normal vector
0054       n[i].Normalize();
0055       // set stored normal n_ in mixed precision
0056       n_[i].Set(static_cast<Real_t>(n[i].x()), static_cast<Real_t>(n[i].y()));
0057     }
0058     // static cast to mixed precision after normal calculation
0059     p_[0].Set(static_cast<Real_t>(x1), static_cast<Real_t>(y1));
0060     p_[1].Set(static_cast<Real_t>(x2), static_cast<Real_t>(y2));
0061     p_[2].Set(static_cast<Real_t>(x3), static_cast<Real_t>(y3));
0062     p_[3].Set(static_cast<Real_t>(x4), static_cast<Real_t>(y4));
0063   }
0064   template <typename Real_i>
0065   QuadrilateralMask(const QuadrilateralMask<Real_i> &other)
0066   {
0067     for (int i = 0; i < 4; ++i) {
0068       p_[i] = Point2D<Real_t>(other.p_[i]);
0069       n_[i] = Point2D<Real_t>(other.n_[i]);
0070     }
0071   }
0072
0073   /// @brief Fills the 3D extent of the quadrilateral
0074   /// @param window Extent window to be filled
0075   /// @param aMin Bottom extent corner
0076   /// @param aMax Top extent corner
0077   void Extent3D(Vector3D<Real_t> &aMin, Vector3D<Real_t> &aMax) const
0078   {
0079     aMin.Set(vecgeom::InfinityLength<Real_t>());
0080     aMax.Set(-vecgeom::InfinityLength<Real_t>());
0081     aMin.z() = aMax.z() = Real_t(0);
0082     for (int i = 0; i < 4; ++i) {
0083       aMin.x() = vecCore::math::Min(aMin.x(), p_[i].x());
0084       aMax.x() = vecCore::math::Max(aMax.x(), p_[i].x());
0085       aMin.y() = vecCore::math::Min(aMin.y(), p_[i].y());
0086       aMax.y() = vecCore::math::Max(aMax.y(), p_[i].y());
0087     }
0088   }
0089
0090   /// @brief Returns the extent of the quadrilateral
0091   /// @param window Extent window to be filled
0092   void GetExtent(WindowMask<Real_t> &window) const
0093   {
0094     window.rangeU.Set(p_[0].x());
0095     window.rangeV.Set(p_[0].y());
0096     for (int i = 0; i < 4; ++i) {
0097       window.rangeU.Set(vecCore::math::Min(window.rangeU[0], p_[i].x()),
0098                         vecCore::math::Max(window.rangeU[1], p_[i].x()));
0099       window.rangeV.Set(vecCore::math::Min(window.rangeV[0], p_[i].y()),
0100                         vecCore::math::Max(window.rangeV[1], p_[i].y()));
0101     }
0102   }
0103
0104   /// @brief Checks if the point is within mask.
0105   /// @details The point is within the quadrilateral if all dot products of point
0106   /// position relative to each vertex and corresponding segment normal are negative.
0107   /// @param local Local coordinates of the point
0108   /// @return true if the point is inside the mask.
0109   VECCORE_ATT_HOST_DEVICE
0110   bool Inside(Vector3D<Real_t> const &local) const
0111   {
0112     Vector2D<Real_t> const local2D(local.x(), local.y());
0113     return (n_[0].Dot(local2D - p_[0]) < vecgeom::kToleranceDist<Real_t> &&
0114             n_[1].Dot(local2D - p_[1]) < vecgeom::kToleranceDist<Real_t> &&
0115             n_[2].Dot(local2D - p_[2]) < vecgeom::kToleranceDist<Real_t> &&
0116             n_[3].Dot(local2D - p_[3]) < vecgeom::kToleranceDist<Real_t>);
0117   }
0118
0119   /// @brief Computes the closest distance from a point in XY plane and the triangle.
0120   /// @param local Local coordinates of the point
0121   /// @return Safety from point to triangle.
0122   VECCORE_ATT_HOST_DEVICE
0123   Real_t Safety(Vector3D<Real_t> const &local, Real_t safetySurf, bool &valid) const
0124   {
0125     valid = true;
0126     Vector2D<Real_t> const local2D(local.x(), local.y());
0127     Real_t safety = safetySurf;
0128     bool withinBound[4];
0129     for (int i = 0; i < 4; ++i) {
0130       withinBound[i] = n_[i].Dot(local2D - p_[i]) <= 0;
0131     }
0132     if (withinBound[0] && withinBound[1] && withinBound[2] && withinBound[3]) return safetySurf;
0133
0134     Real_t dseg_squared = vecgeom::InfinityLength<Real_t>();
0135     for (int i = 0; i < 4; ++i) {
0136       if (!withinBound[i])
0137         dseg_squared = vecCore::math::Min(dseg_squared, DistanceToSegmentSquared(local2D, p_[i], p_[(i + 1) % 4]));
0138     }
0139     safety = vecCore::math::Sqrt(dseg_squared + safetySurf * safetySurf);
0140
0141     return safety;
0142   }
0143
0144   /// @brief Safe distance from a point assumed inside the quadrilateral
0145   /// @details Used on host only for frame checks
0146   /// @param local Projected point in local coordinates
0147   /// @return Safe distance
0148   Real_t SafetyInside(Vector3D<Real_t> const &local) const
0149   {
0150     if (Inside(local) == false) return Real_t(0);
0151     Vector2D<Real_t> const local2D(local.x(), local.y());
0152     Real_t safety_squared = vecgeom::InfinityLength<Real_t>();
0153     for (int i = 0; i < 4; ++i) {
0154       safety_squared = vecCore::math::Min(safety_squared, DistanceToSegmentSquared(local2D, p_[i], p_[(i + 1) % 4]));
0155     }
0156     return vecCore::math::Sqrt(safety_squared);
0157   }
0158
0159   /// @brief Safe distance from a point assumed outside the quadrilateral
0160   /// @details Used on host only for frame checks
0161   /// @param local Projected point in local coordinates
0162   /// @return Safe distance
0163   Real_t SafetyOutside(Vector3D<Real_t> const &local) const
0164   {
0165     if (Inside(local) == true) return Real_t(0);
0166     Vector2D<Real_t> const local2D(local.x(), local.y());
0167     Real_t safety_squared = vecgeom::InfinityLength<Real_t>();
0168     for (int i = 0; i < 3; ++i) {
0169       safety_squared = vecCore::math::Min(safety_squared, DistanceToSegmentSquared(local2D, p_[i], p_[(i + 1) % 4]));
0170     }
0171     return vecCore::math::Sqrt(safety_squared);
0172   }
0173 };
0174
0175 } // namespace vgbrep
0176
0177 #endif