Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-19 09:23:34

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008 
0009 #pragma once
0010 
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Tolerance.hpp"
0013 
0014 #include <algorithm>
0015 #include <cmath>
0016 #include <utility>
0017 #include <vector>
0018 
0019 /// Helper methods for polyhedron vertices drawing and inside/outside checks.
0020 namespace Acts::detail::VerticesHelper {
0021 
0022 /// A method that inserts the cartesian extrema points and segments
0023 /// a curved segment into sub segments
0024 ///
0025 /// @param phiMin the minimum Phi of the bounds object
0026 /// @param phiMax the maximum Phi of the bounds object
0027 /// @param phiRef is a vector of reference phi values to be included as well
0028 /// @param phiTolerance is the tolerance for reference phi insertion
0029 /// @return a vector
0030 std::vector<ActsScalar> phiSegments(ActsScalar phiMin = -M_PI,
0031                                     ActsScalar phiMax = M_PI,
0032                                     const std::vector<ActsScalar>& phiRefs = {},
0033                                     ActsScalar phiTolerance = 1e-6);
0034 
0035 /// Helper method to create a regular 2 or 3 D segment
0036 ///  between two phi values
0037 ///
0038 /// @tparam vertex_t Type of vertex to be applied
0039 /// @tparam transform_t Optional transform
0040 ///
0041 /// @param vertices [in,out] The 3D vertices to be filled
0042 /// @param rxy The radius description if first +/= second: ellipse
0043 /// @param phi1 The first phi value
0044 /// @param phi2 The second phi value
0045 /// @param lseg The number of segments for full 2*PI
0046 /// @param addon The additional segments to be built
0047 /// @param offset The out of plane offset position of the bow
0048 /// @param transform The transform applied (optional)
0049 template <typename vertex_t, typename transform_t>
0050 void createSegment(std::vector<vertex_t>& vertices,
0051                    std::pair<ActsScalar, ActsScalar> rxy, ActsScalar phi1,
0052                    ActsScalar phi2, unsigned int lseg, int addon = 0,
0053                    const vertex_t& offset = vertex_t::Zero(),
0054                    const transform_t& transform = transform_t::Identity()) {
0055   // Calculate the number of segments - 1 is the minimum
0056   unsigned int segs =
0057       static_cast<unsigned int>(std::abs(phi2 - phi1) / (2 * M_PI) * lseg);
0058   segs = segs > 0 ? segs : 1;
0059   ActsScalar phistep = (phi2 - phi1) / segs;
0060   // Create the segments
0061   for (unsigned int iphi = 0; iphi < segs + addon; ++iphi) {
0062     ActsScalar phi = phi1 + iphi * phistep;
0063     vertex_t vertex = vertex_t::Zero();
0064     vertex(0) = rxy.first * std::cos(phi);
0065     vertex(1) = rxy.second * std::sin(phi);
0066 
0067     vertex = vertex + offset;
0068     vertices.push_back(transform * vertex);
0069   }
0070 }
0071 
0072 /// Construct vertices on an ellipse-like bound object.
0073 ///
0074 /// @param innerRx The radius of the inner ellipse (in x), 0 if sector
0075 /// @param innerRy The radius of the inner ellipse (in y), 0 if sector
0076 /// @param outerRx The radius of the outer ellipse (in x)
0077 /// @param outerRy The radius of the outer ellipse (in y)
0078 /// @param avgPhi The phi direction of the center if sector
0079 /// @param halfPhi The half phi sector if sector
0080 /// @param lseg The number of segments for for a full 2*pi segment
0081 /// @return a vector of 2d-vectors
0082 std::vector<Vector2> ellipsoidVertices(ActsScalar innerRx, ActsScalar innerRy,
0083                                        ActsScalar outerRx, ActsScalar outerRy,
0084                                        ActsScalar avgPhi = 0.,
0085                                        ActsScalar halfPhi = M_PI,
0086                                        unsigned int lseg = 1);
0087 
0088 /// Construct vertices on an disc/wheel-like bound object.
0089 ///
0090 /// @param innerR The radius of the inner circle (sector)
0091 /// @param outerR The radius of the outer circle (sector)
0092 /// @param avgPhi The phi direction of the center if sector
0093 /// @param halfPhi The half phi sector if sector
0094 /// @param lseg The number of segments for for a full 2*pi segment
0095 /// @return a vector of 2d-vectors
0096 std::vector<Vector2> circularVertices(ActsScalar innerR, ActsScalar outerR,
0097                                       ActsScalar avgPhi = 0.,
0098                                       ActsScalar halfPhi = M_PI,
0099                                       unsigned int lseg = 1);
0100 /// Check if the point is inside the polygon w/o any tolerances.
0101 ///
0102 /// @tparam vertex_container_t is an iterable container
0103 ///
0104 /// @param point is the Vector2Type to check
0105 /// @param vertices Forward iterable container of convex polygon vertices.
0106 ///                 Calling `std::begin`/ `std::end` on the container must
0107 ///                 return an iterator where `*it` must be convertible to
0108 ///                 an `Vector2Type`.
0109 /// @return bool for inside/outside
0110 template <typename vertex_t, typename vertex_container_t>
0111 bool isInsidePolygon(const vertex_t& point,
0112                      const vertex_container_t& vertices) {
0113   // when we move along the edges of a convex polygon, a point on the inside of
0114   // the polygon will always appear on the same side of each edge.
0115   // a point on the outside will switch sides at least once.
0116 
0117   // returns which side of the connecting line between `ll0` and `ll1` the point
0118   // `p` is on. computes the sign of the z-component of the cross-product
0119   // between the line normal vector and the vector from `ll0` to `p`.
0120   auto lineSide = [&](auto&& ll0, auto&& ll1) {
0121     auto normal = ll1 - ll0;
0122     auto delta = point - ll0;
0123     return std::signbit((normal[0] * delta[1]) - (normal[1] * delta[0]));
0124   };
0125 
0126   auto iv = std::begin(vertices);
0127   auto l0 = *iv;
0128   auto l1 = *(++iv);
0129   // use vertex0 to vertex1 to define reference sign and compare w/ all edges
0130   auto reference = lineSide(l0, l1);
0131   for (++iv; iv != std::end(vertices); ++iv) {
0132     l0 = l1;
0133     l1 = *iv;
0134     if (lineSide(l0, l1) != reference) {
0135       return false;
0136     }
0137   }
0138   // manual check for last edge from last vertex back to the first vertex
0139   if (lineSide(l1, *std::begin(vertices)) != reference) {
0140     return false;
0141   }
0142   // point was always on the same side. point must be inside.
0143   return true;
0144 }
0145 
0146 /// Check if the point is inside the rectangle.
0147 ///
0148 /// @tparam vertex_t is vector with [0],[1] access
0149 ///
0150 /// @param point is the Vector2Type to check
0151 /// @param vertices Forward iterable container of convex polygon vertices.
0152 ///                 Calling `std::begin`/ `std::end` on the container must
0153 ///                 return an iterator where `*it` must be convertible to
0154 ///                 an `Vector2Type`.
0155 /// @return bool for inside/outside
0156 template <typename vertex_t>
0157 bool isInsideRectangle(const vertex_t& point, const vertex_t& lowerLeft,
0158                        const vertex_t& upperRight) {
0159   return (lowerLeft[0] <= point[0]) && (point[0] < upperRight[0]) &&
0160          (lowerLeft[1] <= point[1]) && (point[1] < upperRight[1]);
0161 }
0162 
0163 /// This method checks if a cloud of points are on 2D hyper-plane in 3D space.
0164 ///
0165 /// @param vertices The list of vertices to test
0166 /// @param tolerance The allowed out of plane tolerance
0167 /// @return boolean to indicate if all points are inside/outside
0168 bool onHyperPlane(const std::vector<Vector3>& vertices,
0169                   ActsScalar tolerance = s_onSurfaceTolerance);
0170 
0171 /// Calculate the closest point on the polygon.
0172 template <typename Vector2Container>
0173 inline Vector2 computeClosestPointOnPolygon(const Vector2& point,
0174                                             const Vector2Container& vertices,
0175                                             const SquareMatrix2& metric) {
0176   auto squaredNorm = [&](const Vector2& x) {
0177     return (x.transpose() * metric * x).value();
0178   };
0179 
0180   // calculate the closest position on the segment between `ll0` and `ll1` to
0181   // the point as measured by the metric induced by the metric matrix
0182   auto closestOnSegment = [&](auto&& ll0, auto&& ll1) {
0183     // normal vector and position of the closest point along the normal
0184     auto n = ll1 - ll0;
0185     auto n_transformed = metric * n;
0186     auto f = n.dot(n_transformed);
0187     auto u = std::isnormal(f)
0188                  ? (point - ll0).dot(n_transformed) / f
0189                  : 0.5;  // ll0 and ll1 are so close it doesn't matter
0190     // u must be in [0, 1] to still be on the polygon segment
0191     return ll0 + std::clamp(u, 0.0, 1.0) * n;
0192   };
0193 
0194   auto iv = std::begin(vertices);
0195   Vector2 l0 = *iv;
0196   Vector2 l1 = *(++iv);
0197   Vector2 closest = closestOnSegment(l0, l1);
0198   auto closestDist = squaredNorm(closest - point);
0199   // Calculate the closest point on other connecting lines and compare distances
0200   for (++iv; iv != std::end(vertices); ++iv) {
0201     l0 = l1;
0202     l1 = *iv;
0203     Vector2 current = closestOnSegment(l0, l1);
0204     auto currentDist = squaredNorm(current - point);
0205     if (currentDist < closestDist) {
0206       closest = current;
0207       closestDist = currentDist;
0208     }
0209   }
0210   // final edge from last vertex back to the first vertex
0211   Vector2 last = closestOnSegment(l1, *std::begin(vertices));
0212   if (squaredNorm(last - point) < closestDist) {
0213     closest = last;
0214   }
0215   return closest;
0216 }
0217 
0218 /// Calculate the closest point on the box
0219 inline Vector2 computeEuclideanClosestPointOnRectangle(
0220     const Vector2& point, const Vector2& lowerLeft, const Vector2& upperRight) {
0221   /*
0222    *
0223    *        |                 |
0224    *   IV   |       V         | I
0225    *        |                 |
0226    *  ------------------------------
0227    *        |                 |
0228    *        |                 |
0229    *   VIII |     INSIDE      | VI
0230    *        |                 |
0231    *        |                 |
0232    *  ------------------------------
0233    *        |                 |
0234    *   III  |      VII        | II
0235    *        |                 |
0236    *
0237    */
0238 
0239   double l0 = point[0], l1 = point[1];
0240   double loc0Min = lowerLeft[0], loc0Max = upperRight[0];
0241   double loc1Min = lowerLeft[1], loc1Max = upperRight[1];
0242 
0243   // check if inside
0244   if (loc0Min <= l0 && l0 < loc0Max && loc1Min <= l1 && l1 < loc1Max) {
0245     // INSIDE
0246     double dist = std::abs(loc0Max - l0);
0247     Vector2 cls(loc0Max, l1);
0248 
0249     double test = std::abs(loc0Min - l0);
0250     if (test <= dist) {
0251       dist = test;
0252       cls = {loc0Min, l1};
0253     }
0254 
0255     test = std::abs(loc1Max - l1);
0256     if (test <= dist) {
0257       dist = test;
0258       cls = {l0, loc1Max};
0259     }
0260 
0261     test = std::abs(loc1Min - l1);
0262     if (test <= dist) {
0263       return {l0, loc1Min};
0264     }
0265     return cls;
0266   } else {
0267     // OUTSIDE, check sectors
0268     if (l0 > loc0Max) {
0269       if (l1 > loc1Max) {  // I
0270         return {loc0Max, loc1Max};
0271       } else if (l1 <= loc1Min) {  // II
0272         return {loc0Max, loc1Min};
0273       } else {  // VI
0274         return {loc0Max, l1};
0275       }
0276     } else if (l0 < loc0Min) {
0277       if (l1 > loc1Max) {  // IV
0278         return {loc0Min, loc1Max};
0279       } else if (l1 <= loc1Min) {  // III
0280         return {loc0Min, loc1Min};
0281       } else {  // VIII
0282         return {loc0Min, l1};
0283       }
0284     } else {
0285       if (l1 > loc1Max) {  // V
0286         return {l0, loc1Max};
0287       } else {  // l1 <= loc1Min # VII
0288         return {l0, loc1Min};
0289       }
0290       // third case not necessary, see INSIDE above
0291     }
0292   }
0293 }
0294 
0295 }  // namespace Acts::detail::VerticesHelper