computational-cell-analytics
diff --git a/‎MIGRATION_GUIDE.md‎
Lines changed: 12 additions & 5 deletions b/‎MIGRATION_GUIDE.md‎
Lines changed: 12 additions & 5 deletions
diff --git a/‎include/bioimage_cpp/detail/grid.hxx‎
Lines changed: 95 additions & 0 deletions b/‎include/bioimage_cpp/detail/grid.hxx‎
Lines changed: 95 additions & 0 deletions
diff --git a/‎include/bioimage_cpp/detail/threading.hxx‎
Lines changed: 21 additions & 11 deletions b/‎include/bioimage_cpp/detail/threading.hxx‎
Lines changed: 21 additions & 11 deletions
@@ -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,19 +523,22 @@ 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)
 ```
 
 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.
 
@@ -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
 
@@ -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);