Deduplicate number_of_pixels and fix doc/naming inconsistencies

R1: three byte-identical number_of_pixels helpers (label_accumulation,
feature_accumulation, node_label_projection) — the feature_accumulation copy
was dead. Add detail::number_of_elements to detail/grid.hxx and route the two
live call sites through it; drop the dead copy.

R9: drop the redundant nb::cast wrappers inside nb::make_tuple in
relabel_sequential_t (make_tuple already casts its arguments).

Docs/naming:
- AGENTS.md: correct the union-find path (util/union_find.hxx, namespace
  bioimage_cpp::util, not detail/) and list number_of_elements under grid.hxx.
- agglomeration/detail.hxx: the rekey comment described the opposite of the
  code; reword to explain that a kRejectEdge'd edge is deliberately left out of
  the heap.
- filters: the public kwarg was window_ratio while the Python API documents
  window_size; align the nb::arg keyword and the validation message to
  window_size (calls are positional, so behavior is unchanged).
- feature_accumulation: note that float feature values are bit-reproducible
  only for a fixed thread count.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: constantinpape

AGENTS.md

@@ -24,10 +24,10 @@ Required in the environment: `scikit-build-core`, `cmake>=3.21`, `ninja`, `nanob
 
 Before introducing union-finds, priority queues, edge hashing, stride math, threading helpers, or label relabeling, check `include/bioimage_cpp/detail/`:
 
-- `detail/union_find.hxx` — `UnionFind` (path compression + union by rank). `find`, `merge`, `merge_to`, `unite_roots`.
+- `util/union_find.hxx` — `UnionFind` (namespace `bioimage_cpp::util`, not `detail`; path compression + union by rank). `find`, `merge`, `merge_to`, `unite_roots`.
 - `detail/indexed_heap.hxx` — addressable max-heap with mutable priorities. `DenseIndexedHeap` (vector-backed locator, integer keys in `[0, N)`) and `SparseIndexedHeap` (hashmap-backed, arbitrary keys).
 - `detail/edge_hash.hxx` — `Edge`, `edge_key`, `EdgeHash` for hashing unordered node pairs.
-- `detail/grid.hxx` — `c_order_strides`, `valid_offset_target`, `is_valid_grid_edge` for row-major grid offsets.
+- `detail/grid.hxx` — `number_of_elements`, `c_order_strides`, `valid_offset_target`, `is_valid_grid_edge` for row-major grid shapes/offsets.
 - `detail/relabel.hxx` — `dense_relabel` to map arbitrary labels onto `[0, k)` preserving first-occurrence order.
 - `detail/threading.hxx` — `normalize_thread_count`, `parallel_for_chunks(n_threads, n_items, chunk)`.
 

include/bioimage_cpp/detail/grid.hxx

@@ -6,6 +6,16 @@
 
 namespace bioimage_cpp::detail {
 
+// Total number of elements in a row-major array of the given shape (the product
+// of the per-axis extents). Shape entries are assumed non-negative.
+inline std::size_t number_of_elements(const std::vector<std::ptrdiff_t> &shape) {
+    std::size_t total = 1;
+    for (const auto extent : shape) {
+        total *= static_cast<std::size_t>(extent);
+    }
+    return total;
+}
+
 // C-order strides for a row-major array of the given shape, in units of array
 // elements (not bytes). The innermost (last) axis has stride 1.
 inline std::vector<std::ptrdiff_t> c_order_strides(const std::vector<std::ptrdiff_t> &shape) {

include/bioimage_cpp/graph/agglomeration/detail.hxx

@@ -93,8 +93,9 @@ inline std::size_t agglo_merge_dynamic_nodes(
             if (new_priority != current_priority) {
                 edge.weight = new_priority;
                 // The edge may have been previously popped via kRejectEdge
-                // (GASP cannot-link), in which case it is no longer in the
-                // heap. Use push_or_change to handle both cases.
+                // (GASP cannot-link). In that case it must stay out of the heap
+                // (re-adding it would let it be re-popped and merged), so only
+                // update the priority when the edge is still present.
                 if (heap.contains(removed_edge_id)) {
                     heap.change(removed_edge_id, new_priority);
                 }

include/bioimage_cpp/graph/feature_accumulation.hxx

@@ -88,15 +88,6 @@ inline double percentile(std::vector<double> &values, const double percentile) {
     return (1.0 - weight) * lower + weight * upper;
 }
 
-inline std::size_t number_of_pixels(const std::vector<std::ptrdiff_t> &shape) {
-    return static_cast<std::size_t>(std::accumulate(
-        shape.begin(),
-        shape.end(),
-        std::ptrdiff_t{1},
-        [](const std::ptrdiff_t a, const std::ptrdiff_t b) { return a * b; }
-    ));
-}
-
 inline void require_same_spatial_shape(
     const std::vector<std::ptrdiff_t> &labels_shape,
     const std::vector<std::ptrdiff_t> &data_shape,
@@ -370,6 +361,11 @@ void scan_affinity_3d_chunk(
     }
 }
 
+// Combine per-thread partial statistics. count/min/max/percentile results are
+// exact regardless of thread count; the floating-point sum / sum-of-squares
+// (and hence mean / std) can differ in the last ULP between thread counts,
+// because each thread sums a different subset of pixels. Float feature values
+// are therefore bit-reproducible only for a fixed thread count.
 template <class Stats>
 std::vector<Stats> merge_stats(
     std::vector<std::vector<Stats>> &per_thread_stats,

include/bioimage_cpp/graph/label_accumulation.hxx

@@ -1,13 +1,13 @@
 #pragma once
 
 #include "bioimage_cpp/array_view.hxx"
+#include "bioimage_cpp/detail/grid.hxx"
 #include "bioimage_cpp/detail/threading.hxx"
 #include "bioimage_cpp/graph/node_label_projection.hxx"
 #include "bioimage_cpp/graph/region_adjacency_graph.hxx"
 
 #include <cstddef>
 #include <cstdint>
-#include <numeric>
 #include <stdexcept>
 #include <string>
 #include <unordered_map>
@@ -17,15 +17,6 @@ namespace bioimage_cpp::graph {
 
 namespace detail_label_accumulation {
 
-inline std::size_t number_of_pixels(const std::vector<std::ptrdiff_t> &shape) {
-    return static_cast<std::size_t>(std::accumulate(
-        shape.begin(),
-        shape.end(),
-        std::ptrdiff_t{1},
-        [](const std::ptrdiff_t a, const std::ptrdiff_t b) { return a * b; }
-    ));
-}
-
 // Combined (node, other) histogram key. A single map per thread keyed by this
 // pair avoids the n_threads * n_nodes map-of-maps (one std::unordered_map per
 // node per thread), which dominates allocation for RAGs with many nodes.
@@ -104,7 +95,7 @@ void accumulate_labels(
         throw std::invalid_argument("labels contain a node id outside the rag");
     }
 
-    const auto n_pixels = detail_label_accumulation::number_of_pixels(labels.shape);
+    const auto n_pixels = bioimage_cpp::detail::number_of_elements(labels.shape);
     const auto n_threads = detail::normalize_thread_count(number_of_threads, n_pixels);
 
     std::vector<detail_label_accumulation::NodeOtherHistogram<OtherT>> per_thread(n_threads);

include/bioimage_cpp/graph/node_label_projection.hxx

@@ -1,12 +1,12 @@
 #pragma once
 
 #include "bioimage_cpp/array_view.hxx"
+#include "bioimage_cpp/detail/grid.hxx"
 #include "bioimage_cpp/detail/threading.hxx"
 #include "bioimage_cpp/graph/region_adjacency_graph.hxx"
 
 #include <cstddef>
 #include <cstdint>
-#include <numeric>
 #include <stdexcept>
 #include <string>
 #include <vector>
@@ -15,15 +15,6 @@ namespace bioimage_cpp::graph {
 
 namespace detail_projection {
 
-inline std::size_t number_of_pixels(const std::vector<std::ptrdiff_t> &shape) {
-    return static_cast<std::size_t>(std::accumulate(
-        shape.begin(),
-        shape.end(),
-        std::ptrdiff_t{1},
-        [](const std::ptrdiff_t a, const std::ptrdiff_t b) { return a * b; }
-    ));
-}
-
 inline void require_rag_shape_matches_labels(
     const RegionAdjacencyGraph &rag,
     const std::vector<std::ptrdiff_t> &labels_shape
@@ -69,7 +60,7 @@ void project_node_labels_to_pixels(
         throw std::invalid_argument("labels contain a node id outside node_labels");
     }
 
-    const auto n_pixels = detail_projection::number_of_pixels(labels.shape);
+    const auto n_pixels = bioimage_cpp::detail::number_of_elements(labels.shape);
     const auto n_threads = detail::normalize_thread_count(number_of_threads, n_pixels);
     bioimage_cpp::detail::parallel_for_chunks(
         n_threads,

src/bindings/filters.cxx

@@ -43,11 +43,11 @@ void require_order(int order, const char *name, const char *function) {
     }
 }
 
-void require_non_negative_window(double window_ratio, const char *function) {
-    if (window_ratio < 0.0) {
+void require_non_negative_window(double window_size, const char *function) {
+    if (window_size < 0.0) {
         throw std::invalid_argument(
-            std::string(function) + ": window_ratio must be >= 0 (0 selects the "
-            "default), got " + std::to_string(window_ratio)
+            std::string(function) + ": window_size must be >= 0 (0 selects the "
+            "default), got " + std::to_string(window_size)
         );
     }
 }
@@ -476,69 +476,69 @@ void bind_filters(nb::module_ &m) {
     m.def(
         "_gaussian_smoothing_2d_float32", &gaussian_smoothing_2d,
         nb::arg("image"), nb::arg("sigma_y"), nb::arg("sigma_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "2D Gaussian smoothing on a float32 (ny, nx) image with anisotropic sigma."
     );
     m.def(
         "_gaussian_smoothing_3d_float32", &gaussian_smoothing_3d,
         nb::arg("image"),
         nb::arg("sigma_z"), nb::arg("sigma_y"), nb::arg("sigma_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "3D Gaussian smoothing on a float32 (nz, ny, nx) image with anisotropic sigma."
     );
     m.def(
         "_gaussian_derivative_2d_float32", &gaussian_derivative_2d,
         nb::arg("image"), nb::arg("sigma_y"), nb::arg("sigma_x"),
         nb::arg("order_y"), nb::arg("order_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "2D Gaussian derivative on a float32 (ny, nx) image with per-axis order."
     );
     m.def(
         "_gaussian_derivative_3d_float32", &gaussian_derivative_3d,
         nb::arg("image"),
         nb::arg("sigma_z"), nb::arg("sigma_y"), nb::arg("sigma_x"),
         nb::arg("order_z"), nb::arg("order_y"), nb::arg("order_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "3D Gaussian derivative on a float32 (nz, ny, nx) image with per-axis order."
     );
     m.def(
         "_gaussian_gradient_magnitude_2d_float32", &gaussian_gradient_magnitude_2d,
         nb::arg("image"), nb::arg("sigma_y"), nb::arg("sigma_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "Gradient magnitude of a Gaussian-smoothed 2D float32 image."
     );
     m.def(
         "_gaussian_gradient_magnitude_3d_float32", &gaussian_gradient_magnitude_3d,
         nb::arg("image"),
         nb::arg("sigma_z"), nb::arg("sigma_y"), nb::arg("sigma_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "Gradient magnitude of a Gaussian-smoothed 3D float32 image."
     );
     m.def(
         "_laplacian_of_gaussian_2d_float32", &laplacian_of_gaussian_2d,
         nb::arg("image"), nb::arg("sigma_y"), nb::arg("sigma_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "Laplacian of Gaussian on a 2D float32 image."
     );
     m.def(
         "_laplacian_of_gaussian_3d_float32", &laplacian_of_gaussian_3d,
         nb::arg("image"),
         nb::arg("sigma_z"), nb::arg("sigma_y"), nb::arg("sigma_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "Laplacian of Gaussian on a 3D float32 image."
     );
     m.def(
         "_hessian_of_gaussian_eigenvalues_2d_float32", &hessian_of_gaussian_eigenvalues_2d,
         nb::arg("image"), nb::arg("sigma_y"), nb::arg("sigma_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "Eigenvalues of the Hessian of Gaussian on a 2D float32 image. "
         "Output shape: (ny, nx, 2), sorted descending along the trailing axis."
     );
     m.def(
         "_hessian_of_gaussian_eigenvalues_3d_float32", &hessian_of_gaussian_eigenvalues_3d,
         nb::arg("image"),
         nb::arg("sigma_z"), nb::arg("sigma_y"), nb::arg("sigma_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "Eigenvalues of the Hessian of Gaussian on a 3D float32 image. "
         "Output shape: (nz, ny, nx, 3), sorted descending along the trailing axis."
     );
@@ -547,7 +547,7 @@ void bind_filters(nb::module_ &m) {
         nb::arg("image"),
         nb::arg("sigma_inner_y"), nb::arg("sigma_inner_x"),
         nb::arg("sigma_outer_y"), nb::arg("sigma_outer_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "Eigenvalues of the structure tensor on a 2D float32 image. "
         "Output shape: (ny, nx, 2), sorted descending along the trailing axis."
     );
@@ -556,7 +556,7 @@ void bind_filters(nb::module_ &m) {
         nb::arg("image"),
         nb::arg("sigma_inner_z"), nb::arg("sigma_inner_y"), nb::arg("sigma_inner_x"),
         nb::arg("sigma_outer_z"), nb::arg("sigma_outer_y"), nb::arg("sigma_outer_x"),
-        nb::arg("window_ratio") = 0.0,
+        nb::arg("window_size") = 0.0,
         "Eigenvalues of the structure tensor on a 3D float32 image. "
         "Output shape: (nz, ny, nx, 3), sorted descending along the trailing axis."
     );

src/bindings/segmentation.cxx

@@ -392,11 +392,7 @@ nb::tuple relabel_sequential_t(
         inverse_owner
     );
 
-    return nb::make_tuple(
-        nb::cast(relabeled_array),
-        nb::cast(forward_array),
-        nb::cast(inverse_array)
-    );
+    return nb::make_tuple(relabeled_array, forward_array, inverse_array);
 }
 
 template <class InT>