src/Surfaces/IntersectionHelper2D.cpp

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 #include "Acts/Surfaces/detail/IntersectionHelper2D.hpp"
0010
0011 #include "Acts/Definitions/Tolerance.hpp"
0012 #include "Acts/Utilities/VectorHelpers.hpp"
0013 #include "Acts/Utilities/detail/RealQuadraticEquation.hpp"
0014
0015 #include <cmath>
0016
0017 Acts::Intersection2D Acts::detail::IntersectionHelper2D::intersectSegment(
0018     const Vector2& s0, const Vector2& s1, const Vector2& origin,
0019     const Vector2& dir, bool boundCheck) {
0020   using Line = Eigen::ParametrizedLine<ActsScalar, 2>;
0021   using Plane = Eigen::Hyperplane<ActsScalar, 2>;
0022
0023   Vector2 edge(s1 - s0);
0024   ActsScalar det = edge.x() * dir.y() - edge.y() * dir.x();
0025   if (std::abs(det) < s_epsilon) {
0026     return Intersection2D::invalid();
0027   }
0028
0029   auto line = Line(origin, dir);
0030   auto d = line.intersectionParameter(Plane::Through(s0, s1));
0031
0032   Vector2 intersection(origin + d * dir);
0033   Intersection2D::Status status = Intersection2D::Status::reachable;
0034   if (boundCheck) {
0035     auto edgeToSol = intersection - s0;
0036     if (edgeToSol.dot(edge) < 0. || edgeToSol.norm() > (edge).norm()) {
0037       status = Intersection2D::Status::unreachable;
0038     }
0039   }
0040   return Intersection2D(intersection, d, status);
0041 }
0042
0043 std::array<Acts::Intersection2D, 2>
0044 Acts::detail::IntersectionHelper2D::intersectEllipse(ActsScalar Rx,
0045                                                      ActsScalar Ry,
0046                                                      const Vector2& origin,
0047                                                      const Vector2& dir) {
0048   auto createSolution =
0049       [&](const Vector2& sol,
0050           const Vector2& alt) -> std::array<Acts::Intersection2D, 2> {
0051     Vector2 toSolD(sol - origin);
0052     Vector2 toAltD(alt - origin);
0053
0054     ActsScalar solD = std::copysign(toSolD.norm(), toSolD.dot(dir));
0055     ActsScalar altD = std::copysign(toAltD.norm(), toAltD.dot(dir));
0056
0057     if (std::abs(solD) < std::abs(altD)) {
0058       return {Intersection2D(sol, solD, Intersection2D::Status::reachable),
0059               Intersection2D(alt, altD, Intersection2D::Status::reachable)};
0060     }
0061     return {Intersection2D(alt, altD, Intersection2D::Status::reachable),
0062             Intersection2D(sol, solD, Intersection2D::Status::reachable)};
0063   };
0064
0065   // Special cases first
0066   if (std::abs(dir.x()) < s_epsilon) {
0067     ActsScalar solx = origin.x();
0068     ActsScalar D = 1. - solx * solx / (Rx * Rx);
0069     if (D > 0.) {
0070       ActsScalar sqrtD = std::sqrt(D);
0071       Vector2 sol(solx, Ry * sqrtD);
0072       Vector2 alt(solx, -Ry * sqrtD);
0073       return createSolution(sol, alt);
0074     } else if (std::abs(D) < s_epsilon) {
0075       return {Intersection2D(Vector2(solx, 0.), -origin.y(),
0076                              Intersection2D::Status::reachable),
0077               Intersection2D::invalid()};
0078     }
0079     return {Intersection2D::invalid(), Intersection2D::invalid()};
0080   } else if (std::abs(dir.y()) < s_epsilon) {
0081     ActsScalar soly = origin.y();
0082     ActsScalar D = 1. - soly * soly / (Ry * Ry);
0083     if (D > 0.) {
0084       ActsScalar sqrtD = std::sqrt(D);
0085       Vector2 sol(Rx * sqrtD, soly);
0086       Vector2 alt(-Rx * sqrtD, soly);
0087       return createSolution(sol, alt);
0088     } else if (std::abs(D) < s_epsilon) {
0089       return {Intersection2D(Vector2(0., soly), -origin.x(),
0090                              Intersection2D::Status::reachable),
0091               Intersection2D::invalid()};
0092     }
0093     return {Intersection2D::invalid(), Intersection2D::invalid()};
0094   }
0095   // General solution
0096   ActsScalar k = dir.y() / dir.x();
0097   ActsScalar d = origin.y() - k * origin.x();
0098   ActsScalar Ry2 = Ry * Ry;
0099   ActsScalar alpha = 1. / (Rx * Rx) + k * k / Ry2;
0100   ActsScalar beta = 2. * k * d / Ry2;
0101   ActsScalar gamma = d * d / Ry2 - 1;
0102   Acts::detail::RealQuadraticEquation solver(alpha, beta, gamma);
0103   if (solver.solutions == 1) {
0104     ActsScalar x = solver.first;
0105     Vector2 sol(x, k * x + d);
0106     Vector2 toSolD(sol - origin);
0107     ActsScalar solD = std::copysign(toSolD.norm(), toSolD.dot(dir));
0108     return {Intersection2D(sol, solD, Intersection2D::Status::reachable),
0109             Intersection2D::invalid()};
0110   } else if (solver.solutions > 1) {
0111     ActsScalar x0 = solver.first;
0112     ActsScalar x1 = solver.second;
0113     Vector2 sol(x0, k * x0 + d);
0114     Vector2 alt(x1, k * x1 + d);
0115     return createSolution(sol, alt);
0116   }
0117   return {Intersection2D::invalid(), Intersection2D::invalid()};
0118 }
0119
0120 Acts::Intersection2D Acts::detail::IntersectionHelper2D::intersectCircleSegment(
0121     ActsScalar R, ActsScalar phiMin, ActsScalar phiMax, const Vector2& origin,
0122     const Vector2& dir) {
0123   auto intersections = intersectCircle(R, origin, dir);
0124   for (const auto& candidate : intersections) {
0125     if (candidate.pathLength() > 0.) {
0126       ActsScalar phi = Acts::VectorHelpers::phi(candidate.position());
0127       if (phi > phiMin && phi < phiMax) {
0128         return candidate;
0129       }
0130     }
0131   }
0132   return Intersection2D::invalid();
0133 }