Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:11:03

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