VecGeom/volumes/PlanarPolygonUtilities.h

0001 /// @file PlanarPolygonUtilities.h
0002 /// @brief Geometry predicates shared by PlanarPolygon triangulation users.
0003
0004 #ifndef VECGEOM_PLANAR_POLYGON_UTILITIES_H
0005 #define VECGEOM_PLANAR_POLYGON_UTILITIES_H
0006
0007 #include "VecGeom/base/Global.h"
0008 #include "VecGeom/base/Vector3D.h"
0009 #include "VecGeom/volumes/PlanarPolygon.h"
0010
0011 namespace vecgeom {
0012 inline namespace VECGEOM_IMPL_NAMESPACE {
0013 namespace PlanarPolygonUtilities {
0014
0015 /// @brief Return a base polygon vertex as an XY-plane vector.
0016 VECGEOM_FORCE_INLINE
0017 VECCORE_ATT_HOST_DEVICE
0018 Vector3D<Precision> Vertex2D(PlanarPolygon const &polygon, size_t index)
0019 {
0020   auto const &vertices = polygon.GetVertices();
0021   return Vector3D<Precision>(vertices.x()[index], vertices.y()[index], 0.);
0022 }
0023
0024 /// @brief Check whether a polygon vertex is on the line through two other vertices.
0025 /// @details Uses a length-scaled squared cross product to avoid square roots.
0026 VECGEOM_FORCE_INLINE
0027 VECCORE_ATT_HOST_DEVICE
0028 bool IsSameLine(PlanarPolygon const &polygon, size_t point, size_t first, size_t second)
0029 {
0030   const auto *x        = polygon.GetVertices().x();
0031   const auto *y        = polygon.GetVertices().y();
0032   const Precision dx   = x[second] - x[first];
0033   const Precision dy   = y[second] - y[first];
0034   const Precision tol  = static_cast<Precision>(0.5) * kTolerance;
0035   const Precision tol2 = tol * tol;
0036
0037   const Precision relx = x[point] - x[first];
0038   const Precision rely = y[point] - y[first];
0039   const Precision len2 = dx * dx + dy * dy;
0040   if (len2 <= tol2) return relx * relx + rely * rely <= tol2;
0041
0042   const Precision cross = dx * rely - dy * relx;
0043   return cross * cross <= len2 * tol2;
0044 }
0045
0046 /// @brief Check whether a polygon vertex is on a segment within tolerance.
0047 VECGEOM_FORCE_INLINE
0048 VECCORE_ATT_HOST_DEVICE
0049 bool IsSameLineSegment(PlanarPolygon const &polygon, size_t point, size_t first, size_t second)
0050 {
0051   const auto *x        = polygon.GetVertices().x();
0052   const auto *y        = polygon.GetVertices().y();
0053   const Precision tol  = static_cast<Precision>(0.5) * kTolerance;
0054   const Precision px   = x[point];
0055   const Precision py   = y[point];
0056   const Precision minx = Min(x[first], x[second]) - tol;
0057   const Precision maxx = Max(x[first], x[second]) + tol;
0058   const Precision miny = Min(y[first], y[second]) - tol;
0059   const Precision maxy = Max(y[first], y[second]) + tol;
0060   if (px < minx || px > maxx || py < miny || py > maxy) return false;
0061   return IsSameLine(polygon, point, first, second);
0062 }
0063
0064 /// @brief Check whether two vertices are strictly on the same side of an edge.
0065 VECGEOM_FORCE_INLINE
0066 VECCORE_ATT_HOST_DEVICE
0067 bool IsSameSide(PlanarPolygon const &polygon, size_t lhs, size_t rhs, size_t first, size_t second)
0068 {
0069   const auto *x = polygon.GetVertices().x();
0070   const auto *y = polygon.GetVertices().y();
0071
0072   return ((x[lhs] - x[first]) * (y[second] - y[first]) - (x[second] - x[first]) * (y[lhs] - y[first])) *
0073              ((x[rhs] - x[first]) * (y[second] - y[first]) - (x[second] - x[first]) * (y[rhs] - y[first])) >
0074          0.;
0075 }
0076
0077 /// @brief Check whether a vertex lies inside or on a candidate triangle.
0078 VECGEOM_FORCE_INLINE
0079 VECCORE_ATT_HOST_DEVICE
0080 bool IsPointInsideTriangle(PlanarPolygon const &polygon, size_t first, size_t second, size_t third, size_t point)
0081 {
0082   const auto *x = polygon.GetVertices().x();
0083   const auto *y = polygon.GetVertices().y();
0084
0085   if ((x[point] < x[first] && x[point] < x[second] && x[point] < x[third]) ||
0086       (x[point] > x[first] && x[point] > x[second] && x[point] > x[third]) ||
0087       (y[point] < y[first] && y[point] < y[second] && y[point] < y[third]) ||
0088       (y[point] > y[first] && y[point] > y[second] && y[point] > y[third]))
0089     return false;
0090
0091   const bool inside  = IsSameSide(polygon, point, first, second, third) &&
0092                        IsSameSide(polygon, point, second, first, third) &&
0093                        IsSameSide(polygon, point, third, first, second);
0094   const bool on_edge = IsSameLineSegment(polygon, point, first, second) ||
0095                        IsSameLineSegment(polygon, point, second, third) ||
0096                        IsSameLineSegment(polygon, point, third, first);
0097   return inside || on_edge;
0098 }
0099
0100 /// @brief Check whether three consecutive vertices form a convex clockwise corner.
0101 VECGEOM_FORCE_INLINE
0102 VECCORE_ATT_HOST_DEVICE
0103 bool IsConvexCorner(PlanarPolygon const &polygon, size_t previous, size_t current, size_t next)
0104 {
0105   const auto *x = polygon.GetVertices().x();
0106   const auto *y = polygon.GetVertices().y();
0107   const Precision cross =
0108       (x[current] - x[previous]) * (y[next] - y[current]) - (x[next] - x[current]) * (y[current] - y[previous]);
0109   return cross < 0.;
0110 }
0111
0112 } // namespace PlanarPolygonUtilities
0113 } // namespace VECGEOM_IMPL_NAMESPACE
0114 } // namespace vecgeom
0115
0116 #endif