VecGeom/base/Utils3D.h

0001 ///
0002 /// \file Utils3D.h
0003 /// \author Andrei Gheata (andrei.gheata@cern.ch)
0004
0005 #ifndef VECGEOM_BASE_UTILS3D_H_
0006 #define VECGEOM_BASE_UTILS3D_H_
0007
0008 #include "VecGeom/base/Vector3D.h"
0009 #include "VecGeom/base/Transformation3D.h"
0010
0011 #ifndef VECCORE_CUDA
0012 #include <vector>
0013 #else
0014 #include "VecGeom/base/Vector.h"
0015 #endif
0016
0017 /// A set of 3D geometry utilities
0018 namespace vecgeom {
0019 inline namespace VECGEOM_IMPL_NAMESPACE {
0020
0021 namespace Utils3D {
0022
0023 using Vec_t = Vector3D<Precision>;
0024
0025 template <typename T>
0026 #ifndef VECCORE_CUDA
0027 using vector_t = std::vector<T>;
0028 #else
0029 using vector_t = vecgeom::Vector<T>; // this has problems, like in the initializers Vector v = {...};
0030 #endif
0031
0032 enum EPlaneXing_t { kParallel = 0, kIdentical, kIntersecting };
0033
0034 enum EBodyXing_t { kDisjoint = 0, kTouching, kOverlapping }; // do not change order
0035
0036 /// @brief A basic plane in Hessian normal form
0037 struct Plane {
0038   Vector3D<Precision> fNorm; ///< Unit normal vector to plane
0039   Precision fDist = 0.;      ///< Distance to plane (positive if origin on same side as normal)
0040
0041   Plane() : fNorm() {}
0042
0043   Plane(Vector3D<Precision> const &norm, Precision dist)
0044   {
0045     fNorm = norm;
0046     fDist = dist;
0047   }
0048
0049   /// @brief Transform the plane by a general transformation
0050   void Transform(Transformation3D const &tr);
0051 };
0052
0053 /// Structure to hold line-line intersection results
0054 ///
0055 /// In the case of fIntersect:
0056 ///                                    ^ lA.fPts[1]
0057 ///                                   /
0058 ///                               lA /
0059 ///                                 /
0060 ///                                /              lB
0061 ///         <---------------------X----------------------->
0062 ///  lB.fPts[0]                  /                          lB.fPts[1]
0063 ///                             /
0064 ///                            /    X = lA.fPts[0] + fA * (lA.fPts[1] - lA.fPts[0])
0065 ///                           /     X = lB.fPts[0] + fB * (lB.fPts[1] - lB.fPts[0])
0066 ///               lA.fPts[0] v
0067 ///
0068 ///
0069 ///
0070 /// In the case of fOverlap:
0071 ///
0072 ///                     lB.fPts[0]                       lB.fPts[1]
0073 ///         <------------------<------------>------------->
0074 ///  lA.fPts[0]                              lA.fPts[1]
0075 ///
0076 ///
0077 ///                lB.fPts[0] = lA.fPts[0] + fA * (lA.fPts[1] - lA.fPts[0])
0078 ///                lB.fPts[1] = lA.fPts[0] + fB * (lA.fPts[1] - lA.fPts[0])
0079 ///
0080 /// In the case of fNoIntersect and fParallel, fA and fB have no meaning.
0081
0082 struct LineIntersection {
0083   enum { fParallel, fOverlap, fIntersect, fNoIntersect } fType; ///< Intersection type
0084   Precision fA; ///< Used in calculating the intersection point (see the structure explanation)
0085   Precision fB; ///< Used in calculating the intersection point (see the structure explanation)
0086 };
0087
0088 /// @brief A line segment
0089 struct Line {
0090   Vector3D<Precision> fPts[2]; ///< Points defining the line
0091
0092   /// Computes line-line intersection.
0093   LineIntersection *Intersect(const Line &lB);
0094
0095   /// Indicates whether a point lies on the line segment defined by fPts
0096   bool IsPointOnLine(const Vec_t &p);
0097 };
0098
0099 /// @brief A polygon defined by vertices and normal
0100 /* The list of vertices is a reference to an external array. The used vertex indices have to be defined such
0101    that consecutive segments cross product is on the same side as the normal. */
0102 struct Polygon {
0103   size_t fN       = 0;     ///< Number of vertices
0104   bool fConvex    = false; ///< Convexity
0105   bool fHasNorm   = false; ///< Normal is already supplied
0106   bool fValid     = false; ///< Polygon is not degenerate
0107   Precision fDist = 0.;    ///< Distance to plane in the Hessian form
0108   Vec_t fNorm;             ///< Unit normal vector to plane
0109   vector_t<Vec_t> &fVert;  ///< Global list of vertices shared with other polygons
0110   vector_t<size_t> fInd;   ///< [fN] Indices of vertices
0111   vector_t<Vec_t> fSides;  ///< [fN] Side vectors
0112
0113   /// @brief Constructor taking the number of vertices, a reference to a vector of vertices and the convexity
0114   VECCORE_ATT_HOST_DEVICE
0115   Polygon(size_t n, vector_t<Vec_t> &vertices, bool convex = false);
0116
0117   // @brief Fast constructor with no checks in case of convex polygons, providing the normal vector (normalized)
0118   VECCORE_ATT_HOST_DEVICE
0119   Polygon(size_t n, vector_t<Vec_t> &vertices, Vec_t const &norm);
0120
0121   /// Contructor with vertices and indices. Calls CheckAndFixDegenerate() on the polygon.
0122   Polygon(size_t n, vector_t<Vec_t> &vertices, vector_t<size_t> const &indices, bool convex);
0123
0124   /// @brief Copy constructor
0125   VECCORE_ATT_HOST_DEVICE
0126   Polygon(const Polygon &other)
0127       : fN(other.fN), fConvex(other.fConvex), fHasNorm(other.fHasNorm), fValid(other.fValid), fDist(other.fDist),
0128         fNorm(other.fNorm), fVert(other.fVert), fInd(other.fInd), fSides(other.fSides)
0129   {
0130   }
0131
0132   /// @brief Assignment operator
0133   VECGEOM_FORCE_INLINE
0134   Polygon &operator=(const Polygon &other)
0135   {
0136     if (&other == this) return *this;
0137     fN       = other.fN;
0138     fConvex  = other.fConvex;
0139     fHasNorm = other.fHasNorm;
0140     fValid   = other.fValid;
0141     fDist    = other.fDist;
0142     fNorm    = other.fNorm;
0143     fVert    = other.fVert;
0144     fInd     = other.fInd;
0145     fSides   = other.fSides;
0146     return *this;
0147   }
0148
0149   /// @brief Setter for a vertex index
0150   VECGEOM_FORCE_INLINE
0151   void SetVertex(size_t ind, size_t ivert) { fInd[ind] = ivert; }
0152
0153   /// @brief Getter for a vertex
0154   VECGEOM_FORCE_INLINE
0155   Vec_t const &GetVertex(size_t i) const { return fVert[fInd[i]]; }
0156
0157   /// @brief Setter from an array of vertex indices
0158   template <typename T>
0159   void SetVertices(vector_t<T> indices)
0160   {
0161     for (size_t i = 0; i < fN; ++i)
0162       fInd[i] = size_t(indices[i]);
0163   }
0164
0165   /// @brief Initialization is mandatory before first use
0166   void Init();
0167
0168   /// @brief Transform the polygon by a general transformation
0169   void Transform(Transformation3D const &tr);
0170
0171   /// Makes best effort to fix a degenerate polygon. In case it cannot be fixed, marks polygon as invalid.
0172   /// Only repeated vertices, and not collinear vertices are handled.
0173   void CheckAndFixDegenerate();
0174
0175   /// Decomposes a polygon into triangles.
0176   /// @param[out] polys The resulting triangles.
0177   void TriangulatePolygon(vector_t<Polygon> &polys) const;
0178
0179   /// Checks if fVert[i1] is a convex vertex in
0180   /// the polygon given its next neighbor fVert[i2] and previous neighbor fVert[i0].
0181   bool isConvexVertex(size_t i0, size_t i1, size_t i2) const;
0182
0183   /// Checks if a point is inside the triangle formed by fVert[i0]-fVert[i1]-fVert[i2]
0184   /// Point should already be on the plane.
0185   bool isPointInsideTriangle(const Vec_t &point, size_t i0, size_t i1, size_t i2) const;
0186
0187   /// Checks if a point is inside the polygon
0188   /// Point should already be on the plane.
0189   bool IsPointInside(const Vec_t &p) const;
0190
0191 #ifndef VECCORE_CUDA
0192   /// Computes the overlapping part of two polygons
0193   struct PolygonIntersection *Intersect(const Polygon &clipper);
0194 #endif
0195
0196   /// Retrieves the extent of the polygon
0197   /// @param[out] x Polygon boundaries in x axis, x[0] holding the smallest x value.
0198   /// @param[out] y Polygon boundaries in y axis, y[0] holding the smallest y value.
0199   /// @param[out] z Polygon boundaries in z axis, z[0] holding the smallest z value.
0200   void Extent(Precision x[2], Precision y[2], Precision z[2]);
0201 };
0202
0203 /// @brief A simple polyhedron defined by vertices and polygons
0204 struct Polyhedron {
0205   vector_t<Vec_t> fVert;    ///< Vector of vertices
0206   vector_t<Polygon> fPolys; ///< Vector of polygons
0207
0208   /// @brief Constructors
0209   Polyhedron(){};
0210   Polyhedron(size_t nvert, size_t npolys)
0211   {
0212     fVert.reserve(nvert);
0213     fPolys.reserve(npolys);
0214   }
0215
0216   /// @brief Polyhedrons can be re-used
0217   VECGEOM_FORCE_INLINE
0218   void Reset(size_t nvert, size_t npolys)
0219   {
0220     fVert.reserve(nvert);
0221     fVert.clear();
0222     fPolys.reserve(npolys);
0223     fPolys.clear();
0224   }
0225
0226   /// @brief Get number of vertices
0227   VECGEOM_FORCE_INLINE
0228   size_t GetNvertices() const { return fVert.size(); }
0229
0230   /// @brief Get number of polygons
0231   VECGEOM_FORCE_INLINE
0232   size_t GetNpolygons() const { return fPolys.size(); }
0233
0234   /// @brief Get a reference to a specific vertex
0235   VECGEOM_FORCE_INLINE
0236   Vec_t const &GetVertex(size_t i) const { return fVert[i]; }
0237
0238   /// @brief Get a reference to a specific polygon
0239   VECGEOM_FORCE_INLINE
0240   Polygon const &GetPolygon(size_t i) const { return fPolys[i]; }
0241
0242   /// @brief Transform the polygon by a general transformation
0243   void Transform(Transformation3D const &tr);
0244
0245   /// Adds given polygon to this polyhedron
0246   void AddPolygon(Polygon &poly, bool triangulate);
0247 };
0248
0249 /// Structure to hold the result of Polygon-Polygon intersection.
0250 ///
0251 ///
0252 ///              +------+         +------+   +-------+
0253 ///              |      |         |      |   |       |
0254 ///              |      |         |      |   |    xxxx----+
0255 ///              +------x------+  xxxxxxxx   |    x  x    |
0256 ///                     |      |  |      |   |    x  x    |
0257 ///                     |      |  |      |   |    xxxx----+
0258 ///                     +------+  +------+   +-------+
0259 ///                 point           line          polygon
0260
0261 struct PolygonIntersection {
0262   vector_t<Vec_t> fPoints;     ///< Resulting points from the intersection
0263   vector_t<Line> fLines;       ///< Resulting lines from the intersection
0264   vector_t<Polygon> fPolygons; ///< Resulting polygons from the intersection
0265   vector_t<Vec_t> fVertices;   ///< Resulting vertices in case there are polygons
0266 };
0267
0268 #ifndef VECCORE_CUDA
0269 /// @brief Streamers
0270 std::ostream &operator<<(std::ostream &os, Plane const &hpl);
0271 std::ostream &operator<<(std::ostream &os, Polygon const &poly);
0272 std::ostream &operator<<(std::ostream &os, Polyhedron const &polyh);
0273 #endif
0274
0275 /// @brief Function to find the crossing line between two planes.
0276 EPlaneXing_t PlaneXing(Plane const &pl1, Plane const &pl2, Vector3D<Precision> &point, Vector3D<Precision> &direction);
0277
0278 // @brief Function to find if 2 arbitrary polygons cross each other.
0279 EBodyXing_t PolygonXing(Polygon const &poly1, Polygon const &poly2, Line *line = nullptr);
0280
0281 // @brief Function to find if 2 arbitrary polyhedrons cross each other.
0282 EBodyXing_t PolyhedronXing(Polyhedron const &poly1, Polyhedron const &poly2, vector_t<Line> &lines);
0283
0284 /// @brief Function to determine crossing of two arbitrary placed boxes
0285 /** The function takes the box parameters and their transformations in a common frame.
0286     A fast check is performed if both transformations are identity. */
0287 EBodyXing_t BoxCollision(Vector3D<Precision> const &box1, Transformation3D const &tr1, Vector3D<Precision> const &box2,
0288                          Transformation3D const &tr2);
0289
0290 /// @brief Function filling a polyhedron as a box
0291 void FillBoxPolyhedron(Vec_t const &dimensions, Polyhedron &polyh);
0292
0293 } // namespace Utils3D
0294 } // namespace VECGEOM_IMPL_NAMESPACE
0295 } // namespace vecgeom
0296 #endif