File indexing completed on 2025-01-19 09:23:34
0001
0002
0003
0004
0005
0006
0007
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
0020 namespace Acts::detail::VerticesHelper {
0021
0022
0023
0024
0025
0026
0027
0028
0029
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
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
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
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
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
0073
0074
0075
0076
0077
0078
0079
0080
0081
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
0089
0090
0091
0092
0093
0094
0095
0096 std::vector<Vector2> circularVertices(ActsScalar innerR, ActsScalar outerR,
0097 ActsScalar avgPhi = 0.,
0098 ActsScalar halfPhi = M_PI,
0099 unsigned int lseg = 1);
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110 template <typename vertex_t, typename vertex_container_t>
0111 bool isInsidePolygon(const vertex_t& point,
0112 const vertex_container_t& vertices) {
0113
0114
0115
0116
0117
0118
0119
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
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
0139 if (lineSide(l1, *std::begin(vertices)) != reference) {
0140 return false;
0141 }
0142
0143 return true;
0144 }
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
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
0164
0165
0166
0167
0168 bool onHyperPlane(const std::vector<Vector3>& vertices,
0169 ActsScalar tolerance = s_onSurfaceTolerance);
0170
0171
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
0181
0182 auto closestOnSegment = [&](auto&& ll0, auto&& ll1) {
0183
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;
0190
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
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
0211 Vector2 last = closestOnSegment(l1, *std::begin(vertices));
0212 if (squaredNorm(last - point) < closestDist) {
0213 closest = last;
0214 }
0215 return closest;
0216 }
0217
0218
0219 inline Vector2 computeEuclideanClosestPointOnRectangle(
0220 const Vector2& point, const Vector2& lowerLeft, const Vector2& upperRight) {
0221
0222
0223
0224
0225
0226
0227
0228
0229
0230
0231
0232
0233
0234
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
0244 if (loc0Min <= l0 && l0 < loc0Max && loc1Min <= l1 && l1 < loc1Max) {
0245
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
0268 if (l0 > loc0Max) {
0269 if (l1 > loc1Max) {
0270 return {loc0Max, loc1Max};
0271 } else if (l1 <= loc1Min) {
0272 return {loc0Max, loc1Min};
0273 } else {
0274 return {loc0Max, l1};
0275 }
0276 } else if (l0 < loc0Min) {
0277 if (l1 > loc1Max) {
0278 return {loc0Min, loc1Max};
0279 } else if (l1 <= loc1Min) {
0280 return {loc0Min, loc1Min};
0281 } else {
0282 return {loc0Min, l1};
0283 }
0284 } else {
0285 if (l1 > loc1Max) {
0286 return {l0, loc1Max};
0287 } else {
0288 return {l0, loc1Min};
0289 }
0290
0291 }
0292 }
0293 }
0294
0295 }