geometry/shapes/cylinder2D.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/containers.hpp"
0014 #include "detray/definitions/detail/qualifiers.hpp"
0015 #include "detray/definitions/indexing.hpp"
0016 #include "detray/definitions/math.hpp"
0017 #include "detray/definitions/units.hpp"
0018 #include "detray/geometry/coordinates/cartesian3D.hpp"
0019 #include "detray/geometry/coordinates/cylindrical2D.hpp"
0020 #include "detray/geometry/detail/shape_utils.hpp"
0021
0022 // System include(s)
0023 #include <limits>
0024 #include <ostream>
0025 #include <string_view>
0026
0027 namespace detray {
0028
0029 /// @brief Geometrical shape of a 2D cylinder.
0030 ///
0031 /// It is defined by r and the two half lengths rel to the coordinate center.
0032 class cylinder2D {
0033  public:
0034   /// The name for this shape
0035   static constexpr std::string_view name = "cylinder2D";
0036
0037   enum boundaries : unsigned int {
0038     e_r = 0u,
0039     e_lower_z = 1u,
0040     e_upper_z = 2u,
0041     e_size = 3u,
0042   };
0043
0044   /// Container definition for the shape boundary values
0045   template <concepts::scalar scalar_t>
0046   using bounds_type = darray<scalar_t, boundaries::e_size>;
0047
0048   /// Local coordinate frame for boundary checks
0049   template <concepts::algebra algebra_t>
0050   using local_frame_type = cylindrical2D<algebra_t>;
0051
0052   /// Result type of a boundary check
0053   template <typename bool_t>
0054   using result_type = detail::boundary_check_result<bool_t>;
0055
0056   /// Dimension of the local coordinate system
0057   static constexpr std::size_t dim{2u};
0058
0059   /// @brief Find the minimum distance to any boundary.
0060   ///
0061   /// @note the point is expected to be given in local coordinates by the
0062   /// caller.
0063   ///
0064   /// @param bounds the boundary values for this shape
0065   /// @param loc_p the point to be checked in the local coordinate system
0066   ///
0067   /// @return the minimum distance.
0068   template <concepts::scalar scalar_t, concepts::point point_t>
0069   DETRAY_HOST_DEVICE inline scalar_t min_dist_to_boundary(
0070       const bounds_type<scalar_t> &bounds, const point_t &loc_p) const {
0071     return math::min(math::fabs(loc_p[1] - bounds[e_lower_z]),
0072                      math::fabs(bounds[e_upper_z] - loc_p[1]));
0073   }
0074
0075   /// @brief Check boundary values for a local point.
0076   /// @{
0077   /// @param bounds the boundary values for this shape
0078   /// @param trf the surface transform
0079   /// @param glob_p the point to be checked in the global coordinate system
0080   /// @param tol dynamic tolerance determined by caller
0081   ///
0082   /// @return true if the local point lies within the given boundaries.
0083   template <concepts::algebra algebra_t>
0084   DETRAY_HOST_DEVICE constexpr result_type<dbool<algebra_t>> check_boundaries(
0085       const bounds_type<dscalar<algebra_t>> &bounds,
0086       const dtransform3D<algebra_t> &trf, const dpoint3D<algebra_t> &glob_p,
0087       const dscalar<algebra_t> tol =
0088           std::numeric_limits<dscalar<algebra_t>>::epsilon(),
0089       const dscalar<algebra_t> edge_tol = 0.f) const {
0090     // Rotate to local cartesian
0091     const auto loc_p = cartesian3D<algebra_t>::global_to_local(trf, glob_p, {});
0092
0093     // Only need the z-position for the check
0094     return check_boundaries(
0095         bounds,
0096         dpoint2D<algebra_t>{static_cast<dscalar<algebra_t>>(0.f), loc_p[2]},
0097         tol, edge_tol);
0098   }
0099
0100   /// @note the point is expected to be given in local coordinates by the
0101   /// caller. For the conversion from global cartesian coordinates, the
0102   /// nested @c shape struct can be used. The point is assumed to be in
0103   /// the cylinder 2D frame (r * phi, z).
0104   ///
0105   /// @tparam is_rad_check whether the radial bound should be checked in this
0106   /// call.
0107   ///
0108   /// @param bounds the boundary values for this shape
0109   /// @param loc_p the point to be checked in the local coordinate system
0110   /// @param tol dynamic tolerance determined by caller
0111   ///
0112   /// @return true if the local point lies within the given boundaries.
0113   template <concepts::scalar scalar_t, concepts::point point_t>
0114   DETRAY_HOST_DEVICE constexpr auto check_boundaries(
0115       const bounds_type<scalar_t> &bounds, const point_t &loc_p,
0116       const scalar_t tol = std::numeric_limits<scalar_t>::epsilon(),
0117       const scalar_t edge_tol = 0.f) const {
0118     auto inside_mask{((bounds[e_lower_z] - tol) <= loc_p[1] &&
0119                       loc_p[1] <= (bounds[e_upper_z] + tol))};
0120
0121     decltype(inside_mask) inside_edge{false};
0122     if (detail::any_of(edge_tol > 0.f)) {
0123       const scalar_t full_tol{tol + edge_tol};
0124       inside_edge = ((bounds[e_lower_z] - full_tol) <= loc_p[1] &&
0125                      loc_p[1] <= (bounds[e_upper_z] + full_tol));
0126     }
0127
0128     return result_type<decltype(inside_mask)>{inside_mask, inside_edge};
0129   }
0130   /// @}
0131
0132   /// @brief Measure of the shape: Area
0133   ///
0134   /// @param bounds the boundary values for this shape
0135   ///
0136   /// @returns the cylinder area on the cylinder of radius r.
0137   template <concepts::scalar scalar_t>
0138   DETRAY_HOST_DEVICE constexpr scalar_t measure(
0139       const bounds_type<scalar_t> &bounds) const {
0140     return area(bounds);
0141   }
0142
0143   /// @brief The area of a the shape
0144   ///
0145   /// @param bounds the boundary values for this shape
0146   ///
0147   /// @returns the stereo annulus area.
0148   template <concepts::scalar scalar_t>
0149   DETRAY_HOST_DEVICE constexpr scalar_t area(
0150       const bounds_type<scalar_t> &bounds) const {
0151     return 2.f * constant<scalar_t>::pi * bounds[e_r] *
0152            (bounds[e_upper_z] - bounds[e_lower_z]);
0153   }
0154
0155   /// @brief Merge two 2D cylinder shapes
0156   ///
0157   /// @param bounds the boundary values for this shape
0158   /// @param o_bounds the boundary values for the other shape
0159   ///
0160   /// @returns merged bound values
0161   template <concepts::scalar scalar_t>
0162   DETRAY_HOST_DEVICE constexpr bounds_type<scalar_t> merge(
0163       const bounds_type<scalar_t> &bounds,
0164       const bounds_type<scalar_t> &o_bounds) const {
0165     bounds_type<scalar_t> new_bounds{};
0166
0167     new_bounds[e_r] = math::max(bounds[e_r], o_bounds[e_r]);
0168     new_bounds[e_lower_z] = math::min(bounds[e_lower_z], o_bounds[e_lower_z]);
0169     new_bounds[e_upper_z] = math::max(bounds[e_upper_z], o_bounds[e_upper_z]);
0170
0171     return new_bounds;
0172   }
0173
0174   /// @brief Lower and upper point for minimal axis aligned bounding box.
0175   ///
0176   /// Computes the min and max vertices in a local cartesian frame.
0177   ///
0178   /// @param bounds the boundary values for this shape
0179   /// @param env dynamic envelope around the shape
0180   ///
0181   /// @returns an array of coordinates that contains the lower point (first
0182   /// three values) and the upper point (latter three values) .
0183   template <concepts::algebra algebra_t>
0184   DETRAY_HOST_DEVICE inline darray<dscalar<algebra_t>, 6> local_min_bounds(
0185       const bounds_type<dscalar<algebra_t>> &bounds,
0186       const dscalar<algebra_t> env =
0187           std::numeric_limits<dscalar<algebra_t>>::epsilon()) const {
0188     assert(env > 0.f);
0189
0190     const dscalar<algebra_t> xy_bound{bounds[e_r] + env};
0191     return {-xy_bound, -xy_bound, bounds[e_lower_z] - env,
0192             xy_bound,  xy_bound,  bounds[e_upper_z] + env};
0193   }
0194
0195   /// @returns the shapes centroid in local cartesian coordinates
0196   template <concepts::algebra algebra_t>
0197   DETRAY_HOST_DEVICE dpoint3D<algebra_t> centroid(
0198       const bounds_type<dscalar<algebra_t>> &bounds) const {
0199     return {static_cast<dscalar<algebra_t>>(0.f),
0200             static_cast<dscalar<algebra_t>>(0.f),
0201             0.5f * (bounds[e_lower_z] + bounds[e_upper_z])};
0202   }
0203
0204   /// Generate vertices in local cartesian frame
0205   ///
0206   /// @param bounds the boundary values for the cylinder
0207   /// @param n_seg is the number of line segments
0208   ///
0209   /// @return a generated list of vertices
0210   template <concepts::algebra algebra_t>
0211   DETRAY_HOST dvector<dpoint3D<algebra_t>> vertices(
0212       const bounds_type<dscalar<algebra_t>> & /*unused*/,
0213       dindex /*unused*/) const {
0214     throw std::runtime_error(
0215         "Vertex generation for cylinders is not implemented");
0216     return {};
0217   }
0218
0219   /// @brief Check consistency of boundary values.
0220   ///
0221   /// @param bounds the boundary values for this shape
0222   /// @param os output stream for error messages
0223   ///
0224   /// @return true if the bounds are consistent.
0225   template <concepts::scalar scalar_t>
0226   DETRAY_HOST constexpr bool check_consistency(
0227       const bounds_type<scalar_t> &bounds, std::ostream &os) const {
0228     constexpr auto tol{10.f * std::numeric_limits<scalar_t>::epsilon()};
0229
0230     if (bounds[e_r] < tol) {
0231       os << "DETRAY ERROR (HOST): Radius must be in the range (0, "
0232             "numeric_max)"
0233          << std::endl;
0234       return false;
0235     }
0236     if (bounds[e_lower_z] >= bounds[e_upper_z] ||
0237         math::fabs(bounds[e_lower_z] - bounds[e_upper_z]) < tol) {
0238       os << "DETRAY ERROR (HOST): Neg. half length must be smaller than "
0239             "pos. half "
0240             "length.";
0241       return false;
0242     }
0243
0244     return true;
0245   }
0246 };
0247
0248 }  // namespace detray