geometry/shapes/trapezoid2D.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/geometry/coordinates/cartesian2D.hpp"
0018 #include "detray/geometry/coordinates/cartesian3D.hpp"
0019 #include "detray/geometry/detail/shape_utils.hpp"
0020
0021 // System include(s)
0022 #include <limits>
0023 #include <ostream>
0024 #include <string_view>
0025
0026 namespace detray {
0027
0028 /// @brief Geometrical shape of a trapezoid2D.
0029 ///
0030 /// It is defined by half lengths in local0 coordinate bounds[0] and bounds[1]
0031 /// at -/+ half length in the local1 coordinate bounds[2]. bounds[3] contains
0032 /// the precomputed value of 1 / (2 * bounds[2]), which avoids
0033 /// excessive floating point divisions.
0034 class trapezoid2D {
0035  public:
0036   /// The name for this shape
0037   static constexpr std::string_view name = "trapezoid2D";
0038
0039   enum boundaries : unsigned int {
0040     e_half_length_0 = 0u,
0041     e_half_length_1 = 1u,
0042     e_half_length_2 = 2u,
0043     e_divisor = 3u,  // 1 / (2 * bounds[e_half_length_2])
0044     e_size = 4u,
0045   };
0046
0047   /// Container definition for the shape boundary values
0048   template <concepts::scalar scalar_t>
0049   using bounds_type = darray<scalar_t, boundaries::e_size>;
0050
0051   /// Local coordinate frame for boundary checks
0052   template <concepts::algebra algebra_t>
0053   using local_frame_type = cartesian2D<algebra_t>;
0054
0055   /// Result type of a boundary check
0056   template <typename bool_t>
0057   using result_type = detail::boundary_check_result<bool_t>;
0058
0059   /// Dimension of the local coordinate system
0060   static constexpr std::size_t dim{2u};
0061
0062   /// @brief Find the minimum distance to any boundary.
0063   ///
0064   /// @note the point is expected to be given in local coordinates by the
0065   /// caller.
0066   ///
0067   /// @param bounds the boundary values for this shape
0068   /// @param loc_p the point to be checked in the local coordinate system
0069   ///
0070   /// @return the minimum distance.
0071   template <concepts::scalar scalar_t, concepts::point point_t>
0072   DETRAY_HOST_DEVICE inline scalar_t min_dist_to_boundary(
0073       const bounds_type<scalar_t> &bounds, const point_t &loc_p) const {
0074     // Minimum distance between loc point and line between corner points
0075     const scalar_t d_y = -2.f * bounds[e_half_length_2];
0076     const scalar_t d_x = bounds[e_half_length_1] - bounds[e_half_length_0];
0077
0078     const scalar_t denom = math::sqrt(d_y * d_y + d_x * d_x);
0079     const scalar_t dist_to_line =
0080         math::fabs(d_y * loc_p[0] - d_x * loc_p[1] -
0081                    bounds[e_half_length_2] *
0082                        (bounds[e_half_length_1] + bounds[e_half_length_0])) /
0083         denom;
0084
0085     return math::min(dist_to_line, math::fabs(math::fabs(loc_p[1]) -
0086                                               bounds[e_half_length_2]));
0087   }
0088
0089   /// @brief Check boundary values for a local point.
0090   /// @{
0091   /// @param bounds the boundary values for this shape
0092   /// @param trf the surface transform
0093   /// @param glob_p the point to be checked in the global coordinate system
0094   /// @param tol dynamic tolerance determined by caller
0095   ///
0096   /// @return true if the local point lies within the given boundaries.
0097   template <concepts::algebra algebra_t>
0098   DETRAY_HOST_DEVICE constexpr result_type<dbool<algebra_t>> check_boundaries(
0099       const bounds_type<dscalar<algebra_t>> &bounds,
0100       const dtransform3D<algebra_t> &trf, const dpoint3D<algebra_t> &glob_p,
0101       const dscalar<algebra_t> tol =
0102           std::numeric_limits<dscalar<algebra_t>>::epsilon(),
0103       const dscalar<algebra_t> edge_tol = 0.f) const {
0104     // Get the full local position
0105     const dpoint3D<algebra_t> loc_p =
0106         cartesian3D<algebra_t>::global_to_local(trf, glob_p, {});
0107
0108     return check_boundaries(bounds, loc_p, tol, edge_tol);
0109   }
0110
0111   /// @note the point is expected to be given in local coordinates by the
0112   /// caller. For the conversion from global cartesian coordinates, the
0113   /// nested @c shape struct can be used.
0114   ///
0115   /// @param bounds the boundary values for this shape
0116   /// @param loc_p the point to be checked in the local coordinate system
0117   /// @param tol dynamic tolerance determined by caller
0118   ///
0119   /// @return true if the local point lies within the given boundaries.
0120   template <concepts::scalar scalar_t, concepts::point point_t>
0121   DETRAY_HOST_DEVICE constexpr auto check_boundaries(
0122       const bounds_type<scalar_t> &bounds, const point_t &loc_p,
0123       const scalar_t tol = std::numeric_limits<scalar_t>::epsilon(),
0124       const scalar_t edge_tol = 0.f) const {
0125     const scalar_t rel_y =
0126         (bounds[e_half_length_2] + loc_p[1]) * bounds[e_divisor];
0127     const scalar_t bound_x{
0128         bounds[e_half_length_0] +
0129         rel_y * (bounds[e_half_length_1] - bounds[e_half_length_0])};
0130
0131     auto inside_mask{(math::fabs(loc_p[0]) <= (bound_x + tol) &&
0132                       math::fabs(loc_p[1]) <= (bounds[e_half_length_2] + tol))};
0133
0134     decltype(inside_mask) inside_edge{false};
0135     if (detail::any_of(edge_tol > 0.f)) {
0136       const scalar_t full_tol{tol + edge_tol};
0137
0138       inside_edge =
0139           (math::fabs(loc_p[0]) <= (bound_x + full_tol) &&
0140            math::fabs(loc_p[1]) <= (bounds[e_half_length_2] + full_tol));
0141     }
0142
0143     return result_type<decltype(inside_mask)>{inside_mask, inside_edge};
0144   }
0145   /// @}
0146
0147   /// @brief Measure of the shape: Area
0148   ///
0149   /// @param bounds the boundary values for this shape
0150   ///
0151   /// @returns the trapezoid area on the plane
0152   template <concepts::scalar scalar_t>
0153   DETRAY_HOST_DEVICE constexpr scalar_t measure(
0154       const bounds_type<scalar_t> &bounds) const {
0155     return area(bounds);
0156   }
0157
0158   /// @brief The area of a the shape
0159   ///
0160   /// @param bounds the boundary values for this shape
0161   ///
0162   /// @returns the trapezoid area.
0163   template <concepts::scalar scalar_t>
0164   DETRAY_HOST_DEVICE constexpr scalar_t area(
0165       const bounds_type<scalar_t> &bounds) const {
0166     return 2.f * (bounds[e_half_length_0] + bounds[e_half_length_1]) *
0167            bounds[e_half_length_2];
0168   }
0169
0170   /// @brief Merge two trapezoid shapes
0171   ///
0172   /// @param bounds the boundary values for this shape
0173   /// @param o_bounds the boundary values for the other shape
0174   ///
0175   /// @returns merged bound values
0176   template <concepts::scalar scalar_t>
0177   DETRAY_HOST_DEVICE constexpr bounds_type<scalar_t> merge(
0178       const bounds_type<scalar_t> &bounds,
0179       const bounds_type<scalar_t> &o_bounds) const {
0180     bounds_type<scalar_t> new_bounds{};
0181
0182     new_bounds[e_half_length_0] =
0183         math::max(bounds[e_half_length_0], o_bounds[e_half_length_0]);
0184     new_bounds[e_half_length_1] =
0185         math::max(bounds[e_half_length_1], o_bounds[e_half_length_1]);
0186     new_bounds[e_half_length_2] =
0187         math::max(bounds[e_half_length_2], o_bounds[e_half_length_2]);
0188
0189     // Recalculate the memoized divisor
0190     new_bounds[e_divisor] = 1.f / (2.f * new_bounds[e_half_length_2]);
0191
0192     return new_bounds;
0193   }
0194
0195   /// @brief Lower and upper point for minimal axis aligned bounding box.
0196   ///
0197   /// Computes the min and max vertices in a local cartesian frame.
0198   ///
0199   /// @param bounds the boundary values for this shape
0200   /// @param env dynamic envelope around the shape
0201   ///
0202   /// @returns and array of coordinates that contains the lower point (first
0203   /// three values) and the upper point (latter three values) .
0204   template <concepts::algebra algebra_t>
0205   DETRAY_HOST_DEVICE inline darray<dscalar<algebra_t>, 6> local_min_bounds(
0206       const bounds_type<dscalar<algebra_t>> &bounds,
0207       const dscalar<algebra_t> env =
0208           std::numeric_limits<dscalar<algebra_t>>::epsilon()) const {
0209     using scalar_t = dscalar<algebra_t>;
0210
0211     assert(env > 0.f);
0212     const scalar_t x_bound{(bounds[e_half_length_0] > bounds[e_half_length_1]
0213                                 ? bounds[e_half_length_0]
0214                                 : bounds[e_half_length_1]) +
0215                            env};
0216     const scalar_t y_bound{bounds[e_half_length_2] + env};
0217     return {-x_bound, -y_bound, -env, x_bound, y_bound, env};
0218   }
0219
0220   /// @returns the shapes centroid in local cartesian coordinates
0221   template <concepts::algebra algebra_t>
0222   DETRAY_HOST_DEVICE dpoint3D<algebra_t> centroid(
0223       const bounds_type<dscalar<algebra_t>> &bounds) const {
0224     using scalar_t = dscalar<algebra_t>;
0225
0226     const scalar_t h_2{bounds[e_half_length_2]};
0227     const scalar_t a_2{bounds[e_half_length_1]};
0228     const scalar_t b_2{bounds[e_half_length_0]};
0229
0230     const scalar_t y{2.f * h_2 * (2.f * a_2 + b_2) * 1.f / (3.f * (a_2 + b_2))};
0231
0232     return {static_cast<scalar_t>(0), y - h_2, static_cast<scalar_t>(0)};
0233   }
0234
0235   /// Generate vertices in local cartesian frame
0236   ///
0237   /// @param bounds the boundary values for the trapezoid
0238   /// @param ls is the number of line segments
0239   ///
0240   /// @return a generated list of vertices
0241   template <concepts::algebra algebra_t>
0242   DETRAY_HOST dvector<dpoint3D<algebra_t>> vertices(
0243       const bounds_type<dscalar<algebra_t>> &bounds, dindex /*ignored*/) const {
0244     using point3_t = dpoint3D<algebra_t>;
0245
0246     // left hand lower corner
0247     point3_t lh_lc{-bounds[e_half_length_0], -bounds[e_half_length_2],
0248                    static_cast<dscalar<algebra_t>>(0.f)};
0249     // right hand lower corner
0250     point3_t rh_lc{bounds[e_half_length_0], -bounds[e_half_length_2],
0251                    static_cast<dscalar<algebra_t>>(0.f)};
0252     // right hand upper corner
0253     point3_t rh_uc{bounds[e_half_length_1], bounds[e_half_length_2],
0254                    static_cast<dscalar<algebra_t>>(0.f)};
0255     // left hand upper corner
0256     point3_t lh_uc{-bounds[e_half_length_1], bounds[e_half_length_2],
0257                    static_cast<dscalar<algebra_t>>(0.f)};
0258
0259     // Return the confining vertices
0260     return {lh_lc, rh_lc, rh_uc, lh_uc};
0261   }
0262
0263   /// @brief Check consistency of boundary values.
0264   ///
0265   /// @param bounds the boundary values for this shape
0266   /// @param os output stream for error messages
0267   ///
0268   /// @return true if the bounds are consistent.
0269   template <concepts::scalar scalar_t>
0270   DETRAY_HOST constexpr bool check_consistency(
0271       const bounds_type<scalar_t> &bounds, std::ostream &os) const {
0272     constexpr auto tol{10.f * std::numeric_limits<scalar_t>::epsilon()};
0273
0274     if (bounds[e_half_length_0] < tol || bounds[e_half_length_1] < tol) {
0275       os << "DETRAY ERROR (HOST): Half length in x must be in the range "
0276             "(0, "
0277             "numeric_max)"
0278          << std::endl;
0279       return false;
0280     }
0281     if (bounds[e_half_length_2] < tol) {
0282       os << "DETRAY ERROR (HOST): Half length in y must be in the range "
0283             "(0, "
0284             "numeric_max)"
0285          << std::endl;
0286       return false;
0287     }
0288
0289     if (const auto div{1.f / (2.f * bounds[e_half_length_2])};
0290         math::fabs(bounds[e_divisor] - div) > tol) {
0291       os << "DETRAY ERROR (HOST): Divisor incorrect. Should be: " << div
0292          << std::endl;
0293       return false;
0294     }
0295
0296     return true;
0297   }
0298 };
0299
0300 }  // namespace detray