feat: introduce navigation stream (#3538)

This PR introduces a `NavigationStream` object and a corresponding helper. 

The `NavigationStream` should build the backbone of the Gen3 navigator, with having a `TargetStream` and a `GeometryStream` to deal with.

For details of the ongoing discussion, please see #3526.


@paulgessinger @andiwand

Core/CMakeLists.txt

@@ -108,6 +108,7 @@ add_subdirectory(src/Detector)
 add_subdirectory(src/Geometry)
 add_subdirectory(src/MagneticField)
 add_subdirectory(src/Material)
+add_subdirectory(src/Navigation)
 add_subdirectory(src/Propagator)
 add_subdirectory(src/Surfaces)
 add_subdirectory(src/TrackFinding)

Core/include/Acts/Navigation/NavigationStream.hpp

@@ -0,0 +1,178 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Surfaces/BoundaryTolerance.hpp"
+#include "Acts/Utilities/Intersection.hpp"
+
+#include <tuple>
+#include <vector>
+
+namespace Acts {
+
+// To be removed when the namespace Experimental is omitted
+namespace Experimental {
+class Portal;
+}
+using namespace Experimental;
+
+class Surface;
+
+/// The NavigationStream is a container for the navigation candidates that
+/// are currentlu processed in a given context. The context could be local to a
+/// volume, or global to an entire track following.
+///
+/// The current candidates are stored in a vector of candidates, where an index
+/// is used to indicate the current active candidate.
+class NavigationStream {
+ public:
+  /// The query point for the navigation stream
+  ///
+  /// This holds the position and direction from which the navigation stream
+  /// should either be initialized or updated.
+  struct QueryPoint {
+    /// The position of the query point
+    Vector3 position = Vector3::Zero();
+    /// The direction of the query point
+    Vector3 direction = Vector3::Zero();
+  };
+
+  /// This is a candidate object of the navigation stream, it holds:
+  ///
+  /// - a Surface intersection
+  /// - a Portal : set if the surface represents a portal
+  /// - a BoundaryTolerance : the boundary tolerance used for the intersection
+  struct Candidate {
+    /// The intersection
+    ObjectIntersection<Surface> intersection =
+        ObjectIntersection<Surface>::invalid();
+    /// The portal
+    const Portal* portal = nullptr;
+    /// The boundary tolerance
+    BoundaryTolerance bTolerance = BoundaryTolerance::None();
+    /// Convenience access to surface
+    const Surface& surface() const { return *intersection.object(); }
+    /// Cinvencience access to the path length
+    ActsScalar pathLength() const { return intersection.pathLength(); }
+
+    /// Order along the path length
+    ///
+    /// @param aCandidate is the first candidate
+    /// @param bCandidate is the second candidate
+    ///
+    /// @return true if aCandidate is closer to the origin
+    constexpr static bool pathLengthOrder(const Candidate& aCandidate,
+                                          const Candidate& bCandidate) {
+      return ObjectIntersection<Surface>::pathLengthOrder(
+          aCandidate.intersection, bCandidate.intersection);
+    }
+  };
+
+  /// Switch to next next candidate
+  ///
+  /// @return true if a next candidate is available
+  bool switchToNextCandidate() {
+    if (m_currentIndex < m_candidates.size()) {
+      ++m_currentIndex;
+      return true;
+    }
+    return false;
+  }
+
+  /// Const access the current candidate
+  const Candidate& currentCandidate() const {
+    return m_candidates.at(m_currentIndex);
+  }
+
+  /// Current Index
+  std::size_t currentIndex() const { return m_currentIndex; }
+
+  /// Non-cost access the candidate vector
+  std::vector<Candidate>& candidates() { return m_candidates; }
+
+  /// Const access the candidate vector
+  const std::vector<Candidate>& candidates() const { return m_candidates; }
+
+  /// Non-cost access the current candidate
+  ///
+  /// This will throw and out of bounds exception if the stream is not
+  /// valid anymore.
+  Candidate& currentCandidate() { return m_candidates.at(m_currentIndex); }
+
+  /// The number of active candidates
+  std::size_t remainingCandidates() const {
+    return (m_candidates.size() - m_currentIndex);
+  }
+
+  /// Fill one surface into the candidate vector
+  ///
+  /// @param surface the surface to be filled
+  /// @param bTolerance the boundary tolerance used for the intersection
+  void addSurfaceCandidate(const Surface* surface,
+                           const BoundaryTolerance& bTolerance);
+
+  /// Fill n surfaces into the candidate vector
+  ///
+  /// @param surfaces the surfaces that are filled in
+  /// @param bTolerance the boundary tolerance used for the intersection
+  void addSurfaceCandidates(const std::vector<const Surface*>& surfaces,
+                            const BoundaryTolerance& bTolerance);
+
+  /// Fill one portal into the candidate vector
+  ///
+  /// @param portal the portals that are filled in
+  void addPortalCandidate(const Portal* portal);
+
+  /// Fill n portals into the candidate vector
+  ///
+  /// @param portals the portals that are filled in
+  void addPortalCandidates(const std::vector<const Portal*>& portals);
+
+  /// Initialize the stream from a query point
+  ///
+  /// @param gctx is the geometry context
+  /// @param queryPoint holds current position, direction, etc.
+  /// @param cTolerance is the candidate search tolerance
+  /// @param onSurfaceTolerance is the tolerance for on-surface intersections
+  ///
+  /// This method will first de-duplicate the candidates on basis of the surface
+  /// pointer to make sure that the multi-intersections are handled correctly.
+  /// This will allow intializeStream() to be called even as a re-initialization
+  /// and still work correctly with at one time valid candidates.
+  ///
+  /// @return true if the stream is active, false indicates that there are no valid candidates
+  bool initialize(const GeometryContext& gctx,
+                  const NavigationStream::QueryPoint& queryPoint,
+                  const BoundaryTolerance& cTolerance,
+                  ActsScalar onSurfaceTolerance = s_onSurfaceTolerance);
+
+  /// Convenience method to update a stream from a new query point,
+  /// this could be called from navigation delegates that do not require
+  /// a local state or from the navigator on the target stream
+  ///
+  /// @param gctx is the geometry context
+  /// @param queryPoint holds current position, direction, etc.
+  /// @param onSurfaceTolerance is the tolerance for on-surface intersections
+  ///
+  /// @return true if the stream is active, false indicate no valid candidates left
+  bool update(const GeometryContext& gctx,
+              const NavigationStream::QueryPoint& queryPoint,
+              ActsScalar onSurfaceTolerance = s_onSurfaceTolerance);
+
+ private:
+  /// The candidates of this navigation stream
+  std::vector<Candidate> m_candidates;
+
+  /// The currently active candidate
+  std::size_t m_currentIndex = 0u;
+};
+
+}  // namespace Acts

Core/include/Acts/Utilities/Intersection.hpp

@@ -69,12 +69,16 @@ class Intersection {
   /// Returns whether the intersection was successful or not
   constexpr bool isValid() const { return m_status != Status::missed; }
 
+  /// Returns the position of the interseciton
   constexpr const Position& position() const { return m_position; }
 
+  /// Returns the path length to the interseciton
   constexpr ActsScalar pathLength() const { return m_pathLength; }
 
+  /// Returns the intersection status enum
   constexpr Status status() const { return m_status; }
 
+  /// Static factory to creae an invalid instesection
   constexpr static Intersection invalid() { return Intersection(); }
 
   /// Comparison function for path length order i.e. intersection closest to
@@ -155,27 +159,35 @@ class ObjectIntersection {
   /// Returns whether the intersection was successful or not
   constexpr bool isValid() const { return m_intersection.isValid(); }
 
+  /// Returns the intersection
   constexpr const Intersection3D& intersection() const {
     return m_intersection;
   }
 
+  /// Returns the position of the interseciton
   constexpr const Intersection3D::Position& position() const {
     return m_intersection.position();
   }
 
+  /// Returns the path length to the interseciton
   constexpr ActsScalar pathLength() const {
     return m_intersection.pathLength();
   }
 
+  /// Returns the status of the interseciton
   constexpr Intersection3D::Status status() const {
     return m_intersection.status();
   }
 
+  /// Returns the object that has been intersected
   constexpr const object_t* object() const { return m_object; }
 
   constexpr std::uint8_t index() const { return m_index; }
 
-  constexpr static ObjectIntersection invalid() { return ObjectIntersection(); }
+  constexpr static ObjectIntersection invalid(
+      const object_t* object = nullptr) {
+    return ObjectIntersection(Intersection3D::invalid(), object);
+  }
 
   constexpr static bool pathLengthOrder(
       const ObjectIntersection& aIntersection,

Core/src/Navigation/CMakeLists.txt

@@ -0,0 +1 @@
+target_sources(ActsCore PRIVATE NavigationStream.cpp)

Core/src/Navigation/NavigationStream.cpp

@@ -0,0 +1,153 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "Acts/Navigation/NavigationStream.hpp"
+
+#include "Acts/Detector/Portal.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+
+#include <algorithm>
+
+bool Acts::NavigationStream::initialize(const GeometryContext& gctx,
+                                        const QueryPoint& queryPoint,
+                                        const BoundaryTolerance& cTolerance,
+                                        ActsScalar onSurfaceTolerance) {
+  // Position and direction from the query point
+  const Vector3& position = queryPoint.position;
+  const Vector3& direction = queryPoint.direction;
+
+  // De-duplicate first (necessary to deal correctly with multiple
+  // intersections) - sort them by surface pointer
+  std::ranges::sort(m_candidates, [](const Candidate& a, const Candidate& b) {
+    return (&a.surface()) < (&b.surface());
+  });
+  // Remove duplicates on basis of the surface pointer
+  m_candidates.erase(std::unique(m_candidates.begin(), m_candidates.end(),
+                                 [](const Candidate& a, const Candidate& b) {
+                                   return (&a.surface()) == (&b.surface());
+                                 }),
+                     m_candidates.end());
+
+  // A container collecting additional candidates from multiple
+  // valid interseciton
+  std::vector<Candidate> additionalCandidates = {};
+  for (auto& [sIntersection, portal, bTolerance] : m_candidates) {
+    // Get the surface from the object intersection
+    const Surface* surface = sIntersection.object();
+    // Intersect the surface
+    auto multiIntersection = surface->intersect(gctx, position, direction,
+                                                cTolerance, onSurfaceTolerance);
+
+    // Split them into valid intersections, keep track of potentially
+    // additional candidates
+    bool originalCandidateUpdated = false;
+    for (const auto& rsIntersection : multiIntersection.split()) {
+      // Skip negative solutions, respecting the on surface tolerance
+      if (rsIntersection.pathLength() < -onSurfaceTolerance) {
+        continue;
+      }
+      // Valid solution is either on surface or updates the distance
+      if (rsIntersection.isValid()) {
+        if (!originalCandidateUpdated) {
+          sIntersection = rsIntersection;
+          originalCandidateUpdated = true;
+        } else {
+          additionalCandidates.emplace_back(
+              Candidate{rsIntersection, portal, bTolerance});
+        }
+      }
+    }
+  }
+
+  // Append the multi intersection candidates
+  m_candidates.insert(m_candidates.end(), additionalCandidates.begin(),
+                      additionalCandidates.end());
+
+  // Sort the candidates by path length
+  std::ranges::sort(m_candidates, Candidate::pathLengthOrder);
+
+  // The we find the first invalid candidate
+  auto firstInvalid =
+      std::ranges::find_if(m_candidates, [](const Candidate& a) {
+        const auto& [aIntersection, aPortal, aTolerance] = a;
+        return !aIntersection.isValid();
+      });
+
+  // Set the range and initialize
+  m_candidates.resize(std::distance(m_candidates.begin(), firstInvalid));
+
+  m_currentIndex = 0;
+  if (m_candidates.empty()) {
+    return false;
+  }
+  return true;
+}
+
+bool Acts::NavigationStream::update(const GeometryContext& gctx,
+                                    const QueryPoint& queryPoint,
+                                    ActsScalar onSurfaceTolerance) {
+  // Position and direction from the query point
+  const Vector3& position = queryPoint.position;
+  const Vector3& direction = queryPoint.direction;
+
+  // Loop over the (currently valid) candidates and update
+  for (; m_currentIndex < m_candidates.size(); ++m_currentIndex) {
+    // Get the candidate, and resolve the tuple
+    Candidate& candidate = currentCandidate();
+    auto& [sIntersection, portal, bTolerance] = candidate;
+    // Get the surface from the object intersection
+    const Surface* surface = sIntersection.object();
+    // (re-)Intersect the surface
+    auto multiIntersection = surface->intersect(gctx, position, direction,
+                                                bTolerance, onSurfaceTolerance);
+    // Split them into valid intersections
+    for (const auto& rsIntersection : multiIntersection.split()) {
+      // Skip wrong index solution
+      if (rsIntersection.index() != sIntersection.index()) {
+        continue;
+      }
+      // Valid solution is either on surface or updates the distance
+      if (rsIntersection.isValid()) {
+        sIntersection = rsIntersection;
+        return true;
+      }
+    }
+  }
+  // No candidate was reachable
+  return false;
+}
+
+void Acts::NavigationStream::addSurfaceCandidate(
+    const Surface* surface, const BoundaryTolerance& bTolerance) {
+  m_candidates.emplace_back(Candidate{
+      ObjectIntersection<Surface>::invalid(surface), nullptr, bTolerance});
+}
+
+void Acts::NavigationStream::addSurfaceCandidates(
+    const std::vector<const Surface*>& surfaces,
+    const BoundaryTolerance& bTolerance) {
+  std::ranges::for_each(surfaces, [&](const auto* surface) {
+    m_candidates.emplace_back(Candidate{
+        ObjectIntersection<Surface>::invalid(surface), nullptr, bTolerance});
+  });
+}
+
+void Acts::NavigationStream::addPortalCandidate(const Portal* portal) {
+  m_candidates.emplace_back(
+      Candidate{ObjectIntersection<Surface>::invalid(&(portal->surface())),
+                portal, BoundaryTolerance::None()});
+}
+
+void Acts::NavigationStream::addPortalCandidates(
+    const std::vector<const Portal*>& portals) {
+  std::ranges::for_each(portals, [&](const auto& portal) {
+    m_candidates.emplace_back(
+        Candidate{ObjectIntersection<Surface>::invalid(&(portal->surface())),
+                  portal, BoundaryTolerance::None()});
+  });
+}