Move the vertexing functionality to the masks and add support for plo…

…tting 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 lseg) const {
 
-            return m.local_min_bounds(env);
+            return mask_group[index].vertices(lseg);
         }
     };
 };

core/include/detray/geometry/surface.hpp

@@ -237,17 +237,19 @@ class surface {
                                                              pos, dir, dtds);
     }
 
-    /// @returns the vertices in local frame
+    /// @returns the vertices in local frame with @param lseg 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 lseg) const {
+        return visit_mask<typename kernels::vertices>(lseg);
     }
 
-    /// @returns the vertices in global frame.
+    /// @returns the vertices in global frame with @param lseg 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 lseg,
                                    const vector3 &dir) const {
-        auto vertices = local_vertices();
+        auto vertices = local_vertices(lseg);
         for (size_t i = 0; i < vertices.size(); i++) {
             vertices[i] = local_to_global(ctx, vertices[i], dir);
         }
@@ -269,7 +271,7 @@ class surface {
     /// @tparam functor_t the prescription to be applied to the mask
     /// @tparam Args      types of additional arguments to the functor
     template <typename functor_t, typename... Args>
-    DETRAY_HOST_DEVICE constexpr auto visit_mask(Args &&... args) const {
+    DETRAY_HOST_DEVICE constexpr auto visit_mask(Args &&...args) const {
         const auto &masks = m_detector.mask_store();
 
         return masks.template visit<functor_t>(m_desc.mask(),
@@ -281,7 +283,7 @@ class surface {
     /// @tparam functor_t the prescription to be applied to the mask
     /// @tparam Args      types of additional arguments to the functor
     template <typename functor_t, typename... Args>
-    DETRAY_HOST_DEVICE constexpr auto visit_material(Args &&... args) const {
+    DETRAY_HOST_DEVICE constexpr auto visit_material(Args &&...args) const {
         const auto &materials = m_detector.material_store();
 
         return materials.template visit<functor_t>(m_desc.material(),
@@ -370,10 +372,10 @@ class surface {
 
 template <typename detector_t, typename descr_t>
 DETRAY_HOST_DEVICE surface(const detector_t &, const descr_t &)
-    ->surface<detector_t>;
+    -> surface<detector_t>;
 
 template <typename detector_t>
 DETRAY_HOST_DEVICE surface(const detector_t &, const geometry::barcode)
-    ->surface<detector_t>;
+    -> surface<detector_t>;
 
 }  // namespace detray

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,90 @@ class annulus2D {
         return corner_pos;
     }
 
+    /// Generate vertices in local cartesian frame
+    ///
+    /// @param bounds the boundary values for the stereo annulus
+    /// @param lseg 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 lseg) 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), lseg);
+        auto outer_phi = detail::phi_values(
+            getter::phi(ur_xy - origin_m), getter::phi(ul_xy - origin_m), lseg);
+
+        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 lseg 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 lseg 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 lseg 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

core/include/detray/masks/line.hpp

@@ -151,6 +151,25 @@ class line {
         return {-xy_bound, -xy_bound, -z_bound, xy_bound, xy_bound, z_bound};
     }
 
+    /// Generate vertices in local cartesian frame
+    ///
+    /// @param bounds the boundary values for the line
+    /// @param lseg 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 /*ignored*/) const {
+
+        point3_t lc = {0.f, 0.f, -bounds[e_half_z]};
+        point3_t rc = {0.f, 0.f, bounds[e_half_z]};
+
+        return {lc, rc};
+    }
+
     /// @brief Check consistency of boundary values.
     ///
     /// @param bounds the boundary values for this shape

core/include/detray/masks/masks.hpp

@@ -220,6 +220,18 @@ class mask {
         return trf.point_to_global(center);
     }
 
+    /// @brief Vertices of the mask in local cartesian coordinates.
+    ///
+    /// Computes vertices along the mask boundary.
+    ///
+    /// @param lseg the number of segments in for arcs
+    ///
+    /// @returns a vector of vertices.
+    DETRAY_HOST
+    auto vertices(const dindex lseg) const -> dvector<point3_t> {
+        return _shape.template vertices<point2_t, point3_t>(_values, lseg);
+    }
+
     /// @returns true if the mask boundary values are consistent
     DETRAY_HOST
     constexpr bool self_check(std::ostream& os) const {