Address review finding from GPT 5.6 Sol

MIGRATION_GUIDE.md

@@ -2286,9 +2286,10 @@ Important differences:
   per-point radius is dynamic (the distance value at each point), matching
   nifty; there is no fixed-radius mode.
 - The algorithm is otherwise identical to nifty, including its float
-  arithmetic, so results match element-for-element. It uses an O(N²)
-  pairwise distance matrix; threshold the distance map first to keep the
-  candidate count modest.
+  arithmetic, so results match element-for-element. It is O(N²) in time but
+  uses only O(`number_of_threads` · N) auxiliary memory (no dense N×N matrix),
+  and `number_of_threads` parallelizes the pairwise evaluation. Still threshold
+  the distance map first to keep the candidate count modest.
 
 ## scikit-fmm
 

include/bioimage_cpp/affinities/compute_affinities.hxx

@@ -1,6 +1,7 @@
 #pragma once
 
 #include "bioimage_cpp/array_view.hxx"
+#include "bioimage_cpp/detail/grid.hxx"
 #include "bioimage_cpp/detail/threading.hxx"
 
 #include <algorithm>
@@ -66,9 +67,9 @@ void compute_affinities_2d(
                 }
 
                 // Sub-rectangle where (y, x) AND (y+dy, x+dx) are both in bounds.
-                const auto y_begin = std::max<std::ptrdiff_t>(0, -dy);
+                const auto y_begin = ::bioimage_cpp::detail::axis_begin_offset(dy, height);
                 const auto y_end = height - std::max<std::ptrdiff_t>(0, dy);
-                const auto x_begin = std::max<std::ptrdiff_t>(0, -dx);
+                const auto x_begin = ::bioimage_cpp::detail::axis_begin_offset(dx, width);
                 const auto x_end = width - std::max<std::ptrdiff_t>(0, dx);
                 if (y_begin >= y_end || x_begin >= x_end) {
                     continue;
@@ -143,11 +144,11 @@ void compute_affinities_3d(
                     std::fill_n(mask_channel, volume, std::uint8_t{0});
                 }
 
-                const auto z_begin = std::max<std::ptrdiff_t>(0, -dz);
+                const auto z_begin = ::bioimage_cpp::detail::axis_begin_offset(dz, depth);
                 const auto z_end = depth - std::max<std::ptrdiff_t>(0, dz);
-                const auto y_begin = std::max<std::ptrdiff_t>(0, -dy);
+                const auto y_begin = ::bioimage_cpp::detail::axis_begin_offset(dy, height);
                 const auto y_end = height - std::max<std::ptrdiff_t>(0, dy);
-                const auto x_begin = std::max<std::ptrdiff_t>(0, -dx);
+                const auto x_begin = ::bioimage_cpp::detail::axis_begin_offset(dx, width);
                 const auto x_end = width - std::max<std::ptrdiff_t>(0, dx);
                 if (z_begin >= z_end || y_begin >= y_end || x_begin >= x_end) {
                     continue;

include/bioimage_cpp/affinities/compute_embedding_distances.hxx

@@ -1,6 +1,7 @@
 #pragma once
 
 #include "bioimage_cpp/array_view.hxx"
+#include "bioimage_cpp/detail/grid.hxx"
 #include "bioimage_cpp/detail/threading.hxx"
 
 #include <algorithm>
@@ -108,9 +109,9 @@ void compute_embedding_distances_2d(
                     distances.data + static_cast<std::ptrdiff_t>(oi) * plane;
                 std::fill_n(dist_channel, plane, ValueT{0});
 
-                const auto y_begin = std::max<std::ptrdiff_t>(0, -dy);
+                const auto y_begin = ::bioimage_cpp::detail::axis_begin_offset(dy, height);
                 const auto y_end = height - std::max<std::ptrdiff_t>(0, dy);
-                const auto x_begin = std::max<std::ptrdiff_t>(0, -dx);
+                const auto x_begin = ::bioimage_cpp::detail::axis_begin_offset(dx, width);
                 const auto x_end = width - std::max<std::ptrdiff_t>(0, dx);
                 if (y_begin >= y_end || x_begin >= x_end) {
                     continue;
@@ -180,11 +181,11 @@ void compute_embedding_distances_3d(
                     distances.data + static_cast<std::ptrdiff_t>(oi) * volume;
                 std::fill_n(dist_channel, volume, ValueT{0});
 
-                const auto z_begin = std::max<std::ptrdiff_t>(0, -dz);
+                const auto z_begin = ::bioimage_cpp::detail::axis_begin_offset(dz, depth);
                 const auto z_end = depth - std::max<std::ptrdiff_t>(0, dz);
-                const auto y_begin = std::max<std::ptrdiff_t>(0, -dy);
+                const auto y_begin = ::bioimage_cpp::detail::axis_begin_offset(dy, height);
                 const auto y_end = height - std::max<std::ptrdiff_t>(0, dy);
-                const auto x_begin = std::max<std::ptrdiff_t>(0, -dx);
+                const auto x_begin = ::bioimage_cpp::detail::axis_begin_offset(dx, width);
                 const auto x_end = width - std::max<std::ptrdiff_t>(0, dx);
                 if (z_begin >= z_end || y_begin >= y_end || x_begin >= x_end) {
                     continue;

include/bioimage_cpp/detail/checked_arithmetic.hxx

@@ -0,0 +1,53 @@
+#pragma once
+
+#include <cstddef>
+#include <cstdint>
+#include <limits>
+#include <stdexcept>
+#include <string>
+
+namespace bioimage_cpp::detail {
+
+inline std::size_t checked_size_cast(const std::uint64_t value, const char *name) {
+    if constexpr (sizeof(std::size_t) < sizeof(std::uint64_t)) {
+        if (value > static_cast<std::uint64_t>(std::numeric_limits<std::size_t>::max())) {
+            throw std::overflow_error(std::string(name) + " exceeds size_t range");
+        }
+    }
+    return static_cast<std::size_t>(value);
+}
+
+inline std::size_t checked_size_add(
+    const std::size_t a,
+    const std::size_t b,
+    const char *name
+) {
+    if (b > std::numeric_limits<std::size_t>::max() - a) {
+        throw std::overflow_error(std::string(name) + " exceeds size_t range");
+    }
+    return a + b;
+}
+
+inline std::size_t checked_size_multiply(
+    const std::size_t a,
+    const std::size_t b,
+    const char *name
+) {
+    if (a != 0 && b > std::numeric_limits<std::size_t>::max() / a) {
+        throw std::overflow_error(std::string(name) + " exceeds size_t range");
+    }
+    return a * b;
+}
+
+inline std::uint64_t checked_u64_add(
+    const std::uint64_t a,
+    const std::uint64_t b,
+    const char *name
+) {
+    if (b > std::numeric_limits<std::uint64_t>::max() - a) {
+        throw std::overflow_error(std::string(name) + " exceeds uint64 range");
+    }
+    return a + b;
+}
+
+} // namespace bioimage_cpp::detail

include/bioimage_cpp/detail/grid.hxx

@@ -103,10 +103,32 @@ inline void valid_axis_range(
         const auto d = static_cast<std::size_t>(delta);
         hi = (d >= length) ? 0 : (length - d);
     } else {
-        const auto d = static_cast<std::size_t>(-delta);
+        // Avoid negating PTRDIFF_MIN. Converting first and subtracting from
+        // zero computes the magnitude in the unsigned domain.
+        const auto d = std::size_t{0} - static_cast<std::size_t>(delta);
         lo = (d >= length) ? length : d;
         hi = length;
     }
 }
 
+// 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
+// `ptrdiff_t` loop bounds and inner-loop codegen, while avoiding the undefined
+// behaviour of negating `delta == PTRDIFF_MIN`: magnitudes at least as large as
+// the axis (which make the whole offset channel empty) are clamped to `length`,
+// so the negation only runs when `-length < delta < 0` and is always safe.
+inline std::ptrdiff_t axis_begin_offset(
+    const std::ptrdiff_t delta,
+    const std::ptrdiff_t length
+) {
+    if (delta >= 0) {
+        return 0;
+    }
+    if (delta <= -length) {
+        return length;
+    }
+    return -delta;
+}
+
 } // namespace bioimage_cpp::detail

include/bioimage_cpp/detail/threading.hxx

@@ -2,6 +2,9 @@
 
 #include <algorithm>
 #include <cstddef>
+#include <exception>
+#include <mutex>
+#include <stdexcept>
 #include <thread>
 #include <utility>
 #include <vector>
@@ -27,34 +30,66 @@ 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::threads. The caller is
+// thread; threads 1..n_threads-1 run in parallel std::jthreads. The caller is
 // responsible for picking n_threads via normalize_thread_count and for the
 // thread safety of `chunk`.
+//
+// Exceptions thrown by `chunk` are captured rather than allowed to escape a
+// worker thread (which would call std::terminate). The first exception on any
+// thread is stored and rethrown on the calling thread after every worker has
+// been joined, so nanobind can translate it into a Python exception.
 template <class Chunk>
 void parallel_for_chunks(
     const std::size_t n_threads,
     const std::size_t number_of_work_items,
     Chunk &&chunk
 ) {
+    if (n_threads == 0) {
+        throw std::invalid_argument("parallel_for_chunks requires n_threads >= 1");
+    }
     const auto bounds = [&](const std::size_t thread_id) {
         const auto begin = thread_id * number_of_work_items / n_threads;
         const auto end = (thread_id + 1) * number_of_work_items / n_threads;
         return std::pair<std::size_t, std::size_t>{begin, end};
     };
 
-    std::vector<std::thread> threads;
+    std::exception_ptr first_exception;
+    std::mutex exception_mutex;
+    const auto guarded_chunk = [&] (
+        const std::size_t thread_id,
+        const std::size_t begin,
+        const std::size_t end
+    ) {
+        try {
+            chunk(thread_id, begin, end);
+        } catch (...) {
+            const std::lock_guard<std::mutex> lock(exception_mutex);
+            if (!first_exception) {
+                first_exception = std::current_exception();
+            }
+        }
+    };
+
+    // 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;
     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, &chunk]() {
-            chunk(thread_id, begin, end);
+        threads.emplace_back([thread_id, begin, end, &guarded_chunk]() {
+            guarded_chunk(thread_id, begin, end);
         });
     }
     const auto [begin, end] = bounds(0);
-    chunk(0, begin, end);
+    guarded_chunk(0, begin, end);
     for (auto &thread : threads) {
         thread.join();
     }
+    if (first_exception) {
+        std::rethrow_exception(first_exception);
+    }
 }
 
 } // namespace bioimage_cpp::detail