diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md
index 3f6fbd8..87500a9 100644
--- a/MIGRATION_GUIDE.md
+++ b/MIGRATION_GUIDE.md
@@ -500,7 +500,11 @@ per-block subgraphs/features to zarr/N5/HDF5) is intentionally left to the
 caller, since I/O and block scheduling belong in Python.
 
 The primitives assume **globally consistent labels** (a segment has the same id
-in every block — as after a stitched distributed watershed). A block owns the
+in every block — as after a stitched distributed watershed). Labels must also
+be reasonably **dense**: the global graph allocates memory proportional to the
+largest node id (`from_unique_edges(number_of_nodes, ...)` builds a dense CSR
+over ids `0 .. number_of_nodes - 1`), so sparse or very large globally unique
+id ranges need a relabeling pass before building the graph. A block owns the
 pixel-pairs whose reference pixel lies in its inner (non-halo) box, so the
 caller reads each block with a halo (≥1 on the forward faces for the region
 graph / an edge map; ≥ `max |offset|` per side for affinities) and passes the
@@ -519,10 +523,10 @@ block_edges, block_stats = d.block_affinity_stats(labels_block, aff_block, offse
 global_edges = d.merge_edges([edges_block_0, edges_block_1, ...])
 graph = bic.graph.UndirectedGraph.from_unique_edges(number_of_nodes, global_edges)
 
-# fold per-block features onto the global edges, then finalize:
+# fold per-block features onto the global edges (in place), then finalize:
 acc = d.empty_edge_stats(graph.number_of_edges)
 for be, bs in per_block_stats:
-    acc = d.merge_block_edge_stats(graph, acc, be, bs)
+    d.merge_block_edge_stats(graph, acc, be, bs)
 features = d.finalize_edge_features(acc, compute_complex_features=True)
 ```
 
@@ -530,8 +534,11 @@ Notes:
 
 - Blocked results reproduce the whole-volume `region_adjacency_graph` /
   `features.*_features` exactly for `size`, `min` and `max`, and to
-  floating-point tolerance for `mean` / `std` (the running sums depend on thread
-  count and merge order).
+  floating-point tolerance for `mean` / `std` (the running moments depend on
+  thread count and merge order). The `(n, 5)` partial statistics use the
+  numerically stable Welford/Chan representation `[count, mean, M2, min, max]`
+  (`M2` = sum of squared deviations from the mean), so `std` stays accurate
+  for values with a large baseline and small spread.
 - **Median and percentiles are not distributable.** The distributed complex
   output is the moment subset `[mean, std, min, max, size]` — the corresponding
   columns of the in-core 12-column complex features.
diff --git a/include/bioimage_cpp/detail/grid.hxx b/include/bioimage_cpp/detail/grid.hxx
index 954d405..58686ac 100644
--- a/include/bioimage_cpp/detail/grid.hxx
+++ b/include/bioimage_cpp/detail/grid.hxx
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <algorithm>
 #include <cstddef>
 #include <cstdint>
 #include <vector>
@@ -111,6 +112,100 @@ inline void valid_axis_range(
     }
 }
 
+// Sweep every (node, target) pair on a 2D row-major grid for which
+// `node + offset` stays inside the grid and the reference node lies in the
+// half-open clip box `[clip_y_lo, clip_y_hi) x [clip_x_lo, clip_x_hi)`. The
+// body receives flat C-order indices for both endpoints and is expected to
+// inline at -O2 since this is a header-only template with a fully-known
+// callable type at instantiation. Callers fold any additional restriction —
+// a per-thread axis-0 slab, an owned block box — into the clip box; passing
+// the full grid extent sweeps every valid reference node. Shared by the
+// in-core feature accumulation / lifted-edge sweeps and the distributed
+// block extraction.
+template <class Body>
+void sweep_clipped_box_2d(
+    const std::ptrdiff_t dy,
+    const std::ptrdiff_t dx,
+    const std::size_t height,
+    const std::size_t width,
+    const std::size_t clip_y_lo,
+    const std::size_t clip_y_hi,
+    const std::size_t clip_x_lo,
+    const std::size_t clip_x_hi,
+    const Body &body
+) {
+    std::size_t y_lo_v, y_hi_v, x_lo_v, x_hi_v;
+    valid_axis_range(dy, height, y_lo_v, y_hi_v);
+    valid_axis_range(dx, width, x_lo_v, x_hi_v);
+    const auto y_lo = std::max(y_lo_v, clip_y_lo);
+    const auto y_hi = std::min(y_hi_v, clip_y_hi);
+    const auto x_lo = std::max(x_lo_v, clip_x_lo);
+    const auto x_hi = std::min(x_hi_v, clip_x_hi);
+    if (y_lo >= y_hi || x_lo >= x_hi) {
+        return;
+    }
+    const auto offset_stride = dy * static_cast<std::ptrdiff_t>(width) + dx;
+    for (std::size_t y = y_lo; y < y_hi; ++y) {
+        const auto row_offset = y * width;
+        for (std::size_t x = x_lo; x < x_hi; ++x) {
+            const auto node = row_offset + x;
+            const auto target = static_cast<std::uint64_t>(
+                static_cast<std::ptrdiff_t>(node) + offset_stride
+            );
+            body(static_cast<std::uint64_t>(node), target);
+        }
+    }
+}
+
+// 3D variant of `sweep_clipped_box_2d`.
+template <class Body>
+void sweep_clipped_box_3d(
+    const std::ptrdiff_t dz,
+    const std::ptrdiff_t dy,
+    const std::ptrdiff_t dx,
+    const std::size_t depth,
+    const std::size_t height,
+    const std::size_t width,
+    const std::size_t clip_z_lo,
+    const std::size_t clip_z_hi,
+    const std::size_t clip_y_lo,
+    const std::size_t clip_y_hi,
+    const std::size_t clip_x_lo,
+    const std::size_t clip_x_hi,
+    const Body &body
+) {
+    std::size_t z_lo_v, z_hi_v, y_lo_v, y_hi_v, x_lo_v, x_hi_v;
+    valid_axis_range(dz, depth, z_lo_v, z_hi_v);
+    valid_axis_range(dy, height, y_lo_v, y_hi_v);
+    valid_axis_range(dx, width, x_lo_v, x_hi_v);
+    const auto z_lo = std::max(z_lo_v, clip_z_lo);
+    const auto z_hi = std::min(z_hi_v, clip_z_hi);
+    const auto y_lo = std::max(y_lo_v, clip_y_lo);
+    const auto y_hi = std::min(y_hi_v, clip_y_hi);
+    const auto x_lo = std::max(x_lo_v, clip_x_lo);
+    const auto x_hi = std::min(x_hi_v, clip_x_hi);
+    if (z_lo >= z_hi || y_lo >= y_hi || x_lo >= x_hi) {
+        return;
+    }
+    const auto slice_size = height * width;
+    const auto offset_stride =
+        dz * static_cast<std::ptrdiff_t>(slice_size) +
+        dy * static_cast<std::ptrdiff_t>(width) + dx;
+    for (std::size_t z = z_lo; z < z_hi; ++z) {
+        const auto slice_offset = z * slice_size;
+        for (std::size_t y = y_lo; y < y_hi; ++y) {
+            const auto row_offset = slice_offset + y * width;
+            for (std::size_t x = x_lo; x < x_hi; ++x) {
+                const auto node = row_offset + x;
+                const auto target = static_cast<std::uint64_t>(
+                    static_cast<std::ptrdiff_t>(node) + offset_stride
+                );
+                body(static_cast<std::uint64_t>(node), target);
+            }
+        }
+    }
+}
+
 // Number of leading positions to skip along an axis of length `length` for a
 // signed offset `delta` (i.e. `max(0, -delta)` clamped to `length`). Returns a
 // plain `std::ptrdiff_t` so the affinity kernels keep their inline
diff --git a/include/bioimage_cpp/detail/threading.hxx b/include/bioimage_cpp/detail/threading.hxx
index 61028f8..3a5d22d 100644
--- a/include/bioimage_cpp/detail/threading.hxx
+++ b/include/bioimage_cpp/detail/threading.hxx
@@ -30,7 +30,7 @@ inline std::size_t normalize_thread_count(
 
 // Split [0, number_of_work_items) into n_threads contiguous chunks and invoke
 // `chunk(thread_id, begin, end)` on each chunk. Thread 0 runs on the calling
-// thread; threads 1..n_threads-1 run in parallel std::jthreads. The caller is
+// thread; threads 1..n_threads-1 run in parallel std::threads. The caller is
 // responsible for picking n_threads via normalize_thread_count and for the
 // thread safety of `chunk`.
 //
@@ -70,17 +70,27 @@ void parallel_for_chunks(
         }
     };
 
-    // jthread's destructor joins, so a failure while creating a later worker
-    // cannot destroy an earlier still-joinable thread and terminate the
-    // process. We still join explicitly below to rethrow only after all work
-    // has completed.
-    std::vector<std::jthread> threads;
+    // A joinable std::thread whose destructor runs without a prior join calls
+    // std::terminate. If spawning a later worker throws (e.g. the OS refuses a
+    // new thread), the threads vector would otherwise destroy the already
+    // joinable workers and terminate the process, so join them first and then
+    // rethrow. (std::jthread would join on destruction automatically, but it is
+    // not available in AppleClang's libc++, and portable macOS wheels are a
+    // hard requirement.)
+    std::vector<std::thread> threads;
     threads.reserve(n_threads > 0 ? n_threads - 1 : 0);
-    for (std::size_t thread_id = 1; thread_id < n_threads; ++thread_id) {
-        const auto [begin, end] = bounds(thread_id);
-        threads.emplace_back([thread_id, begin, end, &guarded_chunk]() {
-            guarded_chunk(thread_id, begin, end);
-        });
+    try {
+        for (std::size_t thread_id = 1; thread_id < n_threads; ++thread_id) {
+            const auto [begin, end] = bounds(thread_id);
+            threads.emplace_back([thread_id, begin, end, &guarded_chunk]() {
+                guarded_chunk(thread_id, begin, end);
+            });
+        }
+    } catch (...) {
+        for (auto &thread : threads) {
+            thread.join();
+        }
+        throw;
     }
     const auto [begin, end] = bounds(0);
     guarded_chunk(0, begin, end);
diff --git a/include/bioimage_cpp/graph/distributed/block_extraction.hxx b/include/bioimage_cpp/graph/distributed/block_extraction.hxx
index 4270d99..c1bd9c7 100644
--- a/include/bioimage_cpp/graph/distributed/block_extraction.hxx
+++ b/include/bioimage_cpp/graph/distributed/block_extraction.hxx
@@ -29,8 +29,8 @@
 // the block's inner (non-halo) box `[own_begin, own_begin + own_shape)`; the
 // neighbor pixel is read from the passed (outer, haloed) array. Because inner
 // boxes tile the volume, every contributing pair is counted exactly once, so
-// per-block partial statistics (count/sum/sum_of_squares add, min/max reduce)
-// reconstruct the whole-volume statistics exactly. The reference pixel is the
+// per-block partial statistics (count adds, mean/M2 combine via Chan, min/max
+// reduce) reconstruct the whole-volume statistics exactly. The reference pixel is the
 // lower pixel of a forward nearest-neighbor edge (region graph / edge map) or
 // the pixel an affinity value is stored at (affinities).
 //
@@ -42,41 +42,47 @@
 namespace bioimage_cpp::graph::distributed {
 
 // Per-edge partial statistics accumulated over the pixels a block owns. Layout
-// mirrors the serialized `(n_edges, 5)` array: count, sum, sum_of_squares, min,
-// max. `count/min/max` are exact under any thread/block ordering; the
-// floating-point `sum/sum_of_squares` (and hence mean/std after finalization)
-// are reproducible only for a fixed thread count and merge order.
+// mirrors the serialized `(n_edges, 5)` array: count, mean, M2 (sum of squared
+// deviations from the mean), min, max. Mean/M2 use the Welford update and the
+// Chan combine instead of raw sum / sum-of-squares so the merged standard
+// deviation stays accurate for values with a large baseline and small spread
+// (raw moments suffer catastrophic cancellation there). `count/min/max` are
+// exact under any thread/block ordering; the floating-point `mean/M2` (and
+// hence mean/std after finalization) are reproducible only for a fixed thread
+// count and merge order.
 struct PartialStats {
     double count = 0.0;
-    double sum = 0.0;
-    double sum_of_squares = 0.0;
+    double mean = 0.0;
+    double m2 = 0.0;
     double minimum = std::numeric_limits<double>::infinity();
     double maximum = -std::numeric_limits<double>::infinity();
 
     void add(const double value) {
-        sum += value;
-        sum_of_squares += value * value;
+        count += 1.0;
+        const auto delta = value - mean;
+        mean += delta / count;
+        m2 += delta * (value - mean);
         minimum = std::min(minimum, value);
         maximum = std::max(maximum, value);
-        count += 1.0;
     }
 
     void merge(const PartialStats &other) {
         if (other.count == 0.0) {
             return;
         }
-        sum += other.sum;
-        sum_of_squares += other.sum_of_squares;
+        const auto delta = other.mean - mean;
+        const auto total = count + other.count;
+        mean += delta * other.count / total;
+        m2 += other.m2 + delta * delta * count * other.count / total;
+        count = total;
         minimum = std::min(minimum, other.minimum);
         maximum = std::max(maximum, other.maximum);
-        count += other.count;
     }
 };
 
 // A block's owned edges together with their aligned partial statistics. `edges`
 // is sorted-unique with `u < v`; `stats` is row-major `(edges.size(), 5)` with
-// columns [count, sum, sum_of_squares, min, max] and row `i` describing
-// `edges[i]`.
+// columns [count, mean, M2, min, max] and row `i` describing `edges[i]`.
 struct BlockEdgeStats {
     std::vector<bioimage_cpp::detail::Edge> edges;
     std::vector<double> stats;
@@ -88,7 +94,6 @@ using bioimage_cpp::detail::checked_label_to_node;
 using bioimage_cpp::detail::Edge;
 using bioimage_cpp::detail::EdgeHash;
 using bioimage_cpp::detail::edge_key;
-using bioimage_cpp::detail::valid_axis_range;
 
 using EdgeSet = std::unordered_set<Edge, EdgeHash>;
 using EdgeStatsMap = std::unordered_map<Edge, PartialStats, EdgeHash>;
@@ -146,104 +151,14 @@ inline void validate_owned_box(
     }
 }
 
-// Sweep reference nodes over the owned box on a 2D grid, calling
-// `body(node, target)` with flat C-order indices into the outer array for each
-// reference node whose `+offset` neighbor stays inside the outer array. Axis 0
-// is additionally restricted to the absolute slab `[slab_begin, slab_end)`
-// (the caller's thread chunk, already inside the owned box).
-template <class Body>
-void sweep_owned_box_2d(
-    const std::ptrdiff_t dy,
-    const std::ptrdiff_t dx,
-    const std::size_t outer_h,
-    const std::size_t outer_w,
-    const std::int64_t own_begin_y,
-    const std::int64_t own_begin_x,
-    const std::int64_t own_shape_y,
-    const std::int64_t own_shape_x,
-    const std::size_t slab_begin,
-    const std::size_t slab_end,
-    const Body &body
-) {
-    std::size_t y_lo_v, y_hi_v, x_lo_v, x_hi_v;
-    valid_axis_range(dy, outer_h, y_lo_v, y_hi_v);
-    valid_axis_range(dx, outer_w, x_lo_v, x_hi_v);
-
-    const auto y_lo = std::max({y_lo_v, static_cast<std::size_t>(own_begin_y), slab_begin});
-    const auto y_hi = std::min({y_hi_v, static_cast<std::size_t>(own_begin_y + own_shape_y), slab_end});
-    const auto x_lo = std::max(x_lo_v, static_cast<std::size_t>(own_begin_x));
-    const auto x_hi = std::min(x_hi_v, static_cast<std::size_t>(own_begin_x + own_shape_x));
-    if (y_lo >= y_hi || x_lo >= x_hi) {
-        return;
-    }
-
-    const auto offset_stride = dy * static_cast<std::ptrdiff_t>(outer_w) + dx;
-    for (std::size_t y = y_lo; y < y_hi; ++y) {
-        const auto row_offset = y * outer_w;
-        for (std::size_t x = x_lo; x < x_hi; ++x) {
-            const auto node = row_offset + x;
-            const auto target = static_cast<std::uint64_t>(
-                static_cast<std::ptrdiff_t>(node) + offset_stride
-            );
-            body(static_cast<std::uint64_t>(node), target);
-        }
-    }
-}
-
-// 3D variant of `sweep_owned_box_2d`.
-template <class Body>
-void sweep_owned_box_3d(
-    const std::ptrdiff_t dz,
-    const std::ptrdiff_t dy,
-    const std::ptrdiff_t dx,
-    const std::size_t outer_d,
-    const std::size_t outer_h,
-    const std::size_t outer_w,
-    const std::int64_t own_begin_z,
-    const std::int64_t own_begin_y,
-    const std::int64_t own_begin_x,
-    const std::int64_t own_shape_z,
-    const std::int64_t own_shape_y,
-    const std::int64_t own_shape_x,
-    const std::size_t slab_begin,
-    const std::size_t slab_end,
-    const Body &body
-) {
-    std::size_t z_lo_v, z_hi_v, y_lo_v, y_hi_v, x_lo_v, x_hi_v;
-    valid_axis_range(dz, outer_d, z_lo_v, z_hi_v);
-    valid_axis_range(dy, outer_h, y_lo_v, y_hi_v);
-    valid_axis_range(dx, outer_w, x_lo_v, x_hi_v);
-
-    const auto z_lo = std::max({z_lo_v, static_cast<std::size_t>(own_begin_z), slab_begin});
-    const auto z_hi = std::min({z_hi_v, static_cast<std::size_t>(own_begin_z + own_shape_z), slab_end});
-    const auto y_lo = std::max(y_lo_v, static_cast<std::size_t>(own_begin_y));
-    const auto y_hi = std::min(y_hi_v, static_cast<std::size_t>(own_begin_y + own_shape_y));
-    const auto x_lo = std::max(x_lo_v, static_cast<std::size_t>(own_begin_x));
-    const auto x_hi = std::min(x_hi_v, static_cast<std::size_t>(own_begin_x + own_shape_x));
-    if (z_lo >= z_hi || y_lo >= y_hi || x_lo >= x_hi) {
-        return;
-    }
-
-    const auto slice_size = outer_h * outer_w;
-    const auto offset_stride =
-        dz * static_cast<std::ptrdiff_t>(slice_size) +
-        dy * static_cast<std::ptrdiff_t>(outer_w) + dx;
-    for (std::size_t z = z_lo; z < z_hi; ++z) {
-        const auto slice_offset = z * slice_size;
-        for (std::size_t y = y_lo; y < y_hi; ++y) {
-            const auto row_offset = slice_offset + y * outer_w;
-            for (std::size_t x = x_lo; x < x_hi; ++x) {
-                const auto node = row_offset + x;
-                const auto target = static_cast<std::uint64_t>(
-                    static_cast<std::ptrdiff_t>(node) + offset_stride
-                );
-                body(static_cast<std::uint64_t>(node), target);
-            }
-        }
-    }
-}
-
-// Dispatch the owned-box sweep for one offset over the correct grid rank.
+// Dispatch the owned-box sweep for one offset over the correct grid rank:
+// `body(node, target)` is called with flat C-order indices into the outer
+// array for each reference node inside the owned box whose `+offset` neighbor
+// stays inside the outer array. Axis 0 is additionally restricted to the
+// absolute slab `[slab_begin, slab_end)` (the caller's thread chunk, already
+// inside the owned box). The loop itself is the shared
+// `detail/grid.hxx::sweep_clipped_box_{2d,3d}`; the owned box and slab fold
+// into its clip box.
 template <class Body>
 void sweep_owned_box(
     const std::vector<std::ptrdiff_t> &offset,
@@ -254,18 +169,27 @@ void sweep_owned_box(
     const std::size_t slab_end,
     const Body &body
 ) {
+    const auto lo = [&](const std::size_t axis) {
+        return static_cast<std::size_t>(own_begin[axis]);
+    };
+    const auto hi = [&](const std::size_t axis) {
+        return static_cast<std::size_t>(own_begin[axis] + own_shape[axis]);
+    };
     if (outer_dims.size() == 2) {
-        sweep_owned_box_2d(
+        bioimage_cpp::detail::sweep_clipped_box_2d(
             offset[0], offset[1], outer_dims[0], outer_dims[1],
-            own_begin[0], own_begin[1], own_shape[0], own_shape[1],
-            slab_begin, slab_end, body
+            std::max(lo(0), slab_begin), std::min(hi(0), slab_end),
+            lo(1), hi(1),
+            body
         );
     } else {
-        sweep_owned_box_3d(
-            offset[0], offset[1], offset[2], outer_dims[0], outer_dims[1], outer_dims[2],
-            own_begin[0], own_begin[1], own_begin[2],
-            own_shape[0], own_shape[1], own_shape[2],
-            slab_begin, slab_end, body
+        bioimage_cpp::detail::sweep_clipped_box_3d(
+            offset[0], offset[1], offset[2],
+            outer_dims[0], outer_dims[1], outer_dims[2],
+            std::max(lo(0), slab_begin), std::min(hi(0), slab_end),
+            lo(1), hi(1),
+            lo(2), hi(2),
+            body
         );
     }
 }
@@ -357,8 +281,8 @@ inline BlockEdgeStats merge_stats_maps(std::vector<EdgeStatsMap> &per_thread) {
     for (const auto &edge : result.edges) {
         const auto &stats = merged[edge];
         result.stats.push_back(stats.count);
-        result.stats.push_back(stats.sum);
-        result.stats.push_back(stats.sum_of_squares);
+        result.stats.push_back(stats.mean);
+        result.stats.push_back(stats.m2);
         result.stats.push_back(stats.minimum);
         result.stats.push_back(stats.maximum);
     }
diff --git a/include/bioimage_cpp/graph/distributed/merge.hxx b/include/bioimage_cpp/graph/distributed/merge.hxx
index 040d430..9ea05f0 100644
--- a/include/bioimage_cpp/graph/distributed/merge.hxx
+++ b/include/bioimage_cpp/graph/distributed/merge.hxx
@@ -9,7 +9,6 @@
 #include <cstddef>
 #include <cstdint>
 #include <stdexcept>
-#include <unordered_set>
 #include <vector>
 
 // Whole-volume merge primitives for the distributed region-adjacency-graph and
@@ -25,7 +24,9 @@ namespace bioimage_cpp::graph::distributed {
 // Edges are canonicalized to `u < v`, self-edges (`u == v`) are dropped, and the
 // result is deduplicated and sorted ascending by `(u, v)` — exactly the
 // precondition of `UndirectedGraph::from_sorted_unique_edges`, so the merged
-// output can be turned into the global graph directly.
+// output can be turned into the global graph directly. Deduplication uses
+// sort + unique rather than a hash set: peak memory stays close to one copy of
+// the input instead of the ~8x node/bucket overhead of an unordered_set.
 inline std::vector<bioimage_cpp::detail::Edge> merge_edges(
     const ConstArrayView<std::uint64_t> &concatenated
 ) {
@@ -34,36 +35,41 @@ inline std::vector<bioimage_cpp::detail::Edge> merge_edges(
     }
 
     const auto number_of_input_edges = static_cast<std::size_t>(concatenated.shape[0]);
-    std::unordered_set<bioimage_cpp::detail::Edge, bioimage_cpp::detail::EdgeHash> merged;
-    merged.reserve(number_of_input_edges);
+    std::vector<bioimage_cpp::detail::Edge> sorted_edges;
+    sorted_edges.reserve(number_of_input_edges);
     for (std::size_t index = 0; index < number_of_input_edges; ++index) {
         const auto u = concatenated.data[2 * index];
         const auto v = concatenated.data[2 * index + 1];
         if (u == v) {
             continue;
         }
-        merged.insert(bioimage_cpp::detail::edge_key(u, v));
+        sorted_edges.push_back(bioimage_cpp::detail::edge_key(u, v));
     }
 
-    std::vector<bioimage_cpp::detail::Edge> sorted_edges(merged.begin(), merged.end());
     std::sort(sorted_edges.begin(), sorted_edges.end());
+    sorted_edges.erase(
+        std::unique(sorted_edges.begin(), sorted_edges.end()), sorted_edges.end()
+    );
     return sorted_edges;
 }
 
 // Fold one block's partial edge statistics into a running whole-volume
-// accumulator, returning the updated accumulator. `current_stats` and the return
-// value are `(number_of_edges, 5)` rows aligned to `global_graph`'s edge ids;
+// accumulator, updating `current_stats` in place. `current_stats` is
+// `(number_of_edges, 5)` rows aligned to `global_graph`'s edge ids;
 // `block_edges`/`block_stats` are the `(n, 2)` / `(n, 5)` outputs of a block
 // extraction. Each block edge is mapped to its global edge id via
 // `global_graph.find_edge`; edges absent from the global graph (id `-1`) are
-// skipped. `count/sum/sum_of_squares` add; `min/max` reduce, seeded from the
-// first contribution to each edge (so `current_stats` may be zero-initialized).
+// skipped. `count` adds, `mean/M2` combine via the Chan formula, and `min/max`
+// reduce; rows with zero count take the block row verbatim (so `current_stats`
+// may be zero-initialized). Mutating in place keeps one merge O(block edges)
+// instead of O(global edges) — the per-block cost must not scale with the
+// global graph.
 //
 // Block edge endpoints must be valid node ids of `global_graph`
 // (`< number_of_nodes`); `find_edge` throws `std::out_of_range` otherwise.
-inline std::vector<double> merge_block_edge_stats(
+inline void merge_block_edge_stats(
     const UndirectedGraph &global_graph,
-    const ConstArrayView<double> &current_stats,
+    const ArrayView<double> &current_stats,
     const ConstArrayView<std::uint64_t> &block_edges,
     const ConstArrayView<double> &block_stats
 ) {
@@ -87,11 +93,6 @@ inline std::vector<double> merge_block_edge_stats(
         );
     }
 
-    const auto number_of_edges = static_cast<std::size_t>(current_stats.shape[0]);
-    std::vector<double> out(
-        current_stats.data, current_stats.data + number_of_edges * 5
-    );
-
     const auto number_of_block_edges = static_cast<std::size_t>(block_edges.shape[0]);
     for (std::size_t index = 0; index < number_of_block_edges; ++index) {
         const auto u = block_edges.data[2 * index];
@@ -101,20 +102,27 @@ inline std::vector<double> merge_block_edge_stats(
             continue;
         }
 
-        double *const o = out.data() + static_cast<std::size_t>(edge) * 5;
+        double *const o = current_stats.data + static_cast<std::size_t>(edge) * 5;
         const double *const b = block_stats.data + index * 5;
+        if (b[0] == 0.0) {
+            continue;
+        }
         if (o[0] == 0.0) {
+            o[0] = b[0];
+            o[1] = b[1];
+            o[2] = b[2];
             o[3] = b[3];
             o[4] = b[4];
         } else {
+            const auto delta = b[1] - o[1];
+            const auto total = o[0] + b[0];
+            o[1] += delta * b[0] / total;
+            o[2] += b[2] + delta * delta * o[0] * b[0] / total;
+            o[0] = total;
             o[3] = std::min(o[3], b[3]);
             o[4] = std::max(o[4], b[4]);
         }
-        o[0] += b[0];
-        o[1] += b[1];
-        o[2] += b[2];
     }
-    return out;
 }
 
 // Turn accumulated partial statistics `(number_of_edges, 5)` into edge features.
@@ -152,9 +160,9 @@ inline void finalize_edge_features(
             continue;
         }
 
-        const auto mean = s[1] / count;
+        const auto mean = s[1];
         if (compute_complex_features) {
-            const auto variance = std::max(0.0, s[2] / count - mean * mean);
+            const auto variance = std::max(0.0, s[2] / count);
             o[0] = mean;
             o[1] = std::sqrt(variance);
             o[2] = s[3];
diff --git a/include/bioimage_cpp/graph/feature_accumulation.hxx b/include/bioimage_cpp/graph/feature_accumulation.hxx
index 1e13ef2..f05a698 100644
--- a/include/bioimage_cpp/graph/feature_accumulation.hxx
+++ b/include/bioimage_cpp/graph/feature_accumulation.hxx
@@ -38,7 +38,6 @@ struct SimpleStats {
 
 struct ComplexStats {
     double sum = 0.0;
-    double sum_of_squares = 0.0;
     double minimum = std::numeric_limits<double>::infinity();
     double maximum = -std::numeric_limits<double>::infinity();
     std::uint64_t count = 0;
@@ -46,7 +45,6 @@ struct ComplexStats {
 
     void add(const double value) {
         sum += value;
-        sum_of_squares += value * value;
         minimum = std::min(minimum, value);
         maximum = std::max(maximum, value);
         ++count;
@@ -55,7 +53,6 @@ struct ComplexStats {
 
     void merge(ComplexStats &other) {
         sum += other.sum;
-        sum_of_squares += other.sum_of_squares;
         minimum = std::min(minimum, other.minimum);
         maximum = std::max(maximum, other.maximum);
         count += other.count;
@@ -196,16 +193,10 @@ void scan_edge_map_3d_chunk(
     }
 }
 
-// `valid_axis_range` lives in `detail/grid.hxx` (it is pure grid geometry,
-// shared with the distributed block-extraction sweeps). Pull it into this
-// namespace so the offset-box sweeps below can call it unqualified.
-using bioimage_cpp::detail::valid_axis_range;
-
 // Sweep every (node, target) pair on a 2D grid for which `node + offset` stays
-// in bounds, restricted to the half-open y-slab [y_begin, y_end). The body
-// receives flat C-order indices for both endpoints and is expected to inline
-// at -O2 since this is a header-only template with a fully-known callable
-// type at instantiation.
+// in bounds, restricted to the half-open y-slab [y_begin, y_end). Thin wrapper
+// around `detail/grid.hxx::sweep_clipped_box_2d` (the loop shared with the
+// distributed block extraction) with the x-axis unclipped.
 template <class Body>
 void sweep_offset_box_2d(
     const std::ptrdiff_t dy,
@@ -216,25 +207,9 @@ void sweep_offset_box_2d(
     const std::size_t y_end,
     const Body &body
 ) {
-    std::size_t y_lo_full, y_hi_full, x_lo, x_hi;
-    valid_axis_range(dy, height, y_lo_full, y_hi_full);
-    valid_axis_range(dx, width, x_lo, x_hi);
-    const auto y_lo = std::max(y_lo_full, y_begin);
-    const auto y_hi = std::min(y_hi_full, y_end);
-    if (y_lo >= y_hi || x_lo >= x_hi) {
-        return;
-    }
-    const auto offset_stride = dy * static_cast<std::ptrdiff_t>(width) + dx;
-    for (std::size_t y = y_lo; y < y_hi; ++y) {
-        const auto row_offset = y * width;
-        for (std::size_t x = x_lo; x < x_hi; ++x) {
-            const auto node = row_offset + x;
-            const auto target = static_cast<std::uint64_t>(
-                static_cast<std::ptrdiff_t>(node) + offset_stride
-            );
-            body(static_cast<std::uint64_t>(node), target);
-        }
-    }
+    bioimage_cpp::detail::sweep_clipped_box_2d(
+        dy, dx, height, width, y_begin, y_end, 0, width, body
+    );
 }
 
 // 3D variant of `sweep_offset_box_2d`. Restricts the sweep to a z-slab.
@@ -250,32 +225,9 @@ void sweep_offset_box_3d(
     const std::size_t z_end,
     const Body &body
 ) {
-    std::size_t z_lo_full, z_hi_full, y_lo, y_hi, x_lo, x_hi;
-    valid_axis_range(dz, depth, z_lo_full, z_hi_full);
-    valid_axis_range(dy, height, y_lo, y_hi);
-    valid_axis_range(dx, width, x_lo, x_hi);
-    const auto z_lo = std::max(z_lo_full, z_begin);
-    const auto z_hi = std::min(z_hi_full, z_end);
-    if (z_lo >= z_hi || y_lo >= y_hi || x_lo >= x_hi) {
-        return;
-    }
-    const auto slice_size = height * width;
-    const auto offset_stride =
-        dz * static_cast<std::ptrdiff_t>(slice_size) +
-        dy * static_cast<std::ptrdiff_t>(width) + dx;
-    for (std::size_t z = z_lo; z < z_hi; ++z) {
-        const auto slice_offset = z * slice_size;
-        for (std::size_t y = y_lo; y < y_hi; ++y) {
-            const auto row_offset = slice_offset + y * width;
-            for (std::size_t x = x_lo; x < x_hi; ++x) {
-                const auto node = row_offset + x;
-                const auto target = static_cast<std::uint64_t>(
-                    static_cast<std::ptrdiff_t>(node) + offset_stride
-                );
-                body(static_cast<std::uint64_t>(node), target);
-            }
-        }
-    }
+    bioimage_cpp::detail::sweep_clipped_box_3d(
+        dz, dy, dx, depth, height, width, z_begin, z_end, 0, height, 0, width, body
+    );
 }
 
 template <class LabelT, class ValueT, class Stats>
@@ -407,7 +359,16 @@ inline void write_complex_features(
 
         const auto count = static_cast<double>(edge_stats.count);
         const auto mean = edge_stats.sum / count;
-        const auto variance = std::max(0.0, edge_stats.sum_of_squares / count - mean * mean);
+        // Two-pass variance over the stored values: exact, and free of the
+        // catastrophic cancellation that sum_of_squares/count - mean^2 shows
+        // for values with a large baseline and small spread. Computed before
+        // the percentile calls, which reorder `values` via nth_element.
+        double sum_of_squared_deviations = 0.0;
+        for (const auto value : edge_stats.values) {
+            const auto deviation = value - mean;
+            sum_of_squared_deviations += deviation * deviation;
+        }
+        const auto variance = sum_of_squared_deviations / count;
         out.data[offset] = mean;
         out.data[offset + 1] = percentile(edge_stats.values, 50.0);
         out.data[offset + 2] = std::sqrt(variance);
diff --git a/src/bindings/graph.cxx b/src/bindings/graph.cxx
index 26d929b..bbc7139 100644
--- a/src/bindings/graph.cxx
+++ b/src/bindings/graph.cxx
@@ -1441,6 +1441,14 @@ std::vector<std::ptrdiff_t> ndarray_shape(ConstDoubleArray array) {
     return shape;
 }
 
+std::vector<std::ptrdiff_t> ndarray_shape(DoubleArray array) {
+    std::vector<std::ptrdiff_t> shape(array.ndim());
+    for (std::size_t axis = 0; axis < array.ndim(); ++axis) {
+        shape[axis] = static_cast<std::ptrdiff_t>(array.shape(axis));
+    }
+    return shape;
+}
+
 // ---- Distributed region-adjacency-graph + edge-feature primitives ----
 
 template <class T>
@@ -1520,24 +1528,25 @@ UInt64Array distributed_merge_edges(ConstUInt64Array edges) {
     return edges_to_uv_array(merged);
 }
 
-DoubleArray distributed_merge_block_edge_stats(
+// Mutates `current_stats` in place (the Python wrapper hands the caller's
+// accumulator back), so one merge stays O(block edges) instead of copying the
+// whole global accumulator per block.
+void distributed_merge_block_edge_stats(
     const Graph &global_graph,
-    ConstDoubleArray current_stats,
+    DoubleArray current_stats,
     ConstUInt64Array block_edges,
     ConstDoubleArray block_stats
 ) {
-    ConstArrayView<double> current_view{current_stats.data(), ndarray_shape(current_stats), {}};
+    ArrayView<double> current_view{current_stats.data(), ndarray_shape(current_stats), {}};
     ConstArrayView<std::uint64_t> block_edges_view{block_edges.data(), ndarray_shape(block_edges), {}};
     ConstArrayView<double> block_stats_view{block_stats.data(), ndarray_shape(block_stats), {}};
 
-    std::vector<double> out;
     {
         nb::gil_scoped_release release;
-        out = graph::distributed::merge_block_edge_stats(
+        graph::distributed::merge_block_edge_stats(
             global_graph, current_view, block_edges_view, block_stats_view
         );
     }
-    return block_stats_to_array(out);
 }
 
 DoubleArray distributed_finalize_edge_features(
@@ -1901,6 +1910,23 @@ void bind_graph(nb::module_ &m) {
             nb::arg("number_of_nodes") = 0,
             nb::arg("reserve_number_of_edges") = 0
         )
+        // Constructor equivalent of the from_edges / from_unique_edges
+        // statics. Unlike a def_static, an __init__ overload constructs an
+        // instance of the *derived* Python class, so the Python subclass in
+        // bioimage_cpp.graph can build itself from an edge array.
+        .def(
+            "__init__",
+            [](Graph *self, const std::uint64_t number_of_nodes, ConstUInt64Array uvs,
+               const bool unique) {
+                new (self) Graph(
+                    unique ? graph_from_unique_edges(number_of_nodes, uvs)
+                           : graph_from_edges(number_of_nodes, uvs)
+                );
+            },
+            nb::arg("number_of_nodes"),
+            nb::arg("uvs"),
+            nb::arg("unique")
+        )
         .def(
             "assign",
             &Graph::assign,
diff --git a/src/bioimage_cpp/graph/__init__.py b/src/bioimage_cpp/graph/__init__.py
index 3b4fe0a..e6fc841 100644
--- a/src/bioimage_cpp/graph/__init__.py
+++ b/src/bioimage_cpp/graph/__init__.py
@@ -97,9 +97,19 @@ def edgesFromNodeList(self, nodes):
 
     @classmethod
     def from_edges(cls, number_of_nodes: int, uvs):
-        graph = cls(number_of_nodes)
-        graph.insert_edges(uvs)
-        return graph
+        return cls(int(number_of_nodes), _as_uv_array(uvs, "uvs"), False)
+
+    @classmethod
+    def from_unique_edges(cls, number_of_nodes: int, uvs):
+        """Bulk-construct from a pre-deduplicated ``(n, 2)`` edge array.
+
+        The caller asserts that no undirected ``(u, v)`` pair appears twice and
+        that ``u != v`` in every row; edges receive ids matching their position
+        in ``uvs``. This skips the per-edge dedup of :meth:`from_edges` and is
+        the fast path for e.g. the merged edge set of
+        :func:`bioimage_cpp.graph.distributed.merge_edges`.
+        """
+        return cls(int(number_of_nodes), _as_uv_array(uvs, "uvs"), True)
 
     @classmethod
     def deserialize(cls, serialization):
diff --git a/src/bioimage_cpp/graph/distributed.py b/src/bioimage_cpp/graph/distributed.py
index 583f1dd..90893a2 100644
--- a/src/bioimage_cpp/graph/distributed.py
+++ b/src/bioimage_cpp/graph/distributed.py
@@ -21,12 +21,15 @@
 
 2. **Merge the graph**: :func:`merge_edges` unions the per-block edges into the
    whole-volume edge set; build the global graph with
-   :meth:`bioimage_cpp.graph.UndirectedGraph.from_unique_edges`.
+   :meth:`bioimage_cpp.graph.UndirectedGraph.from_unique_edges`. Labels must be
+   globally consistent *and* reasonably dense: the graph allocates memory
+   proportional to the largest node id, so sparse or very large globally unique
+   id ranges need a relabeling pass first.
 
 3. **Merge the features**: fold each block's partial statistics onto the global
    edges with :func:`merge_block_edge_stats` (starting from
-   :func:`empty_edge_stats`), then convert to features with
-   :func:`finalize_edge_features`.
+   :func:`empty_edge_stats`; the accumulator is updated in place), then convert
+   to features with :func:`finalize_edge_features`.
 
 Exactly-recoverable features are ``size``, ``mean``, ``std``, ``min`` and
 ``max``. Median and percentiles cannot be reconstructed from block partials, so
@@ -36,6 +39,8 @@
 
 from __future__ import annotations
 
+import operator
+
 import numpy as np
 
 from .. import _core
@@ -68,6 +73,8 @@ def _as_index_vector(values, ndim: int, name: str) -> list[int]:
     array = np.asarray(values)
     if array.ndim != 1 or array.shape[0] != ndim:
         raise ValueError(f"{name} must be a 1D sequence of length {ndim}")
+    if not np.issubdtype(array.dtype, np.integer):
+        raise TypeError(f"{name} must contain integers, got dtype={array.dtype}")
     return [int(value) for value in array]
 
 
@@ -78,6 +85,18 @@ def _as_stats_array(values, name: str) -> np.ndarray:
     return np.ascontiguousarray(array)
 
 
+def _as_edge_part(part) -> np.ndarray:
+    array = np.asarray(part)
+    if array.size:
+        if not np.issubdtype(array.dtype, np.integer):
+            raise TypeError(
+                f"edges must have an integer dtype, got dtype={array.dtype}"
+            )
+        if np.issubdtype(array.dtype, np.signedinteger) and (array < 0).any():
+            raise ValueError("edges must not contain negative node ids")
+    return array.astype(np.uint64, copy=False)
+
+
 def _dispatch_labels(labels, table, name: str):
     label_array = _normalize_labels(labels)
     try:
@@ -131,8 +150,9 @@ def block_edge_map_stats(
     must have the same shape as ``labels``. Returns ``(edges, stats)`` where
     ``edges`` is ``(n, 2)`` ``uint64`` (matching
     :func:`block_region_adjacency_edges` for the same block) and ``stats`` is
-    ``(n, 5)`` ``float64`` with columns ``[count, sum, sum_of_squares, min,
-    max]`` aligned row-by-row to ``edges``.
+    ``(n, 5)`` ``float64`` with columns ``[count, mean, M2, min, max]`` aligned
+    row-by-row to ``edges`` (``M2`` is the sum of squared deviations from the
+    mean, as in Welford's algorithm — numerically stable when merged).
     """
     label_array, run = _dispatch_labels(
         labels, _BLOCK_EDGE_MAP_STATS_BY_DTYPE, "edge-map"
@@ -184,7 +204,14 @@ def block_affinity_stats(
             f"affinities shape={affinity_array.shape}, labels shape={label_array.shape}"
         )
 
-    normalized_offsets = [tuple(int(value) for value in offset) for offset in offsets]
+    try:
+        # operator.index accepts ints and NumPy integers but rejects floats,
+        # which int() would silently truncate.
+        normalized_offsets = [
+            tuple(operator.index(value) for value in offset) for offset in offsets
+        ]
+    except TypeError as error:
+        raise TypeError(f"offsets must contain integers: {error}") from error
     if len(normalized_offsets) != affinity_array.shape[0]:
         raise ValueError(
             "offsets length must match affinities channel count, got "
@@ -213,16 +240,19 @@ def merge_edges(edges) -> np.ndarray:
     self-edges are dropped, and the result is deduplicated and sorted
     lexicographically. The returned ``(E, 2)`` ``uint64`` array feeds
     :meth:`bioimage_cpp.graph.UndirectedGraph.from_unique_edges` directly.
+
+    Outputs are valid inputs, so peak memory can be bounded by merging
+    hierarchically: merge groups of blocks first, then merge the merged
+    results (tree reduction).
     """
     if isinstance(edges, np.ndarray):
-        stacked = edges
+        stacked = _as_edge_part(edges)
     else:
-        parts = [np.asarray(part, dtype=np.uint64) for part in edges]
+        parts = [_as_edge_part(part) for part in edges]
         if not parts:
             return np.empty((0, 2), dtype=np.uint64)
         stacked = np.concatenate(parts, axis=0)
 
-    stacked = np.asarray(stacked, dtype=np.uint64)
     if stacked.ndim != 2 or stacked.shape[1] != 2:
         raise ValueError("edges must have shape (n_edges, 2)")
     return _core._merge_edges(np.ascontiguousarray(stacked))
@@ -245,17 +275,33 @@ def merge_block_edge_stats(
     """Fold one block's partial statistics onto the global edges.
 
     ``current_stats`` is the running ``(E, 5)`` accumulator (rows aligned to
-    ``global_graph`` edge ids; start from :func:`empty_edge_stats`).
+    ``global_graph`` edge ids; start from :func:`empty_edge_stats`). **It is
+    updated in place and returned**, so one merge costs O(block edges)
+    regardless of the global graph size; it must be a C-contiguous, writable
+    ``float64`` array (as produced by :func:`empty_edge_stats`).
     ``block_edges`` / ``block_stats`` are a block extraction's ``(n, 2)`` /
     ``(n, 5)`` outputs. Each block edge is mapped to its global edge id via
     ``global_graph.find_edge``; edges absent from the graph are skipped.
-    ``count/sum/sum_of_squares`` add and ``min/max`` reduce. Returns the updated
-    ``(E, 5)`` accumulator.
+    ``count`` adds, ``mean/M2`` combine via the Chan formula, and ``min/max``
+    reduce.
 
     Block edge endpoints must be valid node ids of ``global_graph``.
     """
-    current = _as_stats_array(current_stats, "current_stats")
-    if int(current.shape[0]) != int(global_graph.number_of_edges):
+    # Do not coerce/copy `current_stats` — the merge mutates it in place, and a
+    # silent copy would discard the update.
+    if not isinstance(current_stats, np.ndarray) or current_stats.dtype != np.float64:
+        raise TypeError(
+            "current_stats must be a float64 numpy array, got "
+            f"{type(current_stats).__name__}"
+            + (f" with dtype={current_stats.dtype}" if isinstance(current_stats, np.ndarray) else "")
+        )
+    if current_stats.ndim != 2 or current_stats.shape[1] != 5:
+        raise ValueError("current_stats must have shape (number_of_edges, 5)")
+    if not current_stats.flags.c_contiguous:
+        raise ValueError("current_stats must be C-contiguous (it is updated in place)")
+    if not current_stats.flags.writeable:
+        raise ValueError("current_stats must be writable (it is updated in place)")
+    if int(current_stats.shape[0]) != int(global_graph.number_of_edges):
         raise ValueError("current_stats rows must match global_graph number_of_edges")
 
     block_edge_array = np.asarray(block_edges, dtype=np.uint64)
@@ -265,12 +311,13 @@ def merge_block_edge_stats(
     if block_stat_array.shape[0] != block_edge_array.shape[0]:
         raise ValueError("block_edges and block_stats must have the same number of rows")
 
-    return _core._merge_block_edge_stats(
+    _core._merge_block_edge_stats(
         global_graph,
-        current,
+        current_stats,
         np.ascontiguousarray(block_edge_array),
         block_stat_array,
     )
+    return current_stats
 
 
 def finalize_edge_features(
diff --git a/src/bioimage_cpp/label_multiset/__init__.py b/src/bioimage_cpp/label_multiset/__init__.py
index 807f6a6..dc5b96a 100644
--- a/src/bioimage_cpp/label_multiset/__init__.py
+++ b/src/bioimage_cpp/label_multiset/__init__.py
@@ -22,7 +22,7 @@
 from __future__ import annotations
 
 from dataclasses import dataclass
-from typing import Optional, Tuple, Union
+from typing import Tuple
 
 import numpy as np
 
diff --git a/tests/graph/test_distributed_rag.py b/tests/graph/test_distributed_rag.py
index 9ce2cb4..b2f6238 100644
--- a/tests/graph/test_distributed_rag.py
+++ b/tests/graph/test_distributed_rag.py
@@ -284,29 +284,34 @@ def test_merge_block_edge_stats_reduction_and_skipping():
     )
     acc = dist.empty_edge_stats(graph.number_of_edges)
 
-    # Block A touches edge (0,1) with values summarized as count=2,sum=3,sq=5,min=1,max=2.
+    # Block A touches edge (0,1) with values {1, 2}: count=2, mean=1.5, M2=0.5.
     edges_a = np.array([[0, 1]], dtype=np.uint64)
-    stats_a = np.array([[2.0, 3.0, 5.0, 1.0, 2.0]], dtype=np.float64)
+    stats_a = np.array([[2.0, 1.5, 0.5, 1.0, 2.0]], dtype=np.float64)
     acc = dist.merge_block_edge_stats(graph, acc, edges_a, stats_a)
 
-    # Block B touches edge (0,1) again plus (0,2) which is absent -> skipped.
+    # Block B touches edge (0,1) again with value {10} plus (0,2) which is
+    # absent from the graph -> skipped.
     edges_b = np.array([[0, 1], [0, 2]], dtype=np.uint64)
     stats_b = np.array(
-        [[1.0, 10.0, 100.0, 10.0, 10.0], [5.0, 5.0, 5.0, -1.0, 9.0]], dtype=np.float64
+        [[1.0, 10.0, 0.0, 10.0, 10.0], [5.0, 1.0, 5.0, -1.0, 9.0]], dtype=np.float64
     )
     acc = dist.merge_block_edge_stats(graph, acc, edges_b, stats_b)
 
-    # edge (0,1): count 2+1, sum 3+10, sq 5+100, min(1,10)=1, max(2,10)=10
-    np.testing.assert_array_equal(acc[0], [3.0, 13.0, 105.0, 1.0, 10.0])
+    # edge (0,1) now describes {1, 2, 10}: count=3, mean=13/3,
+    # M2 = sum((x - mean)^2) = 438/9, min(1,10)=1, max(2,10)=10
+    values = np.array([1.0, 2.0, 10.0])
+    np.testing.assert_allclose(
+        acc[0], [3.0, values.mean(), ((values - values.mean()) ** 2).sum(), 1.0, 10.0]
+    )
     # edge (1,2): never touched
     np.testing.assert_array_equal(acc[1], [0.0, 0.0, 0.0, 0.0, 0.0])
 
 
 def test_finalize_edge_features_formulas_and_zero_count():
-    # edge 0: count=4, sum=8 -> mean 2; sum_sq=24 -> var = 24/4 - 4 = 2 -> std sqrt(2)
+    # edge 0: count=4, mean=2, M2=8 -> var = 8/4 = 2 -> std sqrt(2)
     # edge 1: empty
     stats = np.array(
-        [[4.0, 8.0, 24.0, -1.0, 5.0], [0.0, 0.0, 0.0, 0.0, 0.0]], dtype=np.float64
+        [[4.0, 2.0, 8.0, -1.0, 5.0], [0.0, 0.0, 0.0, 0.0, 0.0]], dtype=np.float64
     )
     simple = dist.finalize_edge_features(stats, compute_complex_features=False)
     np.testing.assert_allclose(simple[0], [2.0, 4.0])
@@ -317,6 +322,67 @@ def test_finalize_edge_features_formulas_and_zero_count():
     np.testing.assert_array_equal(complex_[1], [0.0, 0.0, 0.0, 0.0, 0.0])
 
 
+def test_merge_block_edge_stats_updates_in_place():
+    graph = bic.graph.UndirectedGraph.from_unique_edges(
+        2, np.array([[0, 1]], dtype=np.uint64)
+    )
+    acc = dist.empty_edge_stats(graph.number_of_edges)
+    edges = np.array([[0, 1]], dtype=np.uint64)
+    stats = np.array([[1.0, 2.0, 0.0, 2.0, 2.0]], dtype=np.float64)
+
+    result = dist.merge_block_edge_stats(graph, acc, edges, stats)
+    assert result is acc
+    np.testing.assert_array_equal(acc[0], [1.0, 2.0, 0.0, 2.0, 2.0])
+
+
+def test_merge_block_edge_stats_rejects_bad_accumulator():
+    graph = bic.graph.UndirectedGraph.from_unique_edges(
+        3, np.array([[0, 1], [1, 2]], dtype=np.uint64)
+    )
+    edges = np.array([[0, 1]], dtype=np.uint64)
+    stats = np.array([[1.0, 2.0, 0.0, 2.0, 2.0]], dtype=np.float64)
+
+    with pytest.raises(TypeError, match="float64"):
+        dist.merge_block_edge_stats(
+            graph, np.zeros((2, 5), dtype=np.float32), edges, stats
+        )
+    with pytest.raises(TypeError, match="float64"):
+        dist.merge_block_edge_stats(graph, [[0.0] * 5] * 2, edges, stats)
+    with pytest.raises(ValueError, match="C-contiguous"):
+        dist.merge_block_edge_stats(
+            graph, np.zeros((2, 10), dtype=np.float64)[:, ::2], edges, stats
+        )
+    read_only = dist.empty_edge_stats(graph.number_of_edges)
+    read_only.setflags(write=False)
+    with pytest.raises(ValueError, match="writable"):
+        dist.merge_block_edge_stats(graph, read_only, edges, stats)
+
+
+def test_std_is_stable_for_large_baseline():
+    # Values 1e8 and 1e8 + 1 have std 0.5; the naive sum-of-squares formula
+    # returns 0.0 here due to catastrophic cancellation.
+    labels = np.array([[0, 1], [0, 1]], dtype=np.uint32)
+    edge_map = np.array([[1e8, 1e8], [1e8 + 1, 1e8 + 1]], dtype=np.float64)
+
+    graph = bic.graph.UndirectedGraph.from_unique_edges(
+        2, np.array([[0, 1]], dtype=np.uint64)
+    )
+
+    # Whole array as one block, and split into two one-row blocks so the
+    # cross-block Chan combine is exercised too.
+    single = dist.block_edge_map_stats(labels, edge_map, [0, 0], [2, 2])
+    blocked = [
+        dist.block_edge_map_stats(labels, edge_map, [0, 0], [1, 2]),
+        dist.block_edge_map_stats(labels, edge_map, [1, 0], [1, 2]),
+    ]
+    for per_block in ([single], blocked):
+        acc = dist.empty_edge_stats(graph.number_of_edges)
+        for be, bs in per_block:
+            acc = dist.merge_block_edge_stats(graph, acc, be, bs)
+        features = dist.finalize_edge_features(acc, compute_complex_features=True)
+        np.testing.assert_allclose(features[0, 1], 0.5)
+
+
 # --------------------------------------------------------------------------
 # Validation
 # --------------------------------------------------------------------------
@@ -377,3 +443,56 @@ def test_unsupported_label_dtype_raises():
     labels = np.zeros((5, 5), dtype=np.uint8)
     with pytest.raises(TypeError):
         dist.block_region_adjacency_edges(labels, [0, 0], [5, 5])
+
+
+def test_non_integer_owned_box_raises():
+    labels = np.zeros((6, 6), dtype=np.uint32)
+    with pytest.raises(TypeError, match="own_begin must contain integers"):
+        dist.block_region_adjacency_edges(labels, [0.9, 0], [3, 3])
+    with pytest.raises(TypeError, match="own_shape must contain integers"):
+        dist.block_region_adjacency_edges(labels, [0, 0], [1.9, 3])
+
+
+def test_non_integer_offsets_raise():
+    labels = np.zeros((6, 6), dtype=np.uint32)
+    affinities = np.zeros((1, 6, 6), dtype=np.float64)
+    with pytest.raises(TypeError, match="offsets must contain integers"):
+        dist.block_affinity_stats(labels, affinities, [[0.9, 0]], [0, 0], [6, 6])
+
+
+def test_merge_edges_rejects_negative_and_float_edges():
+    with pytest.raises(ValueError, match="negative"):
+        dist.merge_edges(np.array([[-1, 2]], dtype=np.int64))
+    with pytest.raises(ValueError, match="negative"):
+        dist.merge_edges([np.array([[0, 1]], dtype=np.uint64), [[-1, 2]]])
+    with pytest.raises(TypeError, match="integer dtype"):
+        dist.merge_edges(np.array([[0.5, 2.0]], dtype=np.float64))
+    # non-negative signed integers are fine
+    merged = dist.merge_edges(np.array([[1, 0]], dtype=np.int64))
+    np.testing.assert_array_equal(merged, np.array([[0, 1]], dtype=np.uint64))
+
+
+@pytest.mark.parametrize("dtype", [np.int32, np.int64])
+@pytest.mark.parametrize("number_of_threads", [1, 4])
+def test_negative_labels_raise(dtype, number_of_threads):
+    # A negative label must raise a normal Python exception from every block
+    # scanner, also when the check fires inside a worker thread (this used to
+    # terminate the process via an unhandled exception in std::thread).
+    labels = np.zeros((8, 8), dtype=dtype)
+    labels[3, 3] = -1
+    edge_map = np.zeros((8, 8), dtype=np.float64)
+    affinities = np.zeros((2, 8, 8), dtype=np.float64)
+
+    with pytest.raises(ValueError, match="negative"):
+        dist.block_region_adjacency_edges(
+            labels, [0, 0], [8, 8], number_of_threads=number_of_threads
+        )
+    with pytest.raises(ValueError, match="negative"):
+        dist.block_edge_map_stats(
+            labels, edge_map, [0, 0], [8, 8], number_of_threads=number_of_threads
+        )
+    with pytest.raises(ValueError, match="negative"):
+        dist.block_affinity_stats(
+            labels, affinities, [[1, 0], [0, 1]], [0, 0], [8, 8],
+            number_of_threads=number_of_threads,
+        )
diff --git a/tests/graph/test_rag_features.py b/tests/graph/test_rag_features.py
index 77a2bd9..d672ff1 100644
--- a/tests/graph/test_rag_features.py
+++ b/tests/graph/test_rag_features.py
@@ -66,6 +66,16 @@ def test_edge_map_features_complex():
     )
 
 
+def test_edge_map_features_complex_std_stable_for_large_baseline():
+    # Values 1e8 and 1e8 + 1 have std 0.5; the naive sum-of-squares formula
+    # returns 0.0 here due to catastrophic cancellation.
+    labels = np.array([[0, 1], [0, 1]], dtype=np.uint32)
+    edge_map = np.array([[1e8, 1e8], [1e8 + 1, 1e8 + 1]], dtype=np.float64)
+    rag = bic.graph.region_adjacency_graph(labels)
+    features = bic.graph.features.edge_map_features_complex(rag, labels, edge_map)
+    np.testing.assert_allclose(features[0, 2], 0.5)
+
+
 def test_affinity_features_simple():
     labels = np.array([[1, 1, 2], [1, 3, 2]], dtype=np.int64)
     rag = bic.graph.region_adjacency_graph(labels)
diff --git a/tests/graph/test_rag_lifted_features.py b/tests/graph/test_rag_lifted_features.py
index 24707b7..0e98760 100644
--- a/tests/graph/test_rag_lifted_features.py
+++ b/tests/graph/test_rag_lifted_features.py
@@ -219,7 +219,6 @@ def test_lifted_affinity_features_skips_local_hits():
         ],
         dtype=np.uint32,
     )
-    rag = bic.graph.region_adjacency_graph(labels)
     # Offset (2, 2): in a 3x3 grid only (0, 0) -> (2, 2), both label 0.
     # So no diff, no accumulation. We pass a lifted edge that doesn't exist
     # in the data to make sure no accumulation happens.
diff --git a/tests/graph/test_undirected_graph.py b/tests/graph/test_undirected_graph.py
index bdd489a..51286ee 100644
--- a/tests/graph/test_undirected_graph.py
+++ b/tests/graph/test_undirected_graph.py
@@ -157,6 +157,37 @@ def test_undirected_graph_freeze_is_callable_after_inserts_and_after_reads():
     )
 
 
+def test_from_edges_and_from_unique_edges_return_python_subclass():
+    # from_unique_edges requires lexicographically sorted, canonical edges.
+    uvs = np.array([[0, 1], [0, 2], [1, 2]], dtype=np.uint64)
+    core = bic._core.UndirectedGraph.from_unique_edges(3, uvs)
+
+    for graph in (
+        bic.graph.UndirectedGraph.from_unique_edges(3, uvs),
+        bic.graph.UndirectedGraph.from_edges(3, uvs),
+    ):
+        # The classmethods construct the Python subclass (the core statics
+        # return the raw _core type), so the convenience wrappers work.
+        assert isinstance(graph, bic.graph.UndirectedGraph)
+        assert graph.number_of_nodes == core.number_of_nodes
+        assert graph.number_of_edges == core.number_of_edges
+        np.testing.assert_array_equal(graph.uv_ids(), core.uv_ids())
+        np.testing.assert_array_equal(
+            graph.find_edges([[0, 1], [2, 1]]), np.array([0, 2], dtype=np.int64)
+        )
+
+
+def test_from_edges_deduplicates_but_from_unique_edges_requires_unique():
+    uvs_with_duplicate = np.array([[0, 1], [1, 0], [1, 2]], dtype=np.uint64)
+    graph = bic.graph.UndirectedGraph.from_edges(3, uvs_with_duplicate)
+    assert graph.number_of_edges == 2
+
+    with pytest.raises(ValueError):
+        bic.graph.UndirectedGraph.from_unique_edges(
+            3, np.array([[0, 0]], dtype=np.uint64)
+        )
+
+
 def test_undirected_graph_rejects_unrepresentable_node_count():
     with pytest.raises(OverflowError, match="CSR offsets"):
         bic.graph.UndirectedGraph(np.iinfo(np.uint64).max)
diff --git a/tests/label_multiset/test_merger.py b/tests/label_multiset/test_merger.py
index 58e568e..cc964b2 100644
--- a/tests/label_multiset/test_merger.py
+++ b/tests/label_multiset/test_merger.py
@@ -1,9 +1,7 @@
 import numpy as np
 
-from bioimage_cpp._core import Blocking
 from bioimage_cpp.label_multiset import (
     MultisetMerger,
-    downsample_multiset,
     multiset_from_labels,
 )
 
diff --git a/tests/label_multiset/test_read_subset.py b/tests/label_multiset/test_read_subset.py
index fad1e61..ac841d9 100644
--- a/tests/label_multiset/test_read_subset.py
+++ b/tests/label_multiset/test_read_subset.py
@@ -1,5 +1,4 @@
 import numpy as np
-import pytest
 
 from bioimage_cpp.label_multiset import read_subset