navigation/detail/intersection_kernel.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 // Project include(s)
0012 #include "detray/definitions/algebra.hpp"
0013 #include "detray/definitions/algorithms.hpp"
0014 #include "detray/definitions/detail/qualifiers.hpp"
0015 #include "detray/definitions/units.hpp"
0016 #include "detray/geometry/concepts.hpp"
0017 #include "detray/navigation/intersection/intersection.hpp"
0018 #include "detray/navigation/intersection/intersection_config.hpp"
0019 #include "detray/tracks/ray.hpp"
0020 #include "detray/utils/ranges.hpp"
0021
0022 namespace detray::detail {
0023
0024 /// A functor to add all valid intersections between the trajectory and surface
0025 template <template <typename, typename, bool> class intersector_t>
0026 struct intersection_initialize {
0027   /// Operator function to initialize intersections
0028   ///
0029   /// @tparam mask_group_t is the input mask group type found by variadic
0030   /// unrolling
0031   /// @tparam is_container_t is the intersection container type
0032   /// @tparam traj_t is the input trajectory type (e.g. ray or helix)
0033   /// @tparam surface_t is the input surface type
0034   /// @tparam transform_container_t is the input transform store type
0035   ///
0036   /// @param mask_group is the input mask group
0037   /// @param is_container is the intersection container to be filled
0038   /// @param traj is the input trajectory
0039   /// @param surface is the input surface
0040   /// @param contextual_transforms is the input transform container
0041   /// @param mask_tolerance is the tolerance for mask size
0042   /// @param overstep_tol negative cutoff for the path
0043   ///
0044   /// @return the number of valid intersections
0045   template <typename mask_group_t, typename mask_range_t,
0046             typename is_container_t, typename traj_t, typename surface_t,
0047             typename transform_container_t, concepts::scalar scalar_t>
0048   DETRAY_HOST_DEVICE inline void operator()(
0049       const mask_group_t &mask_group, const mask_range_t &mask_range,
0050       is_container_t &is_container, const traj_t &traj,
0051       const surface_t &sf_desc,
0052       const transform_container_t &contextual_transforms,
0053       const typename transform_container_t::context_type &ctx,
0054       const intersection::config &cfg,
0055       const scalar_t external_mask_tolerance = 0.f) const {
0056     using mask_t = typename mask_group_t::value_type;
0057     using shape_t = typename mask_t::shape;
0058     using algebra_t = typename mask_t::algebra_type;
0059     using intersection_t = typename is_container_t::value_type;
0060
0061     // Find the point of intersection with the underlying geometry
0062     const auto &ctf = contextual_transforms.at(sf_desc.transform(), ctx);
0063
0064     constexpr intersector_t<shape_t, algebra_t, intersection_t::contains_pos()>
0065         intersector{};
0066
0067     constexpr std::uint8_t n_sol{decltype(intersector)::n_solutions};
0068
0069     typename decltype(intersector)::result_type result{};
0070
0071     if constexpr (concepts::cylindrical<mask_t>) {
0072       std::size_t mask_idx{detail::invalid_value<std::size_t>()};
0073       if constexpr (concepts::interval<mask_range_t>) {
0074         mask_idx = mask_range.lower();
0075       } else {
0076         mask_idx = mask_range;
0077       }
0078       assert(mask_idx < mask_group.size());
0079
0080       result = intersector.point_of_intersection(
0081           traj, ctf, mask_group[mask_idx], cfg.overstep_tolerance);
0082     } else {
0083       result =
0084           intersector.point_of_intersection(traj, ctf, cfg.overstep_tolerance);
0085     }
0086
0087     // Check if any valid solutions were found
0088     if constexpr (n_sol > 1) {
0089       bool found_any{false};
0090       for (const auto &ip : result) {
0091         if (ip.is_valid()) {
0092           found_any = true;
0093         }
0094       }
0095       if (!found_any) [[unlikely]] {
0096         return;
0097       }
0098     } else {
0099       if (!result.is_valid()) [[unlikely]] {
0100         return;
0101       }
0102     }
0103
0104     // Resolve the masks that belong to the surface
0105     for (const auto &mask : detray::ranges::subrange(mask_group, mask_range)) {
0106       intersection_t is{};
0107
0108       // Build the resulting intersecion(s) from the intersection point
0109       if constexpr (n_sol > 1) {
0110         std::uint8_t n_found{0u};
0111
0112         for (std::size_t i = 0u; i < n_sol; ++i) {
0113           resolve_mask(is, traj, result[i], sf_desc, mask, ctf, cfg,
0114                        external_mask_tolerance);
0115
0116           if (is.is_probably_inside()) {
0117             insert_sorted(is, is_container);
0118             ++n_found;
0119           }
0120           if (n_found == n_sol) {
0121             return;
0122           }
0123         }
0124       } else {
0125         resolve_mask(is, traj, result, sf_desc, mask, ctf, cfg,
0126                      external_mask_tolerance);
0127
0128         if (is.is_probably_inside()) {
0129           insert_sorted(is, is_container);
0130           return;
0131         }
0132       }
0133     }
0134   }
0135
0136   template <typename intersection_t>
0137   DETRAY_HOST_DEVICE void insert_sorted(
0138       const intersection_t &sfi,
0139       std::vector<intersection_t> &intersections) const {
0140     auto itr_pos =
0141         detray::upper_bound(intersections.cbegin(), intersections.cend(), sfi);
0142
0143     intersections.insert(itr_pos, sfi);
0144   }
0145
0146   /// Specialization for the navigation state cache
0147   template <typename nav_state_t>
0148   DETRAY_HOST_DEVICE void insert_sorted(
0149       const typename nav_state_t::value_type &sfi,
0150       nav_state_t &intersections) const {
0151     auto itr_pos{intersections.cbegin()};
0152
0153     // For just two candidates int the cache, the navigation state keeps
0154     // the first as the previously visited candidate -> no sorting needed
0155     if constexpr (nav_state_t::capacity() > 2u) {
0156       itr_pos = detray::upper_bound(intersections.cbegin(),
0157                                     intersections.cend(), sfi);
0158     }
0159
0160     intersections.insert(itr_pos, sfi);
0161   }
0162 };
0163
0164 /// A functor to update the closest intersection between the trajectory and
0165 /// surface
0166 template <template <typename, typename, bool> class intersector_t>
0167 struct intersection_update {
0168   /// Operator function to update the intersection
0169   ///
0170   /// @tparam mask_group_t is the input mask group type found by variadic
0171   /// unrolling
0172   /// @tparam traj_t is the input trajectory type (e.g. ray or helix)
0173   /// @tparam surface_t is the input surface type
0174   /// @tparam transform_container_t is the input transform store type
0175   ///
0176   /// @param mask_group is the input mask group
0177   /// @param mask_range is the range of masks in the group that belong to the
0178   ///                   surface
0179   /// @param traj is the input trajectory
0180   /// @param surface is the input surface
0181   /// @param contextual_transforms is the input transform container
0182   /// @param mask_tolerance is the tolerance for mask size
0183   /// @param overstep_tol negative cutoff for the path
0184   ///
0185   /// @return the intersection
0186   template <typename mask_group_t, typename mask_range_t, typename traj_t,
0187             typename intersection_t, typename transform_container_t,
0188             concepts::scalar scalar_t>
0189   DETRAY_HOST_DEVICE inline bool operator()(
0190       const mask_group_t &mask_group, const mask_range_t &mask_range,
0191       const traj_t &traj, intersection_t &sfi,
0192       const transform_container_t &contextual_transforms,
0193       const typename transform_container_t::context_type &ctx,
0194       const intersection::config &cfg,
0195       const scalar_t external_mask_tolerance = 0.f) const {
0196     using mask_t = typename mask_group_t::value_type;
0197     using shape_t = typename mask_t::shape;
0198     using algebra_t = typename mask_t::algebra_type;
0199
0200     // Find the point of intersection with the underlying geometry
0201     const auto &ctf = contextual_transforms.at(sfi.surface().transform(), ctx);
0202
0203     constexpr intersector_t<shape_t, algebra_t, intersection_t::contains_pos()>
0204         intersector{};
0205     constexpr std::uint8_t n_sol{decltype(intersector)::n_solutions};
0206
0207     typename decltype(intersector)::result_type result{};
0208
0209     if constexpr (concepts::cylindrical<mask_t>) {
0210       std::size_t mask_idx{detail::invalid_value<std::size_t>()};
0211       if constexpr (concepts::interval<mask_range_t>) {
0212         mask_idx = mask_range.lower();
0213       } else {
0214         mask_idx = mask_range;
0215       }
0216       assert(mask_idx < mask_group.size());
0217
0218       result = intersector.point_of_intersection(
0219           traj, ctf, mask_group[mask_idx], cfg.overstep_tolerance);
0220     } else {
0221       result =
0222           intersector.point_of_intersection(traj, ctf, cfg.overstep_tolerance);
0223     }
0224
0225     // Check if any valid solutions were found
0226     if constexpr (n_sol > 1) {
0227       bool found_any{false};
0228       for (const auto &ip : result) {
0229         if (ip.is_valid()) {
0230           found_any = true;
0231         }
0232       }
0233       if (!found_any) [[unlikely]] {
0234         return false;
0235       }
0236     } else {
0237       if (!result.is_valid()) [[unlikely]] {
0238         return false;
0239       }
0240     }
0241
0242     // Run over the masks that belong to the surface
0243     for (const auto &mask : detray::ranges::subrange(mask_group, mask_range)) {
0244       // Build the resulting intersecion(s) from the intersection point
0245       if constexpr (n_sol > 1) {
0246         resolve_mask(sfi, traj, result[0], sfi.surface(), mask, ctf, cfg,
0247                      external_mask_tolerance);
0248       } else {
0249         resolve_mask(sfi, traj, result, sfi.surface(), mask, ctf, cfg,
0250                      external_mask_tolerance);
0251       }
0252
0253       if (sfi.is_probably_inside()) {
0254         return true;
0255       }
0256     }
0257
0258     return false;
0259   }
0260 };
0261
0262 }  // namespace detray::detail