Acts/Utilities/Intersection.hpp

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/Utilities/Logger.hpp"
0013
0014 #include <algorithm>
0015 #include <array>
0016 #include <cstddef>
0017 #include <cstdint>
0018 #include <limits>
0019 #include <span>
0020 #include <type_traits>
0021
0022 namespace Acts {
0023
0024 /// Status enum
0025 enum class IntersectionStatus : int {
0026   unreachable = 0,
0027   reachable = 1,
0028   onSurface = 2
0029 };
0030
0031 /// Ostream-operator for the IntersectionStatus enum
0032 /// @param os Output stream
0033 /// @param status IntersectionStatus to output
0034 /// @return Reference to output stream
0035 inline std::ostream& operator<<(std::ostream& os, IntersectionStatus status) {
0036   constexpr static std::array<const char*, 3> names = {
0037       {"missed/unreachable", "reachable", "onSurface"}};
0038
0039   os << names[static_cast<std::size_t>(status)];
0040   return os;
0041 }
0042
0043 /// Intersection struct containing the position, path length and status of an
0044 /// intersection.
0045 template <unsigned int DIM>
0046 class Intersection {
0047  public:
0048   /// Position type
0049   using Position = Eigen::Map<const Vector<DIM>>;
0050
0051   /// Constructor with arguments
0052   ///
0053   /// @param position is the position of the intersection
0054   /// @param pathLength is the path length to the intersection
0055   /// @param status is an enum indicating the status of the intersection
0056   constexpr Intersection(const Vector<DIM>& position, double pathLength,
0057                          IntersectionStatus status) noexcept
0058       : Intersection(std::span<const double, DIM>{position.data(), DIM},
0059                      pathLength, status) {}
0060
0061   /// Constructor from position vector, path length, and status
0062   /// @param position The intersection position
0063   /// @param pathLength The path length to the intersection
0064   /// @param status The intersection status
0065   constexpr Intersection(const Position& position, double pathLength,
0066                          IntersectionStatus status) noexcept
0067       : Intersection(std::span<const double, DIM>{position.data(), DIM},
0068                      pathLength, status) {}
0069
0070   /// Constructor from position span, path length, and status
0071   /// @param position Span of position coordinates
0072   /// @param pathLength The path length to the intersection
0073   /// @param status The intersection status
0074   constexpr Intersection(std::span<const double, DIM> position,
0075                          double pathLength, IntersectionStatus status) noexcept
0076       : m_pathLength(pathLength), m_status(status) {
0077     std::ranges::copy(position, m_position.begin());
0078   }
0079
0080   /// Copy constructor
0081   constexpr Intersection(const Intersection&) noexcept = default;
0082   /// Move constructor
0083   constexpr Intersection(Intersection&&) noexcept = default;
0084   /// Copy assignment operator
0085   /// @return Reference to this intersection for chaining
0086   constexpr Intersection& operator=(const Intersection&) noexcept = default;
0087   /// Move assignment operator
0088   /// @return Reference to this intersection for chaining
0089   constexpr Intersection& operator=(Intersection&&) noexcept = default;
0090
0091   /// Returns whether the intersection was successful or not
0092   /// @return True if intersection is reachable or on surface, false if unreachable
0093   constexpr bool isValid() const noexcept {
0094     return m_status != IntersectionStatus::unreachable;
0095   }
0096
0097   /// Returns the position of the interseciton
0098   /// @return Position vector of the intersection point
0099   Position position() const noexcept { return Position{m_position.data()}; }
0100
0101   /// Returns the path length to the intersection
0102   /// @return Signed path length from origin to intersection point
0103   constexpr double pathLength() const noexcept { return m_pathLength; }
0104
0105   /// Returns the intersection status enum
0106   /// @return Status indicating if intersection is unreachable, reachable, or on surface
0107   constexpr IntersectionStatus status() const noexcept { return m_status; }
0108
0109   /// Static factory to create an invalid intersection
0110   /// @return Invalid intersection with unreachable status
0111   constexpr static Intersection Invalid() noexcept { return Intersection(); }
0112
0113   /// Comparison function for path length order i.e. intersection closest to
0114   /// -inf will be first.
0115   /// @param aIntersection First intersection to compare
0116   /// @param bIntersection Second intersection to compare
0117   /// @return True if first intersection has smaller path length than second
0118   constexpr static bool pathLengthOrder(
0119       const Intersection& aIntersection,
0120       const Intersection& bIntersection) noexcept {
0121     auto a = aIntersection.pathLength();
0122     auto b = bIntersection.pathLength();
0123     return a < b;
0124   }
0125
0126   /// Comparison function for closest order i.e. intersection closest to 0 will
0127   /// be first.
0128   /// @param aIntersection First intersection to compare
0129   /// @param bIntersection Second intersection to compare
0130   /// @return True if first intersection is closer to zero path length than second
0131   constexpr static bool closestOrder(
0132       const Intersection& aIntersection,
0133       const Intersection& bIntersection) noexcept {
0134     using enum IntersectionStatus;
0135
0136     if ((aIntersection.status() == unreachable) &&
0137         (bIntersection.status() != unreachable)) {
0138       return false;
0139     }
0140     if ((aIntersection.status() != unreachable) &&
0141         (bIntersection.status() == unreachable)) {
0142       return true;
0143     }
0144     // both are reachable or onSurface now
0145     auto a = aIntersection.pathLength();
0146     auto b = bIntersection.pathLength();
0147     return std::abs(a) < std::abs(b);
0148   }
0149
0150   /// Comparison function for closest forward order i.e. intersection closest to
0151   /// 0 with positive path length will be first.
0152   /// @param aIntersection First intersection to compare
0153   /// @param bIntersection Second intersection to compare
0154   /// @return True if first intersection is closer to zero with preference for forward direction
0155   constexpr static bool closestForwardOrder(
0156       const Intersection& aIntersection,
0157       const Intersection& bIntersection) noexcept {
0158     auto a = aIntersection.pathLength();
0159     auto b = bIntersection.pathLength();
0160     return std::signbit(a) == std::signbit(b) ? std::abs(a) < std::abs(b)
0161                                               : a > b;
0162   }
0163
0164  private:
0165   /// Position of the intersection
0166   std::array<double, DIM> m_position{};
0167   /// Signed path length to the intersection (if valid)
0168   double m_pathLength = std::numeric_limits<double>::infinity();
0169   /// The Status of the intersection
0170   IntersectionStatus m_status = IntersectionStatus::unreachable;
0171
0172   constexpr Intersection() noexcept = default;
0173 };
0174
0175 /// Type alias for 2D intersection
0176 using Intersection2D = Intersection<2>;
0177 /// Type alias for 3D intersection
0178 using Intersection3D = Intersection<3>;
0179
0180 static_assert(std::is_trivially_copy_constructible_v<Intersection2D>);
0181 static_assert(std::is_trivially_move_constructible_v<Intersection2D>);
0182 static_assert(std::is_trivially_move_assignable_v<Intersection2D>);
0183
0184 /// Index type for intersections
0185 using IntersectionIndex = std::uint8_t;
0186 /// Maximum number of intersections that can be stored
0187 static constexpr IntersectionIndex s_maximumNumberOfIntersections = 2;
0188
0189 /// Container for up to two intersections in a given dimension.
0190 template <unsigned int DIM>
0191 class MultiIntersection {
0192  public:
0193   /// Intersection type for this dimension
0194   using IntersectionType = Intersection<DIM>;
0195   /// Pair of intersection and its index
0196   using IndexedIntersection = std::pair<IntersectionType, IntersectionIndex>;
0197
0198   /// Container type for storing intersections
0199   using Container =
0200       std::array<IntersectionType, s_maximumNumberOfIntersections>;
0201
0202   /// Size type for indexing
0203   using size_type = IntersectionIndex;
0204
0205   /// Construct from single intersection
0206   /// @param intersection The intersection
0207   constexpr explicit MultiIntersection(
0208       const IntersectionType& intersection) noexcept
0209       : m_intersections{intersection, IntersectionType::Invalid()}, m_size{1} {}
0210   /// Construct from two intersections
0211   /// @param intersection1 The first intersection
0212   /// @param intersection2 The second intersection
0213   constexpr MultiIntersection(const IntersectionType& intersection1,
0214                               const IntersectionType& intersection2) noexcept
0215       : m_intersections{intersection1, intersection2}, m_size{2} {}
0216
0217   /// Copy constructor
0218   constexpr MultiIntersection(const MultiIntersection&) noexcept = default;
0219   /// Move constructor
0220   constexpr MultiIntersection(MultiIntersection&&) noexcept = default;
0221   /// Copy assignment operator
0222   /// @return Reference to this object
0223   constexpr MultiIntersection& operator=(const MultiIntersection&) noexcept =
0224       default;
0225   /// Move assignment operator
0226   /// @return Reference to this object
0227   constexpr MultiIntersection& operator=(MultiIntersection&&) noexcept =
0228       default;
0229
0230   /// Access intersection by index
0231   /// @param index The index of the intersection
0232   /// @return Reference to the intersection
0233   constexpr const IntersectionType& operator[](IntersectionIndex index) const {
0234     return m_intersections[index];
0235   }
0236
0237   /// Access intersection at index with bounds checking
0238   /// @param index The index of the intersection
0239   /// @return Reference to the intersection
0240   constexpr const IntersectionType& at(IntersectionIndex index) const {
0241     return m_intersections.at(index);
0242   }
0243
0244   /// Get the number of intersections
0245   /// @return The number of intersections
0246   constexpr IntersectionIndex size() const noexcept { return m_size; }
0247
0248   /// Get begin iterator
0249   /// @return Iterator to the beginning
0250   constexpr auto begin() const noexcept {
0251     return std::span(m_intersections.data(), m_size).begin();
0252   }
0253   /// Get end iterator
0254   /// @return Iterator to the end
0255   constexpr auto end() const noexcept {
0256     return std::span(m_intersections.data(), m_size).end();
0257   }
0258
0259   /// Get closest intersection
0260   /// @return The closest intersection
0261   constexpr IntersectionType closest() const noexcept {
0262     return closestWithIndex().first;
0263   }
0264   /// Get closest intersection with its index
0265   /// @return Pair of intersection and its index
0266   constexpr IndexedIntersection closestWithIndex() const noexcept {
0267     auto min = std::ranges::min_element(m_intersections,
0268                                         IntersectionType::closestOrder);
0269     return {*min, static_cast<IntersectionIndex>(
0270                       std::distance(m_intersections.begin(), min))};
0271   }
0272
0273   /// Get closest forward intersection
0274   /// @return The closest forward intersection
0275   constexpr IntersectionType closestForward() const noexcept {
0276     return closestForwardWithIndex().first;
0277   }
0278   /// Get closest forward intersection with its index
0279   /// @return Pair of intersection and its index
0280   constexpr IndexedIntersection closestForwardWithIndex() const noexcept {
0281     auto min = std::ranges::min_element(m_intersections,
0282                                         IntersectionType::closestForwardOrder);
0283     return {*min, static_cast<IntersectionIndex>(
0284                       std::distance(m_intersections.begin(), min))};
0285   }
0286
0287  private:
0288   Container m_intersections{};
0289   IntersectionIndex m_size{};
0290 };
0291
0292 /// Container for up to two 2D intersections
0293 using MultiIntersection2D = MultiIntersection<2>;
0294 /// Container for up to two 3D intersections
0295 using MultiIntersection3D = MultiIntersection<3>;
0296
0297 static_assert(std::is_trivially_copy_constructible_v<MultiIntersection2D>);
0298 static_assert(std::is_trivially_move_constructible_v<MultiIntersection2D>);
0299 static_assert(std::is_trivially_move_assignable_v<MultiIntersection2D>);
0300
0301 namespace detail {
0302
0303 /// This function checks if an intersection path length is valid for the
0304 /// specified near-limit and far-limit
0305 ///
0306 /// @param pathLength The path length of the intersection
0307 /// @param nearLimit The minimum path length for an intersection to be considered
0308 /// @param farLimit The maximum path length for an intersection to be considered
0309 /// @param logger A optionally supplied logger which prints out a lot of infos
0310 ///               at VERBOSE level
0311 bool checkPathLength(double pathLength, double nearLimit, double farLimit,
0312                      const Logger& logger = getDummyLogger());
0313
0314 }  // namespace detail
0315
0316 }  // namespace Acts