Add the calculation of the centroid position to the shapes, using an …

…approximationx for the annulus shape

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

@@ -79,6 +79,16 @@ struct surface_kernels {
         }
     };
 
+    /// A functor get the mask centroid in local cartesian coordinates
+    struct centroid {
+        template <typename mask_group_t, typename index_t>
+        DETRAY_HOST_DEVICE inline point3 operator()(
+            const mask_group_t& mask_group, const index_t& index) const {
+
+            return mask_group[index].centroid();
+        }
+    };
+
     /// A functor to perform global to local bound transformation
     struct global_to_bound {
         template <typename mask_group_t, typename index_t>

core/include/detray/geometry/surface.hpp

@@ -139,7 +139,15 @@ class surface {
         return m_detector.transform_store(ctx)[m_desc.transform()];
     }
 
-    /// @returns the center position on the surface in global coordinates
+    /// @returns the centroid of the surface mask in local cartesian coordinates
+    DETRAY_HOST_DEVICE
+    constexpr auto centroid() const -> const point3 {
+        return visit_mask<typename kernels::centroid>();
+    }
+
+    /// @returns the center position of the surface in global coordinates
+    /// @note for shapes like the annulus this is not synonymous to the controid
+    /// but the focal point of the strip system instead
     DETRAY_HOST_DEVICE
     constexpr auto center(const context &ctx) const -> const point3 {
         return transform(ctx).translation();
@@ -188,16 +196,20 @@ class surface {
                                                              global, dir);
     }
 
-    /// @returns the global position to the given local/bound position @param p
+    /// @returns the global position to the given local position @param local
     /// for a given geometry context @param ctx
-    template <typename point_t,
-              std::enable_if_t<std::is_same_v<point_t, point3> or
-                                   std::is_same_v<point_t, point2>,
-                               bool> = true>
     DETRAY_HOST_DEVICE constexpr point3 local_to_global(
-        const context &ctx, const point_t &p, const vector3 &dir) const {
-        return visit_mask<typename kernels::local_to_global>(transform(ctx), p,
-                                                             dir);
+        const context &ctx, const point3 &local, const vector3 &dir) const {
+        return visit_mask<typename kernels::local_to_global>(transform(ctx),
+                                                             local, dir);
+    }
+
+    /// @returns the global position to the given bound position @param bound
+    /// for a given geometry context @param ctx
+    DETRAY_HOST_DEVICE constexpr point3 bound_to_global(
+        const context &ctx, const point2 &bound, const vector3 &dir) const {
+        return visit_mask<typename kernels::local_to_global>(transform(ctx),
+                                                             bound, dir);
     }
 
     /// @returns the track parametrization projected onto the surface (bound)
@@ -253,11 +265,11 @@ class surface {
     /// @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, const dindex n_seg,
-                                   const vector3 &dir) const {
+    constexpr auto global_vertices(const context &ctx,
+                                   const dindex n_seg) const {
         auto vertices = local_vertices(n_seg);
         for (size_t i = 0; i < vertices.size(); i++) {
-            vertices[i] = local_to_global(ctx, vertices[i], dir);
+            vertices[i] = transform(ctx).point_to_global(vertices[i]);
         }
         return vertices;
     }

core/include/detray/masks/annulus2D.hpp

@@ -264,6 +264,27 @@ class annulus2D {
         return corner_pos;
     }
 
+    /// @returns the shapes centroid in local cartesian coordinates
+    /// @note the caluculated centroid position is only an approximation
+    /// (centroid of the four corner points)!
+    template <typename algebra_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_DEVICE typename algebra_t::point3 centroid(
+        const bounds_t<scalar_t, kDIM> &bounds) const {
+
+        // Strip polar system
+        const auto crns = corners(bounds);
+
+        // Coordinates of the centroid position in strip system
+        const scalar_t r{0.25f * (crns[0] + crns[2] + crns[4] + crns[6])};
+        const scalar_t phi{bounds[e_average_phi]};
+
+        return r * typename algebra_t::point3{math::cos(phi), math_ns::sin(phi),
+                                              0.f};
+    }
+
     /// Generate vertices in local cartesian frame
     ///
     /// @param bounds the boundary values for the stereo annulus

core/include/detray/masks/cuboid3D.hpp

@@ -128,6 +128,20 @@ class cuboid3D {
         return o_bounds;
     }
 
+    /// @returns the shapes centroid in local cartesian coordinates
+    template <typename algebra_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_DEVICE typename algebra_t::point3 centroid(
+        const bounds_t<scalar_t, kDIM> &bounds) const {
+
+        return 0.5f *
+               typename algebra_t::point3{bounds[e_min_x] + bounds[e_max_x],
+                                          bounds[e_min_y] + bounds[e_max_y],
+                                          bounds[e_min_z] + bounds[e_max_z]};
+    }
+
     /// Generate vertices in local cartesian frame
     ///
     /// @param bounds the boundary values for the cuboid

core/include/detray/masks/cylinder2D.hpp

@@ -131,6 +131,17 @@ class cylinder2D {
                 xy_bound,  xy_bound,  bounds[e_p_half_z] + env};
     }
 
+    /// @returns the shapes centroid in local cartesian coordinates
+    template <typename algebra_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_DEVICE typename algebra_t::point3 centroid(
+        const bounds_t<scalar_t, kDIM> &bounds) const {
+
+        return {0.f, 0.f, 0.5f * (bounds[e_n_half_z] + bounds[e_p_half_z])};
+    }
+
     /// Generate vertices in local cartesian frame
     ///
     /// @param bounds the boundary values for the cylinder

core/include/detray/masks/cylinder3D.hpp

@@ -124,6 +124,19 @@ class cylinder3D {
                 r_bound,  r_bound,  bounds[e_max_z] + env};
     }
 
+    /// @returns the shapes centroid in local cartesian coordinates
+    template <typename algebra_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_DEVICE typename algebra_t::point3 centroid(
+        const bounds_t<scalar_t, kDIM> &bounds) const {
+
+        return 0.5f * typename algebra_t::point3{
+                          0.f, (bounds[e_min_phi] + bounds[e_max_phi]),
+                          (bounds[e_min_z] + bounds[e_max_z])};
+    }
+
     /// Generate vertices in local cartesian frame
     ///
     /// @param bounds the boundary values for the cylinder

core/include/detray/masks/line.hpp

@@ -151,6 +151,17 @@ class line {
         return {-xy_bound, -xy_bound, -z_bound, xy_bound, xy_bound, z_bound};
     }
 
+    /// @returns the shapes centroid in local cartesian coordinates
+    template <typename algebra_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_DEVICE typename algebra_t::point3 centroid(
+        const bounds_t<scalar_t, kDIM> &) const {
+
+        return {0.f, 0.f, 0.f};
+    }
+
     /// Generate vertices in local cartesian frame
     ///
     /// @param bounds the boundary values for the line

core/include/detray/masks/masks.hpp

@@ -189,6 +189,11 @@ class mask {
     DETRAY_HOST_DEVICE
     auto volume_link() -> links_type& { return _volume_link; }
 
+    /// @returns the masks centroid in local cartesian coordinates
+    DETRAY_HOST_DEVICE auto centroid() const {
+        return _shape.template centroid<algebra_t>(_values);
+    }
+
     /// @brief Lower and upper point for minimum axis aligned bounding box.
     ///
     /// Computes the min and max vertices in a local 3 dim cartesian frame.

core/include/detray/masks/rectangle2D.hpp

@@ -116,6 +116,17 @@ class rectangle2D {
         return {-x_bound, -y_bound, -env, x_bound, y_bound, env};
     }
 
+    /// @returns the shapes centroid in local cartesian coordinates
+    template <typename algebra_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_DEVICE typename algebra_t::point3 centroid(
+        const bounds_t<scalar_t, kDIM> &) const {
+
+        return {0.f, 0.f, 0.f};
+    }
+
     /// Generate vertices in local cartesian frame
     ///
     /// @param bounds the boundary values for the stereo annulus

core/include/detray/masks/ring2D.hpp

@@ -117,6 +117,17 @@ class ring2D {
         return {-r_bound, -r_bound, -env, r_bound, r_bound, env};
     }
 
+    /// @returns the shapes centroid in local cartesian coordinates
+    template <typename algebra_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_DEVICE typename algebra_t::point3 centroid(
+        const bounds_t<scalar_t, kDIM> &) const {
+
+        return {0.f, 0.f, 0.f};
+    }
+
     /// Generate vertices in local cartesian frame
     ///
     /// @param bounds the boundary values for the stereo annulus

core/include/detray/masks/single3D.hpp

@@ -117,6 +117,22 @@ class single3D {
         return o_bounds;
     }
 
+    /// @returns the shapes centroid in local cartesian coordinates
+    template <typename algebra_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_DEVICE auto centroid(
+        const bounds_t<scalar_t, kDIM> &bounds) const {
+
+        using point3_t = typename algebra_t::point3;
+
+        point3_t centr{0.f, 0.f, 0.f};
+        centr[kCheckIndex] = 0.5f * (bounds[e_lower] + bounds[e_upper]);
+
+        return centr;
+    }
+
     /// Generate vertices in local cartesian frame
     ///
     /// @param bounds the boundary values for the single value

core/include/detray/masks/trapezoid2D.hpp

@@ -130,6 +130,24 @@ class trapezoid2D {
         return {-x_bound, -y_bound, -env, x_bound, y_bound, env};
     }
 
+    /// @returns the shapes centroid in local cartesian coordinates
+    template <typename algebra_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_DEVICE typename algebra_t::point3 centroid(
+        const bounds_t<scalar_t, kDIM> &bounds) const {
+
+        const scalar_t h_2{bounds[e_half_length_2]};
+        const scalar_t a_2{bounds[e_half_length_0]};
+        const scalar_t b_2{bounds[e_half_length_1]};
+
+        const scalar_t y{2.f * h_2 * (2.f * a_2 + b_2) * 1.f /
+                         (3.f * (a_2 + b_2))};
+
+        return {0.f, y - h_2, 0.f};
+    }
+
     /// Generate vertices in local cartesian frame
     ///
     /// @param bounds the boundary values for the trapezoid

core/include/detray/masks/unbounded.hpp

@@ -79,6 +79,16 @@ class unbounded {
         return shape{}.template local_min_bounds<algebra_t>(bounds, env);
     }
 
+    /// @returns the shapes centroid in local cartesian coordinates
+    template <
+        typename algebra_t, template <typename, std::size_t> class bounds_t,
+        typename scalar_t, std::size_t kDIM,
+        typename std::enable_if_t<kDIM == boundaries::e_size, bool> = true>
+    DETRAY_HOST_DEVICE auto centroid(
+        const bounds_t<scalar_t, kDIM>& bounds) const {
+        return shape{}.template centroid<algebra_t>(bounds);
+    }
+
     /// Generate vertices in local cartesian frame
     ///
     /// @param bounds the boundary values for the underlying shape

core/include/detray/masks/unmasked.hpp

@@ -94,6 +94,16 @@ class unmasked {
         return {-inf, -inf, -inf, inf, inf, inf};
     }
 
+    /// @returns the shapes centroid in global cartesian coordinates
+    template <typename algebra_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_DEVICE typename algebra_t::point3 centroid(
+        const bounds_t<scalar_t, kDIM>&) const {
+        return {0.f, 0.f, 0.f};
+    }
+
     /// Generate vertices in local cartesian frame
     ///
     /// @param bounds the boundary values

plugins/svgtools/include/detray/plugins/svgtools/conversion/surface.hpp

@@ -146,8 +146,8 @@ auto inline surface(const transform_t& transform, const mask<annulus2D<>>& m) {
 
     auto ri = static_cast<scalar_t>(m[shape_t::e_min_r]);
     auto ro = static_cast<scalar_t>(m[shape_t::e_max_r]);
-    auto center =
-        svgtools::conversion::point<point3_t>(transform.translation());
+    auto center = svgtools::conversion::point<point3_t>(
+        transform.point_to_global(m.centroid()));
 
     p_surface._type = p_surface_t::type::e_annulus;
     p_surface._radii = {ri, ro};

plugins/svgtools/include/detray/plugins/svgtools/utils/surface_kernels.hpp

@@ -94,7 +94,7 @@ struct link_start_getter {
     template <typename mask_t, typename transform_t>
     auto inline link_start(const mask_t& mask,
                            const transform_t& transform) const {
-        return mask.global_min_bounds_center(transform);
+        return transform.point_to_global(mask.centroid());
     }
 
     // Calculates the (optimal) link starting point for rings.
@@ -134,7 +134,7 @@ struct link_start_getter {
         // Polar coordinate system.
         const typename mask_t::local_frame_type frame{};
 
-        const auto true_center = mask.global_min_bounds_center(transform);
+        const auto true_center = transform.point_to_global(mask.centroid());
         const auto rel_point =
             frame.local_to_global(transform, point3_t{r, phi, z}) -
             transform.translation();
@@ -257,4 +257,4 @@ struct link_end_getter {
     }
 };
 
-}  // namespace detray::svgtools::utils
+}  // namespace detray::svgtools::utils

tests/common/include/tests/common/tools/intersectors/helix_cylinder_intersector.hpp

@@ -47,17 +47,18 @@ struct helix_cylinder_intersector
     /// Operator function to find intersections between helix and cylinder mask
     ///
     /// @tparam mask_t is the input mask type
-    /// @tparam surface_t is the input transform type
+    /// @tparam surface_desc_t is the input transform type
     ///
     /// @param h is the input helix trajectory
+    /// @param sf_desc is the surface descriptor
     /// @param mask is the input mask
     /// @param trf is the transform
     /// @param mask_tolerance is the tolerance for mask edges
     ///
     /// @return the intersection
-    template <typename mask_t, typename surface_t>
+    template <typename mask_t, typename surface_desc_t>
     DETRAY_HOST_DEVICE inline std::array<intersection_t, 2> operator()(
-        const helix_type &h, const surface_t &sf, const mask_t &mask,
+        const helix_type &h, const surface_desc_t &sf_desc, const mask_t &mask,
         const transform3_type &trf,
         const scalar_type mask_tolerance = 0.f) const {
 
@@ -153,7 +154,7 @@ struct helix_cylinder_intersector
 
             // Compute some additional information if the intersection is valid
             if (is.status == intersection::status::e_inside) {
-                is.sf_desc = sf;
+                is.sf_desc = sf_desc;
                 is.direction = std::signbit(s)
                                    ? intersection::direction::e_opposite
                                    : intersection::direction::e_along;