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{nullptr}; ///< 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   Polygon() = default;
0114   /// @brief Constructor taking the number of vertices, a reference to a vector of vertices and the convexity
0115   VECCORE_ATT_HOST_DEVICE
0116   Polygon(size_t n, vector_t<Vec_t> &vertices, bool convex = false);
0117
0118   // @brief Fast constructor with no checks in case of convex polygons, providing the normal vector (normalized)
0119   VECCORE_ATT_HOST_DEVICE
0120   Polygon(size_t n, vector_t<Vec_t> &vertices, Vec_t const &norm);
0121
0122   /// Contructor with vertices and indices. Calls CheckAndFixDegenerate() on the polygon.
0123   Polygon(size_t n, vector_t<Vec_t> &vertices, vector_t<size_t> const &indices, bool convex);
0124
0125   /// @brief Copy constructor
0126   Polygon(const Polygon &other) = default;
0127
0128   /// @brief Assignment operator
0129
0130   Polygon &operator=(const Polygon &other) = default;
0131
0132   /// @brief Setter for a vertex index
0133   VECGEOM_FORCE_INLINE
0134   void SetVertex(size_t ind, size_t ivert) { fInd[ind] = ivert; }
0135
0136   /// @brief Getter for a vertex
0137   VECGEOM_FORCE_INLINE
0138   Vec_t const &GetVertex(size_t i) const { return (*fVert)[fInd[i]]; }
0139
0140   /// @brief Setter from an array of vertex indices
0141   template <typename T>
0142   void SetVertices(vector_t<T> indices)
0143   {
0144     for (size_t i = 0; i < fN; ++i)
0145       fInd[i] = size_t(indices[i]);
0146   }
0147
0148   /// @brief Initialization is mandatory before first use
0149   void Init();
0150
0151   /// @brief Transform the polygon by a general transformation
0152   void Transform(Transformation3D const &tr);
0153
0154   /// Makes best effort to fix a degenerate polygon. In case it cannot be fixed, marks polygon as invalid.
0155   /// Only repeated vertices, and not collinear vertices are handled.
0156   void CheckAndFixDegenerate();
0157
0158   /// Decomposes a polygon into triangles.
0159   /// @param[out] polys The resulting triangles.
0160   void TriangulatePolygon(vector_t<Polygon> &polys) const;
0161
0162   /// Checks if fVert[i1] is a convex vertex in
0163   /// the polygon given its next neighbor fVert[i2] and previous neighbor fVert[i0].
0164   bool isConvexVertex(size_t i0, size_t i1, size_t i2) const;
0165
0166   /// Checks if a point is inside the triangle formed by fVert[i0]-fVert[i1]-fVert[i2]
0167   /// Point should already be on the plane.
0168   bool isPointInsideTriangle(const Vec_t &point, size_t i0, size_t i1, size_t i2) const;
0169
0170   /// Checks if a point is inside the polygon
0171   /// Point should already be on the plane.
0172   bool IsPointInside(const Vec_t &p) const;
0173
0174 #ifndef VECCORE_CUDA
0175   /// Computes the overlapping part of two polygons
0176   struct PolygonIntersection *Intersect(const Polygon &clipper);
0177 #endif
0178
0179   /// Retrieves the extent of the polygon
0180   /// @param[out] x Polygon boundaries in x axis, x[0] holding the smallest x value.
0181   /// @param[out] y Polygon boundaries in y axis, y[0] holding the smallest y value.
0182   /// @param[out] z Polygon boundaries in z axis, z[0] holding the smallest z value.
0183   void Extent(Precision x[2], Precision y[2], Precision z[2]);
0184 };
0185
0186 /// @brief A simple polyhedron defined by vertices and polygons
0187 struct Polyhedron {
0188   vector_t<Vec_t> fVert;    ///< Vector of vertices
0189   vector_t<Polygon> fPolys; ///< Vector of polygons
0190
0191   /// @brief Constructors
0192   Polyhedron() = default;
0193   Polyhedron(size_t nvert, size_t npolys)
0194   {
0195     fVert.reserve(nvert);
0196     fPolys.reserve(npolys);
0197   }
0198
0199   /// @brief Polyhedrons can be re-used
0200   VECGEOM_FORCE_INLINE
0201   void Reset(size_t nvert, size_t npolys)
0202   {
0203     fVert.reserve(nvert);
0204     fVert.clear();
0205     fPolys.reserve(npolys);
0206     fPolys.clear();
0207   }
0208
0209   /// @brief Get number of vertices
0210   VECGEOM_FORCE_INLINE
0211   size_t GetNvertices() const { return fVert.size(); }
0212
0213   /// @brief Get number of polygons
0214   VECGEOM_FORCE_INLINE
0215   size_t GetNpolygons() const { return fPolys.size(); }
0216
0217   /// @brief Get a reference to a specific vertex
0218   VECGEOM_FORCE_INLINE
0219   Vec_t const &GetVertex(size_t i) const { return fVert[i]; }
0220
0221   /// @brief Get a reference to a specific polygon
0222   VECGEOM_FORCE_INLINE
0223   Polygon const &GetPolygon(size_t i) const { return fPolys[i]; }
0224
0225   /// @brief Transform the polygon by a general transformation
0226   void Transform(Transformation3D const &tr);
0227
0228   /// Adds given polygon to this polyhedron
0229   void AddPolygon(Polygon &poly, bool triangulate);
0230 };
0231
0232 /// Structure to hold the result of Polygon-Polygon intersection.
0233 ///
0234 ///
0235 ///              +------+         +------+   +-------+
0236 ///              |      |         |      |   |       |
0237 ///              |      |         |      |   |    xxxx----+
0238 ///              +------x------+  xxxxxxxx   |    x  x    |
0239 ///                     |      |  |      |   |    x  x    |
0240 ///                     |      |  |      |   |    xxxx----+
0241 ///                     +------+  +------+   +-------+
0242 ///                 point           line          polygon
0243
0244 struct PolygonIntersection {
0245   vector_t<Vec_t> fPoints;     ///< Resulting points from the intersection
0246   vector_t<Line> fLines;       ///< Resulting lines from the intersection
0247   vector_t<Polygon> fPolygons; ///< Resulting polygons from the intersection
0248   vector_t<Vec_t> fVertices;   ///< Resulting vertices in case there are polygons
0249 };
0250
0251 #ifndef VECCORE_CUDA
0252 /// @brief Streamers
0253 std::ostream &operator<<(std::ostream &os, Plane const &hpl);
0254 std::ostream &operator<<(std::ostream &os, Polygon const &poly);
0255 std::ostream &operator<<(std::ostream &os, Polyhedron const &polyh);
0256 #endif
0257
0258 /// @brief Function to find the crossing line between two planes.
0259 EPlaneXing_t PlaneXing(Plane const &pl1, Plane const &pl2, Vector3D<Precision> &point, Vector3D<Precision> &direction);
0260
0261 // @brief Function to find if 2 arbitrary polygons cross each other.
0262 EBodyXing_t PolygonXing(Polygon const &poly1, Polygon const &poly2, Line *line = nullptr);
0263
0264 // @brief Function to find if 2 arbitrary polyhedrons cross each other.
0265 EBodyXing_t PolyhedronXing(Polyhedron const &poly1, Polyhedron const &poly2, vector_t<Line> &lines);
0266
0267 /// @brief Function to determine crossing of two arbitrary placed boxes
0268 /** The function takes the box parameters and their transformations in a common frame.
0269     A fast check is performed if both transformations are identity. */
0270 EBodyXing_t BoxCollision(Vector3D<Precision> const &box1, Transformation3D const &tr1, Vector3D<Precision> const &box2,
0271                          Transformation3D const &tr2);
0272
0273 /// @brief Function filling a polyhedron as a box
0274 void FillBoxPolyhedron(Vec_t const &dimensions, Polyhedron &polyh);
0275
0276 } // namespace Utils3D
0277 } // namespace VECGEOM_IMPL_NAMESPACE
0278 } // namespace vecgeom
0279 #endif