ref: Move the vertexing to the mask and add proto surface conversion …

…for straw surfaces (#572)

Move the vertexing functionality to the masks and add support for plotting straw surfaces in svgtools

core/include/detray/geometry/detail/surface_kernels.hpp

@@ -43,7 +43,7 @@ struct surface_kernels {
         DETRAY_HOST_DEVICE inline auto operator()(
             const mask_group_t& mask_group, const index_t& index) const {
 
-            return mask_group.at(index).volume_link();
+            return mask_group[index].volume_link();
         }
     };
 
@@ -64,7 +64,7 @@ struct surface_kernels {
         DETRAY_HOST_DEVICE inline point3 operator()(
             const mask_group_t& mask_group, const index_t& index,
             const transform3& trf3, const point2& bound) const {
-            // Get the concrete mask instance for this surface
+
             const auto& m = mask_group[index];
 
             return m.local_frame().normal(trf3, bound, m);
@@ -114,6 +114,7 @@ struct surface_kernels {
             const mask_group_t& mask_group, const index_t& index,
             const transform3& trf3, const point2& bound,
             const vector3& dir) const {
+
             const auto& m = mask_group[index];
 
             return mask_group[index].local_frame().bound_local_to_global(
@@ -138,10 +139,9 @@ struct surface_kernels {
         DETRAY_HOST_DEVICE inline bound_vector_type operator()(
             const mask_group_t& mask_group, const index_t& index,
             const transform3& trf3, const free_vector_type& free_vec) const {
-            // Get the concrete mask instance for this surface
-            const auto& m = mask_group[index];
-            // Call shape-specific behaviour on the mask
-            return m.local_frame().free_to_bound_vector(trf3, free_vec);
+
+            return mask_group[index].local_frame().free_to_bound_vector(
+                trf3, free_vec);
         }
     };
 
@@ -167,9 +167,8 @@ struct surface_kernels {
             const mask_group_t& mask_group, const index_t& index,
             const transform3& trf3, const free_vector_type& free_vec) const {
 
-            const auto& m = mask_group[index];
-
-            return m.local_frame().free_to_bound_jacobian(trf3, free_vec);
+            return mask_group[index].local_frame().free_to_bound_jacobian(
+                trf3, free_vec);
         }
     };
 
@@ -196,9 +195,8 @@ struct surface_kernels {
             const transform3& trf3, const vector3& pos, const vector3& dir,
             const vector3& dtds) const {
 
-            const auto& m = mask_group[index];
-
-            return m.local_frame().path_correction(pos, dir, dtds, trf3);
+            return mask_group[index].local_frame().path_correction(pos, dir,
+                                                                   dtds, trf3);
         }
     };
 
@@ -211,9 +209,19 @@ struct surface_kernels {
             const scalar_t env =
                 std::numeric_limits<scalar_t>::epsilon()) const {
 
-            const auto& m = mask_group[index];
+            return mask_group[index].local_min_bounds(env);
+        }
+    };
+
+    /// A functor to get the vertices in local coordinates.
+    struct local_vertices {
+
+        template <typename mask_group_t, typename index_t, typename scalar_t>
+        DETRAY_HOST_DEVICE inline auto operator()(
+            const mask_group_t& mask_group, const index_t& index,
+            const dindex n_seg) const {
 
-            return m.local_min_bounds(env);
+            return mask_group[index].vertices(n_seg);
         }
     };
 };

core/include/detray/geometry/surface.hpp

@@ -243,17 +243,19 @@ class surface {
                                                              pos, dir, dtds);
     }
 
-    /// @returns the vertices in local frame
+    /// @returns the vertices in local frame with @param n_seg the number of
+    /// segments used along acrs
     DETRAY_HOST
-    constexpr auto local_vertices() const {
-        return visit_mask<typename kernels::local_vertices>();
+    constexpr auto local_vertices(const dindex n_seg) const {
+        return visit_mask<typename kernels::vertices>(n_seg);
     }
 
-    /// @returns the vertices in global frame.
+    /// @returns the vertices in global frame with @param n_seg the number of
+    /// segments used along acrs
     DETRAY_HOST
-    constexpr auto global_vertices(const context &ctx,
+    constexpr auto global_vertices(const context &ctx, const dindex n_seg,
                                    const vector3 &dir) const {
-        auto vertices = local_vertices();
+        auto vertices = local_vertices(n_seg);
         for (size_t i = 0; i < vertices.size(); i++) {
             vertices[i] = local_to_global(ctx, vertices[i], dir);
         }

core/include/detray/masks/annulus2D.hpp

@@ -14,11 +14,11 @@
 #include "detray/definitions/qualifiers.hpp"
 #include "detray/definitions/units.hpp"
 #include "detray/intersection/plane_intersector.hpp"
+#include "detray/masks/detail/vertexing.hpp"
 #include "detray/surface_finders/grid/detail/axis_binning.hpp"
 #include "detray/surface_finders/grid/detail/axis_bounds.hpp"
 
 // System include(s)
-#include <cmath>
 #include <limits>
 #include <ostream>
 #include <string>
@@ -93,24 +93,6 @@ class annulus2D {
         using binning = dtuple<binning_loc0<C, S>, binning_loc1<C, S>>;
     };
 
-    /// Given a local polar point given in the disc frame, @returns the
-    /// correponding point in the focal frame
-    /*DETRAY_HOST_DEVICE
-    template<typename scalar_t>
-    static inline constexpr loc_point_type<scalar_t> to_focal_frame(
-        const loc_point_type<scalar_t> & pc_mod_point) {
-        return {};
-    }
-
-    /// Given a local polar point given in the focal frame, @returns the
-    /// correponding point in the disc frame
-    DETRAY_HOST_DEVICE
-    template<typename scalar_t>
-    static inline constexpr loc_point_type<scalar_t> to_disc_frame(
-        const loc_point_type<scalar_t> & pc_strp_point) {
-        return {};
-    }*/
-
     /// @returns the stereo angle calculated from the mask @param bounds .
     template <template <typename, std::size_t> class bounds_t,
               typename scalar_t, std::size_t kDIM,
@@ -282,6 +264,92 @@ class annulus2D {
         return corner_pos;
     }
 
+    /// Generate vertices in local cartesian frame
+    ///
+    /// @param bounds the boundary values for the stereo annulus
+    /// @param n_seg is the number of line segments
+    ///
+    /// @return a generated list of vertices
+    template <typename point2_t, typename point3_t,
+              template <typename, std::size_t> class bounds_t,
+              typename scalar_t, std::size_t kDIM,
+              typename std::enable_if_t<kDIM == e_size, bool> = true>
+    DETRAY_HOST dvector<point3_t> vertices(
+        const bounds_t<scalar_t, kDIM> &bounds, dindex n_seg) const {
+
+        scalar_t min_r = bounds[e_min_r];
+        scalar_t max_r = bounds[e_max_r];
+        scalar_t min_phi_rel = bounds[e_min_phi_rel];
+        scalar_t max_phi_rel = bounds[e_max_phi_rel];
+        scalar_t origin_x = bounds[e_shift_x];
+        scalar_t origin_y = bounds[e_shift_y];
+
+        point2_t origin_m = {origin_x, origin_y};
+
+        /// Helper method: find inner outer radius at edges in STRIP PC
+        auto circIx = [](scalar_t O_x, scalar_t O_y, scalar_t r,
+                         scalar_t phi) -> point2_t {
+            //                      _____________________________________________
+            //                     /      2  2                    2    2  2    2
+            //     O_x + O_y*m - \/  - O_x *m  + 2*O_x*O_y*m - O_y  + m *r  + r
+            // x =
+            // --------------------------------------------------------------
+            //                                  2
+            //                                 m  + 1
+            //
+            // y = m*x
+            //
+            scalar_t m = std::tan(phi);
+            point2_t dir = {math_ns::cos(phi), std::sin(phi)};
+            scalar_t x1 = (O_x + O_y * m -
+                           std::sqrt(-std::pow(O_x, 2.f) * std::pow(m, 2.f) +
+                                     2.f * O_x * O_y * m - std::pow(O_y, 2.f) +
+                                     std::pow(m, 2.f) * std::pow(r, 2.f) +
+                                     std::pow(r, 2.f))) /
+                          (std::pow(m, 2.f) + 1.f);
+            scalar_t x2 = (O_x + O_y * m +
+                           std::sqrt(-std::pow(O_x, 2.f) * std::pow(m, 2.f) +
+                                     2.f * O_x * O_y * m - std::pow(O_y, 2.f) +
+                                     std::pow(m, 2.f) * std::pow(r, 2.f) +
+                                     std::pow(r, 2.f))) /
+                          (std::pow(m, 2.f) + 1.f);
+
+            point2_t v1 = {x1, m * x1};
+            if (vector::dot(v1, dir) > 0.f)
+                return v1;
+            return {x2, m * x2};
+        };
+
+        // calculate corners in STRIP XY
+        point2_t ul_xy = circIx(origin_x, origin_y, max_r, max_phi_rel);
+        point2_t ll_xy = circIx(origin_x, origin_y, min_r, max_phi_rel);
+        point2_t ur_xy = circIx(origin_x, origin_y, max_r, min_phi_rel);
+        point2_t lr_xy = circIx(origin_x, origin_y, min_r, min_phi_rel);
+
+        auto inner_phi =
+            detail::phi_values(getter::phi(ll_xy - origin_m),
+                               getter::phi(lr_xy - origin_m), n_seg);
+        auto outer_phi =
+            detail::phi_values(getter::phi(ur_xy - origin_m),
+                               getter::phi(ul_xy - origin_m), n_seg);
+
+        dvector<point3_t> annulus_vertices;
+        annulus_vertices.reserve(inner_phi.size() + outer_phi.size());
+        for (auto iphi : inner_phi) {
+            annulus_vertices.push_back(
+                point3_t{min_r * math_ns::cos(iphi) + origin_x,
+                         min_r * std::sin(iphi) + origin_y, 0.f});
+        }
+
+        for (auto ophi : outer_phi) {
+            annulus_vertices.push_back(
+                point3_t{max_r * math_ns::cos(ophi) + origin_x,
+                         max_r * std::sin(ophi) + origin_y, 0.f});
+        }
+
+        return annulus_vertices;
+    }
+
     /// @brief Check consistency of boundary values.
     ///
     /// @param bounds the boundary values for this shape

core/include/detray/masks/cuboid3D.hpp

@@ -128,6 +128,23 @@ class cuboid3D {
         return o_bounds;
     }
 
+    /// Generate vertices in local cartesian frame
+    ///
+    /// @param bounds the boundary values for the cuboid
+    /// @param n_seg is the number of line segments
+    ///
+    /// @return a generated list of vertices
+    template <typename point2_t, typename point3_t,
+              template <typename, std::size_t> class bounds_t,
+              typename scalar_t, std::size_t kDIM,
+              typename std::enable_if_t<kDIM == e_size, bool> = true>
+    DETRAY_HOST dvector<point3_t> vertices(const bounds_t<scalar_t, kDIM> &,
+                                           dindex) const {
+        throw std::runtime_error(
+            "Vertex generation for cuboids is not implemented");
+        return {};
+    }
+
     /// @brief Check consistency of boundary values.
     ///
     /// @param bounds the boundary values for this shape

core/include/detray/masks/cylinder2D.hpp

@@ -131,6 +131,23 @@ class cylinder2D {
                 xy_bound,  xy_bound,  bounds[e_p_half_z] + env};
     }
 
+    /// Generate vertices in local cartesian frame
+    ///
+    /// @param bounds the boundary values for the cylinder
+    /// @param n_seg is the number of line segments
+    ///
+    /// @return a generated list of vertices
+    template <typename point2_t, typename point3_t,
+              template <typename, std::size_t> class bounds_t,
+              typename scalar_t, std::size_t kDIM,
+              typename std::enable_if_t<kDIM == e_size, bool> = true>
+    DETRAY_HOST dvector<point3_t> vertices(const bounds_t<scalar_t, kDIM> &,
+                                           dindex) const {
+        throw std::runtime_error(
+            "Vertex generation for cylinders is not implemented");
+        return {};
+    }
+
     /// @brief Check consistency of boundary values.
     ///
     /// @param bounds the boundary values for this shape

core/include/detray/masks/cylinder3D.hpp

@@ -124,6 +124,23 @@ class cylinder3D {
                 r_bound,  r_bound,  bounds[e_max_z] + env};
     }
 
+    /// Generate vertices in local cartesian frame
+    ///
+    /// @param bounds the boundary values for the cylinder
+    /// @param n_seg is the number of line segments
+    ///
+    /// @return a generated list of vertices
+    template <typename point2_t, typename point3_t,
+              template <typename, std::size_t> class bounds_t,
+              typename scalar_t, std::size_t kDIM,
+              typename std::enable_if_t<kDIM == e_size, bool> = true>
+    DETRAY_HOST dvector<point3_t> vertices(const bounds_t<scalar_t, kDIM> &,
+                                           dindex) const {
+        throw std::runtime_error(
+            "Vertex generation for 3D cylinders is not implemented");
+        return {};
+    }
+
     /// @brief Check consistency of boundary values.
     ///
     /// @param bounds the boundary values for this shape