detray/utils/bounding_volume.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/detail/qualifiers.hpp"
0014 #include "detray/definitions/units.hpp"
0015 #include "detray/geometry/mask.hpp"
0016 #include "detray/geometry/shapes/cuboid3D.hpp"
0017 #include "detray/geometry/shapes/cylinder3D.hpp"
0018 #include "detray/navigation/intersection/bounding_box/cuboid_intersector.hpp"
0019 #include "detray/navigation/intersection/intersection.hpp"
0020 #include "detray/tracks/ray.hpp"
0021
0022 // System include(s)
0023 #include <cassert>
0024 #include <limits>
0025 #include <type_traits>
0026 #include <vector>
0027
0028 namespace detray {
0029
0030 /// An axis aligned bounding box of a given @tparam shape_t
0031 template <typename shape_t, concepts::algebra algebra_t>
0032 class axis_aligned_bounding_volume {
0033  public:
0034   /// Define geometric properties
0035   /// @{
0036   using scalar_t = dscalar<algebra_t>;
0037   using shape = shape_t;
0038   using boundaries = typename shape_t::boundaries;
0039   template <typename A>
0040   using local_frame = typename shape_t::template local_frame_type<A>;
0041
0042   static constexpr std::size_t dim{shape_t::dim};
0043   /// @}
0044
0045   /// Default constructor builds an infinite box
0046   constexpr axis_aligned_bounding_volume() = default;
0047
0048   /// Constructor from mask boundary values
0049   template <typename... Args>
0050   DETRAY_HOST_DEVICE explicit constexpr axis_aligned_bounding_volume(
0051       std::size_t box_id, Args&&... args)
0052       : m_mask(box_id, std::forward<Args>(args)...) {}
0053
0054   /// Construct around an arbitrary surface @param mask
0055   template <typename mask_t, typename S = shape_t>
0056     requires std::is_same_v<S, cuboid3D>
0057   DETRAY_HOST_DEVICE constexpr axis_aligned_bounding_volume(
0058       const mask_t& mask, std::size_t box_id, const scalar_t envelope)
0059       : m_mask{mask.local_min_bounds(envelope).values(), box_id} {
0060     // Make sure the box is actually 'bounding'
0061     assert(envelope >= std::numeric_limits<scalar_t>::epsilon());
0062   }
0063
0064   /// Construct from mask boundary vector
0065   DETRAY_HOST axis_aligned_bounding_volume(const std::vector<scalar_t>& values,
0066                                            std::size_t box_id)
0067       : m_mask(values, box_id) {
0068     assert(values.size() == shape::boundaries::e_size &&
0069            " Given number of boundaries does not match mask shape.");
0070   }
0071
0072   /// Construct a bounding box around a set of boxes
0073   /// @note the given bounding volumes need to be defnined in the same
0074   /// local coordinate system!
0075   template <typename other_shape_t>
0076     requires std::is_same_v<
0077         typename shape::template local_frame_type<algebra_t>,
0078         typename other_shape_t::template local_frame_type<algebra_t>>
0079   DETRAY_HOST constexpr axis_aligned_bounding_volume(
0080       const std::vector<axis_aligned_bounding_volume<other_shape_t, algebra_t>>&
0081           aabbs,
0082       std::size_t box_id, const scalar_t env) {
0083     using loc_point_t = darray<scalar_t, other_shape_t::dim>;
0084
0085     // Find min/max extent of the local aabb in local coordinates
0086     constexpr scalar_t inv{detail::invalid_value<scalar_t>()};
0087     scalar_t min_x{inv};
0088     scalar_t min_y{inv};
0089     scalar_t min_z{inv};
0090     scalar_t max_x{-inv};
0091     scalar_t max_y{-inv};
0092     scalar_t max_z{-inv};
0093     for (const auto& vol : aabbs) {
0094       const auto min_point = vol.template loc_min<loc_point_t>();
0095       const auto max_point = vol.template loc_max<loc_point_t>();
0096
0097       // Check every coordinate of the points
0098       min_x = min_point[0] < min_x ? min_point[0] : min_x;
0099       min_y = min_point[1] < min_y ? min_point[1] : min_y;
0100
0101       max_x = max_point[0] > max_x ? max_point[0] : max_x;
0102       max_y = max_point[1] > max_y ? max_point[1] : max_y;
0103
0104       if constexpr (min_point.size() > 2) {
0105         min_z = min_point[2] < min_z ? min_point[2] : min_z;
0106         max_z = max_point[2] > max_z ? max_point[2] : max_z;
0107       }
0108     }
0109     m_mask = mask<shape, algebra_t, std::size_t>{
0110         box_id,      min_x - env, min_y - env, min_z - env,
0111         max_x + env, max_y + env, max_z + env};
0112   }
0113
0114   /// Subscript operator @returns a single box boundary.
0115   DETRAY_HOST_DEVICE
0116   constexpr auto operator[](const std::size_t i) const -> scalar_t {
0117     return m_mask[i];
0118   }
0119
0120   /// @returns the bounds of the box, depending on its shape
0121   DETRAY_HOST_DEVICE
0122   constexpr auto id() const -> std::size_t { return m_mask.volume_link(); }
0123
0124   /// @returns the bounds of the box, depending on its shape
0125   DETRAY_HOST_DEVICE
0126   constexpr auto bounds() const -> const mask<shape, algebra_t, std::size_t>& {
0127     return m_mask;
0128   }
0129
0130   /// @returns the minimum bounds of the volume in local coordinates
0131   template <typename point_t>
0132   DETRAY_HOST_DEVICE constexpr auto loc_min() const -> point_t {
0133     if constexpr (std::is_same_v<shape, cuboid3D>) {
0134       return {m_mask[cuboid3D::e_min_x], m_mask[cuboid3D::e_min_y],
0135               m_mask[cuboid3D::e_min_z]};
0136     } else if constexpr (std::is_same_v<shape, cylinder3D>) {
0137       return {-m_mask[cylinder3D::e_max_r], m_mask[cylinder3D::e_min_phi],
0138               m_mask[cylinder3D::e_min_z]};
0139     }
0140
0141     // If the volume shape is not supported, return universal minimum
0142     assert(false);
0143     constexpr scalar_t inv{detail::invalid_value<scalar_t>()};
0144     return point_t{-inv, -inv, -inv};
0145   }
0146
0147   /// @returns the maximum bounds of the volume in local coordinates
0148   template <typename point_t>
0149   DETRAY_HOST_DEVICE constexpr auto loc_max() const -> point_t {
0150     if constexpr (std::is_same_v<shape, cuboid3D>) {
0151       return {m_mask[cuboid3D::e_max_x], m_mask[cuboid3D::e_max_y],
0152               m_mask[cuboid3D::e_max_z]};
0153     } else if constexpr (std::is_same_v<shape, cylinder3D>) {
0154       return {m_mask[cylinder3D::e_max_r], m_mask[cylinder3D::e_max_phi],
0155               m_mask[cylinder3D::e_max_z]};
0156     }
0157
0158     // If the volume shape is not supported, return universal minimum
0159     // (or compilation error for 2D point)
0160     assert(false);
0161     constexpr scalar_t inv{detail::invalid_value<scalar_t>()};
0162     return point_t{inv, inv, inv};
0163   }
0164   /// @returns the minimum bounds of the volume in global cartesian
0165   /// coordinates
0166   template <concepts::transform3D transform3_t>
0167   DETRAY_HOST_DEVICE constexpr auto glob_min(const transform3_t& trf) const ->
0168       typename transform3_t::point3 {
0169     using point3_t = typename transform3_t::point3;
0170
0171     if constexpr (std::is_same_v<shape, cuboid3D>) {
0172       return trf.point_to_global(loc_min<point3_t>());
0173     } else if constexpr (std::is_same_v<shape, cylinder3D>) {
0174       return cylindrical3D<transform3_t>{}.local_to_global(trf, m_mask,
0175                                                            loc_min<point3_t>());
0176     }
0177
0178     // If the volume shape is not supported, return universal minimum
0179     // (or compilation error for 2D point)
0180     assert(false);
0181     constexpr scalar_t inv{detail::invalid_value<scalar_t>()};
0182     return point3_t{-inv, -inv, -inv};
0183   }
0184
0185   /// @returns the maximum bounds of the volume in global cartesian
0186   /// coordinates
0187   template <concepts::transform3D transform3_t>
0188   DETRAY_HOST_DEVICE constexpr auto glob_max(const transform3_t& trf) const ->
0189       typename transform3_t::point3 {
0190     using point3_t = typename transform3_t::point3;
0191
0192     if constexpr (std::is_same_v<shape, cuboid3D>) {
0193       return trf.point_to_global(loc_max<point3_t>());
0194     } else if constexpr (std::is_same_v<shape, cylinder3D>) {
0195       return cylindrical3D<transform3_t>{}.local_to_global(trf, m_mask,
0196                                                            loc_max<point3_t>());
0197     }
0198
0199     // If the volume shape is not supported, return universal minimum
0200     assert(false);
0201     constexpr scalar_t inv{detail::invalid_value<scalar_t>()};
0202     return point3_t{inv, inv, inv};
0203   }
0204
0205   /// @returns the geometric center position in global cartesian system
0206   template <concepts::point3D point3_t>
0207   DETRAY_HOST_DEVICE constexpr auto center() const -> point3_t {
0208     const scalar_t center_x{
0209         0.5f * (m_mask[cuboid3D::e_max_x] + m_mask[cuboid3D::e_min_x])};
0210     const scalar_t center_y{
0211         0.5f * (m_mask[cuboid3D::e_max_y] + m_mask[cuboid3D::e_min_y])};
0212     const scalar_t center_z{
0213         0.5f * (m_mask[cuboid3D::e_max_z] + m_mask[cuboid3D::e_min_z])};
0214
0215     return {detail::is_invalid_value(center_x) ? 0.f : center_x,
0216             detail::is_invalid_value(center_y) ? 0.f : center_y,
0217             detail::is_invalid_value(center_z) ? 0.f : center_z};
0218   }
0219
0220   /// @brief Lower and upper point for minimum axis aligned bounding box of
0221   /// cuboid shape.
0222   ///
0223   /// Computes the min and max vertices in global coordinates from the local
0224   /// minimum aabb. The global aabb is not necessarily an minimum aabb.
0225   ///
0226   /// @param trf affine transformation
0227   ///
0228   /// @returns a new, transformed aabb.
0229   template <concepts::transform3D transform3_t, typename S = shape_t>
0230     requires std::is_same_v<S, cuboid3D>
0231   DETRAY_HOST_DEVICE auto transform(const transform3_t& trf) const
0232       -> axis_aligned_bounding_volume {
0233     using point3_t = typename transform3_t::point3;
0234
0235     const scalar_t scalor_x{
0236         (m_mask[cuboid3D::e_max_x] - m_mask[cuboid3D::e_min_x])};
0237     const scalar_t scalor_y{
0238         (m_mask[cuboid3D::e_max_y] - m_mask[cuboid3D::e_min_y])};
0239     const scalar_t scalor_z{
0240         (m_mask[cuboid3D::e_max_z] - m_mask[cuboid3D::e_min_z])};
0241
0242     // Cannot handle 'inv' propagation through the calculation for now
0243     if (detail::is_invalid_value(scalor_x) ||
0244         detail::is_invalid_value(scalor_y) ||
0245         detail::is_invalid_value(scalor_z)) {
0246       // If the box was infinite to begin with, it stays that way
0247       assert(detail::is_invalid_value(scalor_x) &&
0248              detail::is_invalid_value(scalor_y) &&
0249              detail::is_invalid_value(scalor_z));
0250
0251       return *this;
0252     }
0253
0254     // The new axis vectors, scaled to the aabb dimensions
0255     // (e.g. max_x - min_x)
0256     const point3_t new_box_x = scalor_x * trf.x();
0257     const point3_t new_box_y = scalor_y * trf.y();
0258     const point3_t new_box_z = scalor_z * trf.z();
0259
0260     // Transform the old min and max points to the global frame and
0261     // construct all corner points of the local aabb in global coordinates
0262     darray<point3_t, 8> glob_c_points;
0263     glob_c_points[0] = glob_min(trf);
0264     glob_c_points[1] = glob_max(trf);
0265     glob_c_points[2] = glob_c_points[0] + new_box_x;
0266     glob_c_points[3] = glob_c_points[0] + new_box_y;
0267     glob_c_points[4] = glob_c_points[0] + new_box_z;
0268     glob_c_points[5] = glob_c_points[1] - new_box_x;
0269     glob_c_points[6] = glob_c_points[1] - new_box_y;
0270     glob_c_points[7] = glob_c_points[1] - new_box_z;
0271
0272     // Find min/max extent of the local aabb in global coordinates
0273     constexpr scalar_t inv{detail::invalid_value<scalar_t>()};
0274     scalar_t min_x{inv};
0275     scalar_t min_y{inv};
0276     scalar_t min_z{inv};
0277     scalar_t max_x{-inv};
0278     scalar_t max_y{-inv};
0279     scalar_t max_z{-inv};
0280     for (const point3_t& p : glob_c_points) {
0281       // Check every coordinate of the point
0282       min_x = p[0] < min_x ? p[0] : min_x;
0283       max_x = p[0] > max_x ? p[0] : max_x;
0284       min_y = p[1] < min_y ? p[1] : min_y;
0285       max_y = p[1] > max_y ? p[1] : max_y;
0286       min_z = p[2] < min_z ? p[2] : min_z;
0287       max_z = p[2] > max_z ? p[2] : max_z;
0288     }
0289
0290     // Construct the transformed aabb
0291     return axis_aligned_bounding_volume{
0292         m_mask.volume_link(), min_x, min_y, min_z, max_x, max_y, max_z};
0293   }
0294
0295   /// Checks whether a point lies inside the box. The point has to be defined
0296   /// in the coordinate frame that is spanned by the box axes.
0297   template <concepts::point point_t>
0298   DETRAY_HOST_DEVICE constexpr bool is_inside(
0299       const point_t& loc_p,
0300       const scalar_t t = std::numeric_limits<scalar_t>::epsilon()) const {
0301     return m_mask.is_inside(loc_p, t);
0302   }
0303
0304   /// Intersect the box with a ray
0305   DETRAY_HOST_DEVICE
0306   constexpr bool intersect(
0307       const detail::ray<algebra_t>& ray,
0308       const scalar_t t = std::numeric_limits<scalar_t>::epsilon()) const {
0309     static_assert(std::is_same_v<shape, cuboid3D>,
0310                   "aabbs are only implemented in cuboid shape for now");
0311     return cuboid_intersector{}(ray, m_mask, t);
0312   }
0313
0314  private:
0315   /// Print the bounding volume
0316   DETRAY_HOST friend std::ostream& operator<<(
0317       std::ostream& os, const axis_aligned_bounding_volume& aabb) {
0318     return os << aabb.bounds().to_string();
0319   }
0320   /// Keeps the box boundary values and id
0321   mask<shape, algebra_t, std::size_t> m_mask;
0322 };
0323
0324 }  // namespace detray