File indexing completed on 2025-01-18 09:11:03
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 <numbers>
0017 #include <span>
0018 #include <utility>
0019 #include <vector>
0020
0021
0022 namespace Acts::detail::VerticesHelper {
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
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
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
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
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
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
0093
0094
0095
0096
0097
0098
0099
0100
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
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116 template <typename vertex_t, typename vertex_container_t>
0117 bool isInsidePolygon(const vertex_t& point,
0118 const vertex_container_t& vertices) {
0119
0120
0121
0122
0123
0124
0125
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
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
0145 if (lineSide(l1, *std::begin(vertices)) != reference) {
0146 return false;
0147 }
0148
0149 return true;
0150 }
0151
0152
0153
0154
0155
0156
0157
0158
0159
0160
0161
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
0170
0171
0172
0173
0174 bool onHyperPlane(const std::vector<Vector3>& vertices,
0175 double tolerance = s_onSurfaceTolerance);
0176
0177
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
0186
0187 auto closestOnSegment = [&](auto&& ll0, auto&& ll1) {
0188
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;
0195
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
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
0216 Vector2 last = closestOnSegment(l1, *std::begin(vertices));
0217 if (squaredNorm(last - point) < closestDist) {
0218 closest = last;
0219 }
0220 return closest;
0221 }
0222
0223
0224 inline Vector2 computeEuclideanClosestPointOnRectangle(
0225 const Vector2& point, const Vector2& lowerLeft, const Vector2& upperRight) {
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239
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
0249 if (loc0Min <= l0 && l0 < loc0Max && loc1Min <= l1 && l1 < loc1Max) {
0250
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
0273 if (l0 > loc0Max) {
0274 if (l1 > loc1Max) {
0275 return {loc0Max, loc1Max};
0276 } else if (l1 <= loc1Min) {
0277 return {loc0Max, loc1Min};
0278 } else {
0279 return {loc0Max, l1};
0280 }
0281 } else if (l0 < loc0Min) {
0282 if (l1 > loc1Max) {
0283 return {loc0Min, loc1Max};
0284 } else if (l1 <= loc1Min) {
0285 return {loc0Min, loc1Min};
0286 } else {
0287 return {loc0Min, l1};
0288 }
0289 } else {
0290 if (l1 > loc1Max) {
0291 return {l0, loc1Max};
0292 } else {
0293 return {l0, loc1Min};
0294 }
0295
0296 }
0297 }
0298 }
0299
0300 }