surfaces/mask/TriangleMask.h

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