diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md index 2cbf5c3..3f6fbd8 100644 --- a/MIGRATION_GUIDE.md +++ b/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 diff --git a/include/bioimage_cpp/affinities/compute_affinities.hxx b/include/bioimage_cpp/affinities/compute_affinities.hxx index 92950d2..098038b 100644 --- a/include/bioimage_cpp/affinities/compute_affinities.hxx +++ b/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 @@ -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(0, -dy); + const auto y_begin = ::bioimage_cpp::detail::axis_begin_offset(dy, height); const auto y_end = height - std::max(0, dy); - const auto x_begin = std::max(0, -dx); + const auto x_begin = ::bioimage_cpp::detail::axis_begin_offset(dx, width); const auto x_end = width - std::max(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(0, -dz); + const auto z_begin = ::bioimage_cpp::detail::axis_begin_offset(dz, depth); const auto z_end = depth - std::max(0, dz); - const auto y_begin = std::max(0, -dy); + const auto y_begin = ::bioimage_cpp::detail::axis_begin_offset(dy, height); const auto y_end = height - std::max(0, dy); - const auto x_begin = std::max(0, -dx); + const auto x_begin = ::bioimage_cpp::detail::axis_begin_offset(dx, width); const auto x_end = width - std::max(0, dx); if (z_begin >= z_end || y_begin >= y_end || x_begin >= x_end) { continue; diff --git a/include/bioimage_cpp/affinities/compute_embedding_distances.hxx b/include/bioimage_cpp/affinities/compute_embedding_distances.hxx index 87984bf..2650bcd 100644 --- a/include/bioimage_cpp/affinities/compute_embedding_distances.hxx +++ b/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 @@ -108,9 +109,9 @@ void compute_embedding_distances_2d( distances.data + static_cast(oi) * plane; std::fill_n(dist_channel, plane, ValueT{0}); - const auto y_begin = std::max(0, -dy); + const auto y_begin = ::bioimage_cpp::detail::axis_begin_offset(dy, height); const auto y_end = height - std::max(0, dy); - const auto x_begin = std::max(0, -dx); + const auto x_begin = ::bioimage_cpp::detail::axis_begin_offset(dx, width); const auto x_end = width - std::max(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(oi) * volume; std::fill_n(dist_channel, volume, ValueT{0}); - const auto z_begin = std::max(0, -dz); + const auto z_begin = ::bioimage_cpp::detail::axis_begin_offset(dz, depth); const auto z_end = depth - std::max(0, dz); - const auto y_begin = std::max(0, -dy); + const auto y_begin = ::bioimage_cpp::detail::axis_begin_offset(dy, height); const auto y_end = height - std::max(0, dy); - const auto x_begin = std::max(0, -dx); + const auto x_begin = ::bioimage_cpp::detail::axis_begin_offset(dx, width); const auto x_end = width - std::max(0, dx); if (z_begin >= z_end || y_begin >= y_end || x_begin >= x_end) { continue; diff --git a/include/bioimage_cpp/detail/checked_arithmetic.hxx b/include/bioimage_cpp/detail/checked_arithmetic.hxx new file mode 100644 index 0000000..d910e6f --- /dev/null +++ b/include/bioimage_cpp/detail/checked_arithmetic.hxx @@ -0,0 +1,53 @@ +#pragma once + +#include +#include +#include +#include +#include + +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::numeric_limits::max())) { + throw std::overflow_error(std::string(name) + " exceeds size_t range"); + } + } + return static_cast(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::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::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::max() - a) { + throw std::overflow_error(std::string(name) + " exceeds uint64 range"); + } + return a + b; +} + +} // namespace bioimage_cpp::detail diff --git a/include/bioimage_cpp/detail/grid.hxx b/include/bioimage_cpp/detail/grid.hxx index 9e26ba0..954d405 100644 --- a/include/bioimage_cpp/detail/grid.hxx +++ b/include/bioimage_cpp/detail/grid.hxx @@ -103,10 +103,32 @@ inline void valid_axis_range( const auto d = static_cast(delta); hi = (d >= length) ? 0 : (length - d); } else { - const auto d = static_cast(-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(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 diff --git a/include/bioimage_cpp/detail/threading.hxx b/include/bioimage_cpp/detail/threading.hxx index 79df7db..61028f8 100644 --- a/include/bioimage_cpp/detail/threading.hxx +++ b/include/bioimage_cpp/detail/threading.hxx @@ -2,6 +2,9 @@ #include #include +#include +#include +#include #include #include #include @@ -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 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{begin, end}; }; - std::vector 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 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 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 diff --git a/include/bioimage_cpp/distance/detail/mesh_fast_marching.hxx b/include/bioimage_cpp/distance/detail/mesh_fast_marching.hxx index 5262611..c0b3d53 100644 --- a/include/bioimage_cpp/distance/detail/mesh_fast_marching.hxx +++ b/include/bioimage_cpp/distance/detail/mesh_fast_marching.hxx @@ -30,25 +30,74 @@ namespace bioimage_cpp::distance::detail { inline constexpr double kMeshFmmInfinity = std::numeric_limits::infinity(); -class MeshFastMarching { +class MeshTopology { public: - MeshFastMarching( + MeshTopology( const ConstArrayView &vertices, - const ConstArrayView &faces, - const ConstArrayView *speed + const ConstArrayView &faces ) : vertices_(vertices.data), n_vertices_(static_cast(vertices.shape[0])), faces_(faces.data), - n_faces_(static_cast(faces.shape[0])), - speed_(speed != nullptr ? speed->data : nullptr), - dist_(n_vertices_, kMeshFmmInfinity), - state_(n_vertices_, MeshState::Far), - heap_(n_vertices_) { + n_faces_(static_cast(faces.shape[0])) { build_vertex_face_incidence(); } [[nodiscard]] std::size_t size() const { return n_vertices_; } + [[nodiscard]] const double *vertices() const { return vertices_; } + [[nodiscard]] const std::int64_t *faces() const { return faces_; } + [[nodiscard]] const std::vector &face_offsets() const { + return face_offsets_; + } + [[nodiscard]] const std::vector &face_ids() const { return face_ids_; } + +private: + void build_vertex_face_incidence() { + face_offsets_.assign(n_vertices_ + 1, 0); + for (std::size_t f = 0; f < n_faces_; ++f) { + for (std::size_t corner = 0; corner < 3; ++corner) { + const std::int64_t v = faces_[f * 3 + corner]; + if (v < 0 || static_cast(v) >= n_vertices_) { + throw std::invalid_argument( + "faces contains vertex index " + std::to_string(v) + + " out of range [0, " + std::to_string(n_vertices_) + ")" + ); + } + ++face_offsets_[static_cast(v) + 1]; + } + } + for (std::size_t i = 1; i < face_offsets_.size(); ++i) { + face_offsets_[i] += face_offsets_[i - 1]; + } + face_ids_.resize(face_offsets_.back()); + std::vector insert_pos(face_offsets_.begin(), face_offsets_.end() - 1); + for (std::size_t f = 0; f < n_faces_; ++f) { + for (std::size_t corner = 0; corner < 3; ++corner) { + const auto v = static_cast(faces_[f * 3 + corner]); + face_ids_[insert_pos[v]++] = f; + } + } + } + + const double *vertices_; + std::size_t n_vertices_; + const std::int64_t *faces_; + std::size_t n_faces_; + std::vector face_offsets_; + std::vector face_ids_; +}; + +class MeshFastMarching { +public: + MeshFastMarching(const MeshTopology &topology, const ConstArrayView *speed) + : topology_(topology), + speed_(speed != nullptr ? speed->data : nullptr), + dist_(topology.size(), kMeshFmmInfinity), + state_(topology.size(), MeshState::Far), + heap_(topology.size()) { + } + + [[nodiscard]] std::size_t size() const { return topology_.size(); } [[nodiscard]] double distance(std::size_t vertex) const { return dist_[vertex]; } [[nodiscard]] const std::vector &distances() const { return dist_; } @@ -81,10 +130,13 @@ public: // and already reflects every other incident face's contribution from // the earlier freeze of that face's corner, so min(dist_[w], full // rescan) collapses to min(dist_[w], this-face candidate). - for (std::size_t k = face_offsets_[v]; k < face_offsets_[v + 1]; ++k) { - const std::size_t f = face_ids_[k]; + const auto &face_offsets = topology_.face_offsets(); + const auto &face_ids = topology_.face_ids(); + const auto *faces = topology_.faces(); + for (std::size_t k = face_offsets[v]; k < face_offsets[v + 1]; ++k) { + const std::size_t f = face_ids[k]; for (std::size_t corner = 0; corner < 3; ++corner) { - const std::size_t w = static_cast(faces_[f * 3 + corner]); + const std::size_t w = static_cast(faces[f * 3 + corner]); if (w == v || state_[w] == MeshState::Frozen) { continue; } @@ -96,7 +148,7 @@ public: std::array others{}; std::size_t count = 0; for (std::size_t c2 = 0; c2 < 3; ++c2) { - const auto u = static_cast(faces_[f * 3 + c2]); + const auto u = static_cast(faces[f * 3 + c2]); if (u != w) { others[count++] = u; } @@ -133,35 +185,6 @@ private: heap_.clear(); } - // Build the vertex -> incident-faces CSR (offsets + face ids), mirroring the - // count-then-scatter pattern in detail::mesh_smoothing::build_adjacency. - void build_vertex_face_incidence() { - face_offsets_.assign(n_vertices_ + 1, 0); - for (std::size_t f = 0; f < n_faces_; ++f) { - for (std::size_t corner = 0; corner < 3; ++corner) { - const std::int64_t v = faces_[f * 3 + corner]; - if (v < 0 || static_cast(v) >= n_vertices_) { - throw std::invalid_argument( - "faces contains vertex index " + std::to_string(v) + - " out of range [0, " + std::to_string(n_vertices_) + ")" - ); - } - ++face_offsets_[static_cast(v) + 1]; - } - } - for (std::size_t i = 1; i < face_offsets_.size(); ++i) { - face_offsets_[i] += face_offsets_[i - 1]; - } - face_ids_.resize(face_offsets_.back()); - std::vector insert_pos(face_offsets_.begin(), face_offsets_.end() - 1); - for (std::size_t f = 0; f < n_faces_; ++f) { - for (std::size_t corner = 0; corner < 3; ++corner) { - const auto v = static_cast(faces_[f * 3 + corner]); - face_ids_[insert_pos[v]++] = f; - } - } - } - [[nodiscard]] double slowness(std::size_t vertex) const { if (speed_ == nullptr) { return 1.0; @@ -171,8 +194,8 @@ private: } [[nodiscard]] double edge_length(std::size_t u, std::size_t v) const { - const double *pu = vertices_ + u * 3; - const double *pv = vertices_ + v * 3; + const double *pu = topology_.vertices() + u * 3; + const double *pv = topology_.vertices() + v * 3; const double dx = pu[0] - pv[0]; const double dy = pu[1] - pv[1]; const double dz = pu[2] - pv[2]; @@ -192,9 +215,9 @@ private: return edge_cand; } - const double *pc = vertices_ + c * 3; - const double *pa = vertices_ + a * 3; - const double *pb = vertices_ + b * 3; + const double *pc = topology_.vertices() + c * 3; + const double *pa = topology_.vertices() + a * 3; + const double *pb = topology_.vertices() + b * 3; double dot = 0.0; for (std::size_t d = 0; d < 3; ++d) { dot += (pa[d] - pc[d]) * (pb[d] - pc[d]); @@ -229,15 +252,9 @@ private: return edge_cand; } - const double *vertices_; - std::size_t n_vertices_; - const std::int64_t *faces_; - std::size_t n_faces_; + const MeshTopology &topology_; const double *speed_; - std::vector face_offsets_; - std::vector face_ids_; - std::vector dist_; std::vector state_; bioimage_cpp::detail::DenseIndexedHeap> heap_; diff --git a/include/bioimage_cpp/distance/geodesic_mesh.hxx b/include/bioimage_cpp/distance/geodesic_mesh.hxx index fe0b923..0fc596b 100644 --- a/include/bioimage_cpp/distance/geodesic_mesh.hxx +++ b/include/bioimage_cpp/distance/geodesic_mesh.hxx @@ -62,7 +62,8 @@ inline void geodesic_distance_field_mesh( ArrayView &distances, std::size_t /*n_threads*/ = 1 ) { - detail::MeshFastMarching solver(vertices, faces, speed); + const detail::MeshTopology topology(vertices, faces); + detail::MeshFastMarching solver(topology, speed); const auto source_indices = detail::vertex_indices(sources, solver.size(), "sources"); solver.solve(source_indices); @@ -84,9 +85,8 @@ inline void geodesic_distances_mesh( ArrayView &distances, std::size_t n_threads = 1 ) { - // Build one solver up front to validate the mesh and the point indices. - detail::MeshFastMarching probe(vertices, faces, speed); - const auto point_indices = detail::vertex_indices(points, probe.size(), "points"); + const detail::MeshTopology topology(vertices, faces); + const auto point_indices = detail::vertex_indices(points, topology.size(), "points"); const std::size_t n = point_indices.size(); double *out = distances.data; @@ -98,7 +98,7 @@ inline void geodesic_distances_mesh( bioimage_cpp::detail::parallel_for_chunks( threads, n, [&](std::size_t /*thread_id*/, std::size_t begin, std::size_t end) { - detail::MeshFastMarching solver(vertices, faces, speed); + detail::MeshFastMarching solver(topology, speed); for (std::size_t i = begin; i < end; ++i) { solver.solve({point_indices[i]}); for (std::size_t j = 0; j < n; ++j) { diff --git a/include/bioimage_cpp/distance/non_maximum_distance_suppression.hxx b/include/bioimage_cpp/distance/non_maximum_distance_suppression.hxx index 3a0806c..cf88cda 100644 --- a/include/bioimage_cpp/distance/non_maximum_distance_suppression.hxx +++ b/include/bioimage_cpp/distance/non_maximum_distance_suppression.hxx @@ -2,6 +2,7 @@ #include "bioimage_cpp/array_view.hxx" #include "bioimage_cpp/detail/grid.hxx" +#include "bioimage_cpp/detail/threading.hxx" #include #include @@ -46,12 +47,13 @@ inline std::ptrdiff_t point_to_flat( // (rather than comparing squared distances) keeps boundary ties identical to // nifty. // -// Complexity: O(N^2) time and O(N^2) memory for the symmetric distance matrix. +// Complexity: O(N^2) time and O(number_of_threads * N) auxiliary memory. template inline void non_maximum_distance_suppression( const ConstArrayView &distance_map, const ConstArrayView &points, - std::vector &kept_indices + std::vector &kept_indices, + const std::size_t number_of_threads = 1 ) { if (distance_map.ndim() < 1) { throw std::invalid_argument( @@ -92,49 +94,71 @@ inline void non_maximum_distance_suppression( point_dist[i] = distance_map.data[flat]; } - // Pairwise Euclidean distance matrix (symmetric, N x N). The float - // accumulation and float(sqrt(...)) match nifty bit-for-bit so the - // neighborhood test below produces identical boundary decisions. - std::vector pd(n * n, 0.0f); - for (std::size_t i = 0; i < n; ++i) { - const auto *row_i = - points.data + static_cast(i) * coord_ndim; - for (std::size_t j = i + 1; j < n; ++j) { - const auto *row_j = - points.data + static_cast(j) * coord_ndim; - float sum_sq = 0.0f; - for (std::ptrdiff_t d = 0; d < coord_ndim; ++d) { - const float diff = static_cast(row_i[d]) - - static_cast(row_j[d]); - sum_sq += diff * diff; + const auto threads = bioimage_cpp::detail::normalize_thread_count( + number_of_threads, n + ); + constexpr auto no_candidate = std::numeric_limits::max(); + std::vector> local_best( + threads, std::vector(n, no_candidate) + ); + const auto consider = [&](std::size_t ¤t, const std::size_t candidate) { + const float value = point_dist[candidate]; + if (current == no_candidate) { + if (value > -std::numeric_limits::infinity()) { + current = candidate; } - const auto val = static_cast(std::sqrt(static_cast(sum_sq))); - pd[i * n + j] = val; - pd[j * n + i] = val; + return; } - } + const float current_value = point_dist[current]; + if (value > current_value || (value == current_value && candidate < current)) { + current = candidate; + } + }; + + bioimage_cpp::detail::parallel_for_chunks( + threads, threads, + [&](const std::size_t thread_id, const std::size_t, const std::size_t) { + auto &best = local_best[thread_id]; + for (std::size_t i = thread_id; i < n; i += threads) { + const auto *row_i = + points.data + static_cast(i) * coord_ndim; + for (std::size_t j = i + 1; j < n; ++j) { + const auto *row_j = + points.data + static_cast(j) * coord_ndim; + float sum_sq = 0.0f; + for (std::ptrdiff_t d = 0; d < coord_ndim; ++d) { + const float diff = static_cast(row_i[d]) - + static_cast(row_j[d]); + sum_sq += diff * diff; + } + const auto distance = static_cast( + std::sqrt(static_cast(sum_sq)) + ); + if (!(distance > point_dist[i])) { + consider(best[i], j); + } + if (!(distance > point_dist[j])) { + consider(best[j], i); + } + } + } + } + ); - // For each point, scan all neighbors within its dynamic radius and keep - // the one with the largest distance_map value. Strict `>` ensures the - // first index encountered wins ties, matching nifty's behavior. std::vector bests; bests.reserve(n); for (std::size_t i = 0; i < n; ++i) { - const float d_i = point_dist[i]; - float best_val = -std::numeric_limits::infinity(); - std::size_t best_idx = i; - const auto *pd_row = pd.data() + i * n; - for (std::size_t j = 0; j < n; ++j) { - if (pd_row[j] > d_i) { - continue; - } - const float dj = point_dist[j]; - if (dj > best_val) { - best_val = dj; - best_idx = j; + std::size_t best_idx = no_candidate; + if (!(0.0f > point_dist[i])) { + consider(best_idx, i); + } + for (std::size_t thread = 0; thread < threads; ++thread) { + const auto candidate = local_best[thread][i]; + if (candidate != no_candidate) { + consider(best_idx, candidate); } } - bests.push_back(best_idx); + bests.push_back(best_idx == no_candidate ? i : best_idx); } std::sort(bests.begin(), bests.end()); diff --git a/include/bioimage_cpp/flow/flow_density.hxx b/include/bioimage_cpp/flow/flow_density.hxx index 60cbf2d..6eced35 100644 --- a/include/bioimage_cpp/flow/flow_density.hxx +++ b/include/bioimage_cpp/flow/flow_density.hxx @@ -118,6 +118,20 @@ inline std::ptrdiff_t round_to_flat_index( return flat; } +template +inline bool position_is_in_mask( + const std::array &position, + const GridLayout &grid, + const std::uint8_t *mask +) { + for (std::size_t axis = 0; axis < D; ++axis) { + if (position[axis] < 0.0f || position[axis] > grid.upper[axis]) { + return false; + } + } + return mask[round_to_flat_index(position, grid)] != 0; +} + } // namespace detail enum class IntegrationMethod { @@ -213,14 +227,6 @@ void compute_flow_density( auto &position = positions[i]; clip_position(position); - if (restrict_to_mask) { - const auto here = detail::round_to_flat_index(position, grid); - if (fg_mask.data[here] == 0) { - alive[i] = 0; - continue; - } - } - const auto corners = detail::compute_corners(position, grid); std::array step{}; for (std::size_t axis = 0; axis < D; ++axis) { @@ -250,9 +256,24 @@ void compute_flow_density( alive[i] = 0; continue; } + auto proposed = position; for (std::size_t axis = 0; axis < D; ++axis) { - position[axis] += dt * step[axis]; + proposed[axis] += dt * step[axis]; + } + // A particle whose proposed endpoint leaves the + // foreground is frozen at its last in-mask position + // (only the endpoint is mask-tested, not the RK2 + // midpoint). Every alive particle is seeded in the mask + // and only commits in-mask endpoints, so its current + // position is always a valid last position and no + // separate start-of-step mask test is needed. + if (restrict_to_mask && !detail::position_is_in_mask( + proposed, grid, fg_mask.data + )) { + alive[i] = 0; + continue; } + position = proposed; } } ); diff --git a/include/bioimage_cpp/graph/rag_coordinates.hxx b/include/bioimage_cpp/graph/rag_coordinates.hxx index cfe56a8..ea4fcd2 100644 --- a/include/bioimage_cpp/graph/rag_coordinates.hxx +++ b/include/bioimage_cpp/graph/rag_coordinates.hxx @@ -1,6 +1,7 @@ #pragma once #include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/label_cast.hxx" #include "bioimage_cpp/detail/grid.hxx" #include "bioimage_cpp/detail/threading.hxx" #include "bioimage_cpp/graph/region_adjacency_graph.hxx" @@ -33,9 +34,9 @@ void scan_contacts_2d( const auto row_offset = y * width; for (std::size_t x = 0; x < width; ++x) { const auto pixel = row_offset + x; - const auto u = static_cast(data[pixel]); + const auto u = bioimage_cpp::detail::checked_label_to_node(data[pixel]); if (x + 1 < width) { - const auto v = static_cast(data[pixel + 1]); + const auto v = bioimage_cpp::detail::checked_label_to_node(data[pixel + 1]); if (u != v) { const auto edge = rag.find_edge(u, v); if (edge >= 0) { @@ -44,7 +45,7 @@ void scan_contacts_2d( } } if (y + 1 < height) { - const auto v = static_cast(data[pixel + width]); + const auto v = bioimage_cpp::detail::checked_label_to_node(data[pixel + width]); if (u != v) { const auto edge = rag.find_edge(u, v); if (edge >= 0) { @@ -75,9 +76,9 @@ void scan_contacts_3d( const auto row_offset = slice_offset + y * width; for (std::size_t x = 0; x < width; ++x) { const auto pixel = row_offset + x; - const auto u = static_cast(data[pixel]); + const auto u = bioimage_cpp::detail::checked_label_to_node(data[pixel]); if (x + 1 < width) { - const auto v = static_cast(data[pixel + 1]); + const auto v = bioimage_cpp::detail::checked_label_to_node(data[pixel + 1]); if (u != v) { const auto edge = rag.find_edge(u, v); if (edge >= 0) { @@ -86,7 +87,7 @@ void scan_contacts_3d( } } if (y + 1 < height) { - const auto v = static_cast(data[pixel + width]); + const auto v = bioimage_cpp::detail::checked_label_to_node(data[pixel + width]); if (u != v) { const auto edge = rag.find_edge(u, v); if (edge >= 0) { @@ -95,7 +96,7 @@ void scan_contacts_3d( } } if (z + 1 < depth) { - const auto v = static_cast(data[pixel + slice_size]); + const auto v = bioimage_cpp::detail::checked_label_to_node(data[pixel + slice_size]); if (u != v) { const auto edge = rag.find_edge(u, v); if (edge >= 0) { diff --git a/include/bioimage_cpp/graph/undirected_graph.hxx b/include/bioimage_cpp/graph/undirected_graph.hxx index 1fac051..715cccd 100644 --- a/include/bioimage_cpp/graph/undirected_graph.hxx +++ b/include/bioimage_cpp/graph/undirected_graph.hxx @@ -1,10 +1,12 @@ #pragma once +#include "bioimage_cpp/detail/checked_arithmetic.hxx" #include "bioimage_cpp/detail/edge_hash.hxx" #include #include #include +#include #include #include #include @@ -89,7 +91,7 @@ public: const NodeId number_of_nodes = 0, const EdgeId reserve_number_of_edges = 0 ) { - number_of_nodes_ = number_of_nodes; + number_of_nodes_ = checked_node_count(number_of_nodes); edges_.clear(); edges_.reserve(static_cast(reserve_number_of_edges)); edge_lookup_.clear(); @@ -195,7 +197,11 @@ public: } [[nodiscard]] std::uint64_t serialization_size() const { - return 2 + 2 * number_of_edges(); + const auto n_edges = number_of_edges(); + if (n_edges > (std::numeric_limits::max() - 2) / 2) { + throw std::overflow_error("graph serialization size exceeds uint64 range"); + } + return 2 + 2 * n_edges; } // Force the CSR adjacency to be rebuilt if it's currently stale. Use this @@ -313,7 +319,7 @@ protected: const EdgeId reserve_edges, const EdgeId reserve_lookup ) - : number_of_nodes_(number_of_nodes) { + : number_of_nodes_(checked_node_count(number_of_nodes)) { edges_.reserve(static_cast(reserve_edges)); edge_lookup_.reserve(static_cast(reserve_lookup)); } @@ -358,10 +364,13 @@ protected: // typically 4-5x faster than the previous vector-of-vectors fill on // large grids. void rebuild_adjacency_from_edges() { - const auto n_nodes = static_cast(number_of_nodes_); - const auto data_size = 2 * edges_.size(); + const auto n_nodes = detail::checked_size_cast(number_of_nodes_, "number_of_nodes"); + const auto offsets_size = detail::checked_size_add(n_nodes, 1, "CSR offsets size"); + const auto data_size = detail::checked_size_multiply( + edges_.size(), 2, "CSR adjacency size" + ); // Pass 1: count per-node degree into the offsets buffer (shifted by 1). - adjacency_offsets_.assign(n_nodes + 1, 0); + adjacency_offsets_.assign(offsets_size, 0); for (const auto &edge : edges_) { ++adjacency_offsets_[static_cast(edge.first) + 1]; ++adjacency_offsets_[static_cast(edge.second) + 1]; @@ -416,6 +425,12 @@ protected: } private: + static NodeId checked_node_count(const NodeId number_of_nodes) { + const auto count = detail::checked_size_cast(number_of_nodes, "number_of_nodes"); + detail::checked_size_add(count, 1, "CSR offsets size"); + return number_of_nodes; + } + void ensure_adjacency_built() const { if (adjacency_dirty_) { const_cast(this)->rebuild_adjacency_from_edges(); diff --git a/include/bioimage_cpp/segmentation/relabel_sequential.hxx b/include/bioimage_cpp/segmentation/relabel_sequential.hxx index 4b40649..a0b9fc2 100644 --- a/include/bioimage_cpp/segmentation/relabel_sequential.hxx +++ b/include/bioimage_cpp/segmentation/relabel_sequential.hxx @@ -1,8 +1,12 @@ #pragma once #include "bioimage_cpp/array_view.hxx" +#include "bioimage_cpp/detail/checked_arithmetic.hxx" +#include #include +#include +#include #include #include #include @@ -18,6 +22,13 @@ struct RelabelSequentialMaps { // Relabel an integer array so all non-zero labels become consecutive starting // at offset, in sorted order. Label 0 is treated as background and always // maps to 0. Matches skimage.segmentation.relabel_sequential. +// +// The maps are dense NumPy lookup arrays (forward_map indexed by old label, +// inverse_map indexed by new label), so their sizes scale with `max(label)` +// and `offset`. Inputs whose required maps would overflow `size_t`, the label +// dtype, or (via the pre-sized inverse map) available memory are rejected with +// a clear exception rather than wrapping a size to zero (a heap overflow) or +// silently truncating a new label. template RelabelSequentialMaps relabel_sequential( const ConstArrayView &input, @@ -36,9 +47,14 @@ RelabelSequentialMaps relabel_sequential( ); RelabelSequentialMaps maps; - maps.inverse_map.assign(static_cast(offset), T{0}); if (n == 0) { + // No labels: inverse_map is just the `offset`-length background prefix, + // forward_map is empty. + maps.inverse_map.assign( + detail::checked_size_cast(static_cast(offset), "relabel offset"), + T{0} + ); return maps; } @@ -52,9 +68,16 @@ RelabelSequentialMaps relabel_sequential( } // The forward LUT doubles as a presence bitmap in pass 2 (any non-zero - // value means "present"), then gets rewritten with the final new labels - // by the sorted scan below. - maps.forward_map.assign(static_cast(max_value) + 1u, T{0}); + // value means "present"), then gets rewritten with the final new labels by + // the sorted scan below. Its size is `max_value + 1`; the checked add turns + // a maximal label (e.g. UINT64_MAX, whose `+ 1` would wrap the size to + // zero and make pass 2 write out of bounds) into a clean overflow error. + const auto forward_size = detail::checked_size_add( + detail::checked_size_cast(static_cast(max_value), "relabel max label"), + 1, + "relabel forward map size" + ); + maps.forward_map.assign(forward_size, T{0}); // Pass 2: mark each input value as present in the LUT. for (std::ptrdiff_t i = 0; i < n; ++i) { @@ -63,18 +86,46 @@ RelabelSequentialMaps relabel_sequential( // 0 is reserved background and maps to 0 even if it appeared in input. maps.forward_map[0] = T{0}; - // Sorted-order assignment: walk 1..max_value, replace presence sentinel - // with the next sequential label, and push the old label onto the - // inverse map (already pre-sized to `offset` zeros). + // Sorted-order assignment: walk 1..max_value, replace the presence sentinel + // with the next sequential label, and collect the old labels. The new + // labels occupy [offset, max(T)]; `capacity` is how many fit. Guarding the + // count here — before the (potentially large) inverse-map allocation — + // rejects an out-of-range `offset` cleanly instead of wrapping `new_label` + // past the dtype maximum or pre-sizing a multi-gigabyte inverse map. + const std::uint64_t capacity = + static_cast(std::numeric_limits::max()) - + static_cast(offset) + 1; + std::vector old_labels; T new_label = offset; for (std::size_t v = 1; v <= static_cast(max_value); ++v) { if (maps.forward_map[v] != T{0}) { + if (static_cast(old_labels.size()) >= capacity) { + throw std::overflow_error( + "relabel_sequential: offset + number of distinct labels " + "exceeds the label dtype range" + ); + } maps.forward_map[v] = new_label; - maps.inverse_map.push_back(static_cast(v)); + old_labels.push_back(static_cast(v)); ++new_label; } } + // inverse_map = `offset` background zeros followed by the old labels, in + // ascending new-label order. + const auto prefix = detail::checked_size_cast( + static_cast(offset), "relabel offset" + ); + const auto inverse_size = detail::checked_size_add( + prefix, old_labels.size(), "relabel inverse map size" + ); + maps.inverse_map.assign(inverse_size, T{0}); + std::copy( + old_labels.begin(), + old_labels.end(), + maps.inverse_map.begin() + static_cast(prefix) + ); + // Pass 3: apply the forward LUT to write the relabeled output. for (std::ptrdiff_t i = 0; i < n; ++i) { out.data[i] = maps.forward_map[static_cast(input.data[i])]; diff --git a/include/bioimage_cpp/transformation/affine.hxx b/include/bioimage_cpp/transformation/affine.hxx index c7391f8..a5835e6 100644 --- a/include/bioimage_cpp/transformation/affine.hxx +++ b/include/bioimage_cpp/transformation/affine.hxx @@ -498,11 +498,23 @@ inline double cubic_3d( // Dispatch to the per-voxel sampler for the requested interpolation order. Shared by the affine and // map_coordinates kernels so the interpolation backend lives in exactly one place. `order` must be in // 0..5 (validated by the public entry points before the sampling loop). -template +// +// `CheckFinite` gates a per-coordinate `std::isfinite` guard. map_coordinates +// passes user-supplied coordinates that may be non-finite and instantiates it +// with `true`. Affine coordinates are `matrix * output_grid` with a +// finite-validated matrix, so the affine callers use the default `false` and +// the check is elided at compile time — keeping the affine sampling hot loop +// identical to the pre-fix code. +template inline double sample_2d( const T *data, std::ptrdiff_t in_h, std::ptrdiff_t in_w, double cy, double cx, const int order, double fill ) { + if constexpr (CheckFinite) { + if (!std::isfinite(cy) || !std::isfinite(cx)) { + return fill; + } + } switch (order) { case 0: return nearest_2d(data, in_h, in_w, cy, cx, fill); case 1: return linear_2d(data, in_h, in_w, cy, cx, fill); @@ -513,12 +525,17 @@ inline double sample_2d( } } -template +template inline double sample_3d( const T *data, std::ptrdiff_t in_d, std::ptrdiff_t in_h, std::ptrdiff_t in_w, double cz, double cy, double cx, const int order, double fill ) { + if constexpr (CheckFinite) { + if (!std::isfinite(cz) || !std::isfinite(cy) || !std::isfinite(cx)) { + return fill; + } + } switch (order) { case 0: return nearest_3d(data, in_d, in_h, in_w, cz, cy, cx, fill); case 1: return linear_3d(data, in_d, in_h, in_w, cz, cy, cx, fill); diff --git a/include/bioimage_cpp/transformation/coordinate.hxx b/include/bioimage_cpp/transformation/coordinate.hxx index 7b5eec9..104d783 100644 --- a/include/bioimage_cpp/transformation/coordinate.hxx +++ b/include/bioimage_cpp/transformation/coordinate.hxx @@ -77,7 +77,8 @@ void map_coordinates_2d( T *out_ptr = output.data; for (std::ptrdiff_t p = 0; p < n_out; ++p) { - const double value = detail::sample_2d(in_data, in_h, in_w, cy[p], cx[p], order, fill); + const double value = + detail::sample_2d(in_data, in_h, in_w, cy[p], cx[p], order, fill); out_ptr[p] = detail::to_output(value); } } @@ -120,8 +121,8 @@ void map_coordinates_3d( T *out_ptr = output.data; for (std::ptrdiff_t p = 0; p < n_out; ++p) { - const double value = detail::sample_3d(in_data, in_d, in_h, in_w, - cz[p], cy[p], cx[p], order, fill); + const double value = detail::sample_3d(in_data, in_d, in_h, in_w, + cz[p], cy[p], cx[p], order, fill); out_ptr[p] = detail::to_output(value); } } diff --git a/src/bindings/distance.cxx b/src/bindings/distance.cxx index 087fdde..fc3cfeb 100644 --- a/src/bindings/distance.cxx +++ b/src/bindings/distance.cxx @@ -200,8 +200,6 @@ nb::ndarray non_maximum_distance_suppression_im nb::ndarray points, const std::size_t n_threads ) { - (void)n_threads; // Reserved for future parallelization; single-threaded. - if (distance_map.ndim() == 0) { throw std::invalid_argument("distance_map must have ndim >= 1, got ndim=0"); } @@ -253,7 +251,9 @@ nb::ndarray non_maximum_distance_suppression_im std::vector kept_indices; { nb::gil_scoped_release release; - distance::non_maximum_distance_suppression(map_view, points_view, kept_indices); + distance::non_maximum_distance_suppression( + map_view, points_view, kept_indices, n_threads + ); } const std::size_t n_kept = kept_indices.size(); diff --git a/src/bindings/filters.cxx b/src/bindings/filters.cxx index 874d56b..125d022 100644 --- a/src/bindings/filters.cxx +++ b/src/bindings/filters.cxx @@ -4,6 +4,7 @@ #include +#include #include #include #include @@ -26,7 +27,7 @@ void require_ndim(const ConstImage &image, int expected, const char *function) { } void require_positive_sigma(double sigma, const char *name, const char *function) { - if (!(sigma > 0.0)) { + if (!(std::isfinite(sigma) && sigma > 0.0)) { throw std::invalid_argument( std::string(function) + ": " + name + " must be positive, got " + std::to_string(sigma) @@ -44,7 +45,7 @@ void require_order(int order, const char *name, const char *function) { } void require_non_negative_window(double window_size, const char *function) { - if (window_size < 0.0) { + if (!std::isfinite(window_size) || window_size < 0.0) { throw std::invalid_argument( std::string(function) + ": window_size must be >= 0 (0 selects the " "default), got " + std::to_string(window_size) diff --git a/src/bindings/flow.cxx b/src/bindings/flow.cxx index 5b6421d..ffc950f 100644 --- a/src/bindings/flow.cxx +++ b/src/bindings/flow.cxx @@ -99,6 +99,17 @@ DensityArray compute_flow_density_t( if (number_of_threads < 1) { throw std::invalid_argument("number_of_threads must be >= 1"); } + // Single authoritative finiteness check over the (contiguous) flow buffer; + // the Python wrapper deliberately does not repeat this scan. + std::size_t flow_size = D; + for (std::size_t axis = 0; axis < D; ++axis) { + flow_size *= flow.shape(axis + 1); + } + for (std::size_t index = 0; index < flow_size; ++index) { + if (!std::isfinite(flow.data()[index])) { + throw std::invalid_argument("flow must contain only finite values"); + } + } std::vector out_shape(D); std::vector view_shape(D); diff --git a/src/bindings/graph.cxx b/src/bindings/graph.cxx index b47df57..26d929b 100644 --- a/src/bindings/graph.cxx +++ b/src/bindings/graph.cxx @@ -628,12 +628,19 @@ Graph graph_from_unique_edges(const std::uint64_t number_of_nodes, ConstUInt64Ar for (std::size_t index = 0; index < uvs.shape(0); ++index) { const auto u = in[2 * index]; const auto v = in[2 * index + 1]; - if (u == v) { - throw std::invalid_argument("self edges are not supported"); + if (u >= v) { + throw std::invalid_argument( + "uvs must contain canonical edges with u < v" + ); } if (u >= number_of_nodes || v >= number_of_nodes) { throw std::out_of_range("edge endpoint exceeds number_of_nodes"); } + if (!edges.empty() && !(edges.back() < Graph::Edge{u, v})) { + throw std::invalid_argument( + "uvs must be strictly lexicographically sorted and duplicate-free" + ); + } edges.emplace_back(u, v); } return Graph::from_sorted_unique_edges(number_of_nodes, std::move(edges)); diff --git a/src/bindings/transformation.cxx b/src/bindings/transformation.cxx index 716fecb..d2f2deb 100644 --- a/src/bindings/transformation.cxx +++ b/src/bindings/transformation.cxx @@ -7,6 +7,7 @@ #include +#include #include #include #include @@ -61,6 +62,18 @@ OutputArray affine_transform_t( const auto input_shape = shape_of(input); const auto input_strides = detail::c_order_strides(input_shape); const auto matrix_shape = shape_of(matrix); + // The affine sampler assumes finite coordinates (matrix * grid), so the + // matrix is validated finite here, at the single binding boundary, rather + // than per output pixel. + std::size_t matrix_size = 1; + for (const auto extent : matrix_shape) { + matrix_size *= static_cast(extent); + } + for (std::size_t index = 0; index < matrix_size; ++index) { + if (!std::isfinite(matrix.data()[index])) { + throw std::invalid_argument("matrix must contain only finite values"); + } + } const auto matrix_strides = detail::c_order_strides(matrix_shape); const auto starts_shape = shape_of(starts); const auto starts_strides = detail::c_order_strides(starts_shape); diff --git a/src/bioimage_cpp/_validation.py b/src/bioimage_cpp/_validation.py new file mode 100644 index 0000000..4489b17 --- /dev/null +++ b/src/bioimage_cpp/_validation.py @@ -0,0 +1,72 @@ +"""Small, dependency-free validators shared by public Python wrappers.""" + +from __future__ import annotations + +import operator + +import numpy as np + + +def strict_index(value, name: str, *, minimum: int | None = None, + maximum: int | None = None) -> int: + if isinstance(value, (bool, np.bool_)): + raise TypeError(f"{name} must be an integer, got bool") + try: + result = operator.index(value) + except TypeError as error: + raise TypeError(f"{name} must be an integer, got {value!r}") from error + if minimum is not None and result < minimum: + raise ValueError(f"{name} must be >= {minimum}, got {result}") + if maximum is not None and result > maximum: + raise ValueError(f"{name} must be <= {maximum}, got {result}") + return result + + +def strict_integer_array(values, name: str, *, dtype: np.dtype, + ndim: int | None = None, + shape: tuple[int | None, ...] | None = None, + non_negative: bool = False) -> np.ndarray: + array = np.asarray(values) + if not np.issubdtype(array.dtype, np.integer) or np.issubdtype(array.dtype, np.bool_): + raise TypeError(f"{name} must contain integers, got dtype={array.dtype}") + if ndim is not None and array.ndim != ndim: + raise ValueError(f"{name} must have ndim={ndim}, got ndim={array.ndim}") + if shape is not None: + if array.ndim != len(shape) or any( + expected is not None and actual != expected + for actual, expected in zip(array.shape, shape) + ): + raise ValueError(f"{name} must have shape {shape}, got shape={array.shape}") + if non_negative and np.issubdtype(array.dtype, np.signedinteger): + if array.size and np.any(array < 0): + raise ValueError(f"{name} must contain non-negative integers") + + target = np.dtype(dtype) + # Fast path: a lossless (safe) cast is representable for every value by + # definition, so skip the O(n) min()/max() range scan. This keeps hot + # callers such as `objective.energy(labels)` (uint64->uint64, + # uint32->uint64) as cheap as a plain `np.asarray`. The range scan only + # runs for genuinely narrowing or sign-changing casts. + if array.size and not np.can_cast(array.dtype, target): + info = np.iinfo(target) + minimum = int(array.min()) + maximum = int(array.max()) + if minimum < info.min or maximum > info.max: + raise ValueError( + f"{name} values must fit dtype {target}, got range [{minimum}, {maximum}]" + ) + return np.ascontiguousarray(array, dtype=target) + + +def strict_offsets(offsets, ndim: int, name: str = "offsets") -> list[tuple[int, ...]]: + untyped = np.asarray(offsets) + if untyped.size == 0: + raise ValueError(f"{name} must not be empty") + array = strict_integer_array(offsets, name, dtype=np.int64, ndim=2) + if array.shape[1] != ndim: + raise ValueError( + f"{name}[0] must have length {ndim}; each offset must have " + f"length {ndim} to match the spatial ndim={ndim}, " + f"got shape={array.shape}" + ) + return [tuple(int(v) for v in row) for row in array] diff --git a/src/bioimage_cpp/affinities/compute_affinities.py b/src/bioimage_cpp/affinities/compute_affinities.py index 9d2c870..326df44 100644 --- a/src/bioimage_cpp/affinities/compute_affinities.py +++ b/src/bioimage_cpp/affinities/compute_affinities.py @@ -8,6 +8,7 @@ import numpy as np from .. import _core +from .._validation import strict_index, strict_offsets _COMPUTE_AFFINITIES_2D_BY_DTYPE = { @@ -95,9 +96,7 @@ def compute_affinities( f"labels must have one of dtypes ({supported}), got dtype={array.dtype}" ) from error - normalized_offsets = [ - [int(value) for value in offset] for offset in np.asarray(offsets).tolist() - ] + normalized_offsets = strict_offsets(offsets, array.ndim) if len(normalized_offsets) == 0: raise ValueError("offsets must not be empty") if any(len(offset) != array.ndim for offset in normalized_offsets): @@ -106,9 +105,7 @@ def compute_affinities( f"spatial ndim={array.ndim}" ) - n_threads = int(number_of_threads) - if n_threads < 1: - raise ValueError("number_of_threads must be >= 1") + n_threads = strict_index(number_of_threads, "number_of_threads", minimum=1) if ignore_label is None: typed_ignore: int | None = None diff --git a/src/bioimage_cpp/affinities/compute_embedding_distances.py b/src/bioimage_cpp/affinities/compute_embedding_distances.py index 9831ce9..f2303b8 100644 --- a/src/bioimage_cpp/affinities/compute_embedding_distances.py +++ b/src/bioimage_cpp/affinities/compute_embedding_distances.py @@ -7,6 +7,7 @@ import numpy as np from .. import _core +from .._validation import strict_index, strict_offsets _COMPUTE_EMBEDDING_DISTANCES_BY_NDIM = { @@ -67,9 +68,7 @@ def compute_embedding_distances( ) spatial_ndim = array.ndim - 1 - normalized_offsets = [ - [int(value) for value in offset] for offset in np.asarray(offsets).tolist() - ] + normalized_offsets = strict_offsets(offsets, spatial_ndim) if len(normalized_offsets) == 0: raise ValueError("offsets must not be empty") if any(len(offset) != spatial_ndim for offset in normalized_offsets): @@ -85,9 +84,7 @@ def compute_embedding_distances( f"norm must be one of ({supported}), got {norm!r}" ) - n_threads = int(number_of_threads) - if n_threads < 1: - raise ValueError("number_of_threads must be >= 1") + n_threads = strict_index(number_of_threads, "number_of_threads", minimum=1) run = _COMPUTE_EMBEDDING_DISTANCES_BY_NDIM[array.ndim] return run(array, normalized_offsets, norm_str, n_threads) diff --git a/src/bioimage_cpp/distance/_distance.py b/src/bioimage_cpp/distance/_distance.py index 9a42afe..36d0f5a 100644 --- a/src/bioimage_cpp/distance/_distance.py +++ b/src/bioimage_cpp/distance/_distance.py @@ -249,7 +249,7 @@ def non_maximum_distance_suppression( ``distance_map`` in NumPy axis order. Supported dtypes: ``int64``, ``uint64``, ``int32``, ``uint32``. number_of_threads - Reserved for future parallelization; currently single-threaded. + Number of worker threads used for pairwise candidate evaluation. Returns ------- @@ -259,9 +259,7 @@ def non_maximum_distance_suppression( Notes ----- - Uses an ``O(N^2)`` pairwise distance matrix internally; suitable for ``N`` - up to roughly 30k points. For larger candidate sets, threshold the - distance map more aggressively before calling. + Uses ``O(N^2)`` time and ``O(number_of_threads * N)`` auxiliary memory. """ function = "non_maximum_distance_suppression" diff --git a/src/bioimage_cpp/distance/_geodesic.py b/src/bioimage_cpp/distance/_geodesic.py index 3100d59..73eb4e3 100644 --- a/src/bioimage_cpp/distance/_geodesic.py +++ b/src/bioimage_cpp/distance/_geodesic.py @@ -5,15 +5,9 @@ - **masks** (regular grid): distances stay inside the nonzero region and never cross background voxels (fast-marching / Eikonal formulation, à la scikit-fmm). -- **surfaces** (triangle meshes): distances are measured across the mesh - surface (exact MMP geodesics, à la pygeodesic). - -.. note:: - - The C++ solvers are not implemented yet. These wrappers validate their - arguments fully and dispatch to the ``_core`` bindings, which currently - raise ``RuntimeError("... not yet implemented")``. The reference behaviour - lives in ``development/distance/``. +- **surfaces** (triangle meshes): first-order Kimmel--Sethian fast marching + across the mesh surface. This approximates exact MMP geodesics; error grows + on very obtuse triangulations because obtuse-angle unfolding is not used. """ from __future__ import annotations @@ -23,12 +17,13 @@ import numpy as np from .. import _core +from .._validation import strict_integer_array from ._distance import _as_binary_input, _normalize_sampling, _normalize_threads def _as_coordinates(points: np.ndarray, ndim: int, function: str, name: str) -> np.ndarray: """Coerce a point set to a C-contiguous ``(n, ndim)`` int64 array.""" - array = np.ascontiguousarray(points) + array = np.asarray(points) if array.ndim == 1 and array.size == ndim: # A single coordinate may be given as a flat (ndim,) vector. array = array.reshape(1, ndim) @@ -41,17 +36,17 @@ def _as_coordinates(points: np.ndarray, ndim: int, function: str, name: str) -> f"{function}: {name}.shape[1] must equal the mask ndim ({ndim}), " f"got {name}.shape[1]={array.shape[1]}" ) - return np.ascontiguousarray(array, dtype=np.int64) + return strict_integer_array(array, name, dtype=np.int64, ndim=2) def _as_vertex_indices(indices: np.ndarray, function: str, name: str) -> np.ndarray: """Coerce vertex indices to a C-contiguous 1-D int64 array.""" - array = np.atleast_1d(np.ascontiguousarray(indices)) + array = np.atleast_1d(np.asarray(indices)) if array.ndim != 1: raise ValueError( f"{function}: {name} must be 1-D vertex indices, got ndim={array.ndim}" ) - return np.ascontiguousarray(array, dtype=np.int64) + return strict_integer_array(array, name, dtype=np.int64, ndim=1) def _as_mesh(vertices: np.ndarray, faces: np.ndarray, function: str): @@ -62,7 +57,9 @@ def _as_mesh(vertices: np.ndarray, faces: np.ndarray, function: str): f"{function}: vertices must have shape (n_vertices, 3), " f"got shape={vertices.shape}" ) - faces = np.ascontiguousarray(faces, dtype=np.int64) + if not np.all(np.isfinite(vertices)): + raise ValueError(f"{function}: vertices must contain only finite values") + faces = strict_integer_array(faces, "faces", dtype=np.int64, ndim=2) if faces.ndim != 2 or faces.shape[1] != 3: raise ValueError( f"{function}: faces must have shape (n_faces, 3), got shape={faces.shape}" @@ -82,9 +79,23 @@ def _as_speed( f"{function}: speed must have shape {tuple(expected_shape)}, " f"got shape={array.shape}" ) + if not np.all(np.isfinite(array)) or np.any(array <= 0.0): + raise ValueError(f"{function}: speed values must be finite and strictly positive") return array +def _require_foreground_points( + mask: np.ndarray, points: np.ndarray, function: str, name: str +) -> None: + for axis, extent in enumerate(mask.shape): + if points.size and ( + np.any(points[:, axis] < 0) or np.any(points[:, axis] >= extent) + ): + raise ValueError(f"{function}: {name} contains an out of bounds coordinate") + if points.size and np.any(mask[tuple(points.T)] == 0): + raise ValueError(f"{function}: {name} must lie inside the foreground mask") + + def geodesic_distance_field( mask: np.ndarray, sources: np.ndarray, @@ -119,8 +130,7 @@ def geodesic_distance_field( If ``True``, also return the per-axis gradient of the field (see Returns). Analogous to :func:`vector_difference_transform`. number_of_threads - ``0`` uses ``hardware_concurrency``; a positive value pins the thread - count. Default ``1``. + Retained for API consistency; a single field solve is serial. Returns ------- @@ -139,6 +149,7 @@ def geodesic_distance_field( binary = _as_binary_input(mask, function) ndim = binary.ndim sources_arr = _as_coordinates(sources, ndim, function, "sources") + _require_foreground_points(binary, sources_arr, function, "sources") sampling_values = _normalize_sampling(sampling, ndim, function) speed_arr = _as_speed(speed, tuple(binary.shape), function) n_threads = _normalize_threads(number_of_threads, function) @@ -215,6 +226,7 @@ def geodesic_distances( binary = _as_binary_input(mask, function) ndim = binary.ndim points_arr = _as_coordinates(points, ndim, function, "points") + _require_foreground_points(binary, points_arr, function, "points") sampling_values = _normalize_sampling(sampling, ndim, function) speed_arr = _as_speed(speed, tuple(binary.shape), function) n_threads = _normalize_threads(number_of_threads, function) @@ -247,8 +259,7 @@ def geodesic_distance_field_mesh( Optional per-vertex speed of shape ``(n_vertices,)``. ``None`` gives unit-speed geodesic distance. number_of_threads - ``0`` uses ``hardware_concurrency``; a positive value pins the thread - count. Default ``1``. + Retained for API consistency; a single field solve is serial. Returns ------- diff --git a/src/bioimage_cpp/filters/_filters.py b/src/bioimage_cpp/filters/_filters.py index 1e1b56e..813660b 100644 --- a/src/bioimage_cpp/filters/_filters.py +++ b/src/bioimage_cpp/filters/_filters.py @@ -23,6 +23,7 @@ import numpy as np from .. import _core +from .._validation import strict_index _FLOAT_INPUT_DTYPES = (np.float32, np.float64) _INT_INPUT_DTYPES = (np.uint8, np.uint16) @@ -38,15 +39,25 @@ def _prepare_input(image: np.ndarray, function: str) -> tuple[np.ndarray, np.dty f"{function}: image must be 2D or 3D, got ndim={image.ndim}" ) if image.dtype == np.float32: - return np.ascontiguousarray(image), np.dtype(np.float32) - if image.dtype == np.float64: - return np.ascontiguousarray(image, dtype=np.float32), np.dtype(np.float64) - if image.dtype in (np.dtype(t) for t in _INT_INPUT_DTYPES): + prepared, out_dtype = np.ascontiguousarray(image), np.dtype(np.float32) + elif image.dtype == np.float64: + prepared, out_dtype = np.ascontiguousarray(image, dtype=np.float32), np.dtype(np.float64) + elif image.dtype in (np.dtype(t) for t in _INT_INPUT_DTYPES): return np.ascontiguousarray(image, dtype=np.float32), np.dtype(np.float32) - raise TypeError( - f"{function}: image dtype must be one of (float32, float64, uint8, " - f"uint16), got dtype={image.dtype}" - ) + else: + raise TypeError( + f"{function}: image dtype must be one of (float32, float64, uint8, " + f"uint16), got dtype={image.dtype}" + ) + # Reject non-finite float inputs without allocating a whole-array boolean + # temporary: min()/max() are single-pass reductions and propagate NaN/inf + # (any NaN -> NaN, any +/-inf -> +/-inf), so a non-finite extremum flags a + # non-finite sample. Integer inputs return above and skip this entirely. + if prepared.size and not ( + np.isfinite(prepared.min()) and np.isfinite(prepared.max()) + ): + raise ValueError(f"{function}: image must contain only finite values") + return prepared, out_dtype def _broadcast_per_axis( @@ -56,13 +67,16 @@ def _broadcast_per_axis( function: str, ) -> tuple[float, ...]: if np.isscalar(value): - return (float(value),) * ndim - seq = tuple(float(v) for v in value) + seq = (float(value),) * ndim + else: + seq = tuple(float(v) for v in value) if len(seq) != ndim: raise ValueError( f"{function}: {name} must be a scalar or a sequence of length " f"{ndim}, got length {len(seq)}" ) + if not all(np.isfinite(v) and v > 0.0 for v in seq): + raise ValueError(f"{function}: {name} values must be finite and positive") return seq @@ -72,8 +86,8 @@ def _broadcast_order( function: str, ) -> tuple[int, ...]: if np.isscalar(value): - return (int(value),) * ndim - seq = tuple(int(v) for v in value) + return (strict_index(value, "order", minimum=0, maximum=2),) * ndim + seq = tuple(strict_index(v, "order", minimum=0, maximum=2) for v in value) if len(seq) != ndim: raise ValueError( f"{function}: order must be a scalar or a sequence of length " @@ -82,6 +96,13 @@ def _broadcast_order( return seq +def _normalize_window(window_size: float, function: str) -> float: + value = float(window_size) + if not np.isfinite(value) or value < 0.0: + raise ValueError(f"{function}: window_size must be finite and non-negative") + return value + + def _finalise(result: np.ndarray, out_dtype: np.dtype) -> np.ndarray: if result.dtype == out_dtype: return result @@ -98,13 +119,14 @@ def gaussian_smoothing( function = "gaussian_smoothing" prepared, out_dtype = _prepare_input(image, function) sigmas = _broadcast_per_axis(sigma, prepared.ndim, "sigma", function) + window = _normalize_window(window_size, function) if prepared.ndim == 2: result = _core._gaussian_smoothing_2d_float32( - prepared, sigmas[0], sigmas[1], float(window_size) + prepared, sigmas[0], sigmas[1], window ) else: result = _core._gaussian_smoothing_3d_float32( - prepared, sigmas[0], sigmas[1], sigmas[2], float(window_size) + prepared, sigmas[0], sigmas[1], sigmas[2], window ) return _finalise(result, out_dtype) @@ -121,15 +143,16 @@ def gaussian_derivative( prepared, out_dtype = _prepare_input(image, function) sigmas = _broadcast_per_axis(sigma, prepared.ndim, "sigma", function) orders = _broadcast_order(order, prepared.ndim, function) + window = _normalize_window(window_size, function) if prepared.ndim == 2: result = _core._gaussian_derivative_2d_float32( prepared, sigmas[0], sigmas[1], - orders[0], orders[1], float(window_size), + orders[0], orders[1], window, ) else: result = _core._gaussian_derivative_3d_float32( prepared, sigmas[0], sigmas[1], sigmas[2], - orders[0], orders[1], orders[2], float(window_size), + orders[0], orders[1], orders[2], window, ) return _finalise(result, out_dtype) @@ -144,13 +167,14 @@ def gaussian_gradient_magnitude( function = "gaussian_gradient_magnitude" prepared, out_dtype = _prepare_input(image, function) sigmas = _broadcast_per_axis(sigma, prepared.ndim, "sigma", function) + window = _normalize_window(window_size, function) if prepared.ndim == 2: result = _core._gaussian_gradient_magnitude_2d_float32( - prepared, sigmas[0], sigmas[1], float(window_size) + prepared, sigmas[0], sigmas[1], window ) else: result = _core._gaussian_gradient_magnitude_3d_float32( - prepared, sigmas[0], sigmas[1], sigmas[2], float(window_size) + prepared, sigmas[0], sigmas[1], sigmas[2], window ) return _finalise(result, out_dtype) @@ -166,13 +190,14 @@ def laplacian_of_gaussian( function = "laplacian_of_gaussian" prepared, out_dtype = _prepare_input(image, function) sigmas = _broadcast_per_axis(sigma, prepared.ndim, "sigma", function) + window = _normalize_window(window_size, function) if prepared.ndim == 2: result = _core._laplacian_of_gaussian_2d_float32( - prepared, sigmas[0], sigmas[1], float(window_size) + prepared, sigmas[0], sigmas[1], window ) else: result = _core._laplacian_of_gaussian_3d_float32( - prepared, sigmas[0], sigmas[1], sigmas[2], float(window_size) + prepared, sigmas[0], sigmas[1], sigmas[2], window ) return _finalise(result, out_dtype) @@ -192,13 +217,14 @@ def hessian_of_gaussian_eigenvalues( function = "hessian_of_gaussian_eigenvalues" prepared, out_dtype = _prepare_input(image, function) sigmas = _broadcast_per_axis(sigma, prepared.ndim, "sigma", function) + window = _normalize_window(window_size, function) if prepared.ndim == 2: result = _core._hessian_of_gaussian_eigenvalues_2d_float32( - prepared, sigmas[0], sigmas[1], float(window_size) + prepared, sigmas[0], sigmas[1], window ) else: result = _core._hessian_of_gaussian_eigenvalues_3d_float32( - prepared, sigmas[0], sigmas[1], sigmas[2], float(window_size) + prepared, sigmas[0], sigmas[1], sigmas[2], window ) return _finalise(result, out_dtype) @@ -219,18 +245,19 @@ def structure_tensor_eigenvalues( prepared, out_dtype = _prepare_input(image, function) inner = _broadcast_per_axis(inner_sigma, prepared.ndim, "inner_sigma", function) outer = _broadcast_per_axis(outer_sigma, prepared.ndim, "outer_sigma", function) + window = _normalize_window(window_size, function) if prepared.ndim == 2: result = _core._structure_tensor_eigenvalues_2d_float32( prepared, inner[0], inner[1], outer[0], outer[1], - float(window_size), + window, ) else: result = _core._structure_tensor_eigenvalues_3d_float32( prepared, inner[0], inner[1], inner[2], outer[0], outer[1], outer[2], - float(window_size), + window, ) return _finalise(result, out_dtype) diff --git a/src/bioimage_cpp/flow/_flow.py b/src/bioimage_cpp/flow/_flow.py index 43287c0..7cd9989 100644 --- a/src/bioimage_cpp/flow/_flow.py +++ b/src/bioimage_cpp/flow/_flow.py @@ -7,6 +7,7 @@ import numpy as np from .. import _core +from .._validation import strict_index def _normalize_mask(fg_mask: np.ndarray, shape: tuple[int, ...]) -> np.ndarray: @@ -34,8 +35,8 @@ def _normalize_sigma( raise ValueError( f"spacing must be a sequence of length {ndim}, got shape {sp.shape}" ) - if np.any(sp <= 0): - raise ValueError("spacing values must be positive") + if not np.all(np.isfinite(sp)) or np.any(sp <= 0): + raise ValueError("spacing values must be finite and positive") return tuple((sig / sp).tolist()) return sigma @@ -110,9 +111,7 @@ def compute_flow_density( f"flow first axis must match spatial ndim={ndim}, got {array.shape[0]}" ) - n_steps = int(n_iter) - if n_steps < 0: - raise ValueError("n_iter must be >= 0") + n_steps = strict_index(n_iter, "n_iter", minimum=0) step_size = float(dt) if not np.isfinite(step_size) or step_size < 0: raise ValueError("dt must be finite and >= 0") @@ -122,10 +121,10 @@ def compute_flow_density( integration_method = str(method) if integration_method not in ("euler", "rk2"): raise ValueError(f"method must be 'euler' or 'rk2', got '{integration_method}'") - n_threads = int(number_of_threads) - if n_threads < 1: - raise ValueError("number_of_threads must be >= 1") + n_threads = strict_index(number_of_threads, "number_of_threads", minimum=1) + # Finiteness of the flow tensor is validated once in the binding layer + # (bindings/flow.cxx); it is deliberately not re-scanned here. contiguous_flow = np.ascontiguousarray(array, dtype=np.float32) mask = _normalize_mask(fg_mask, tuple(contiguous_flow.shape[1:])) diff --git a/src/bioimage_cpp/graph/__init__.py b/src/bioimage_cpp/graph/__init__.py index 3d01e55..3b4fe0a 100644 --- a/src/bioimage_cpp/graph/__init__.py +++ b/src/bioimage_cpp/graph/__init__.py @@ -43,6 +43,7 @@ import numpy as np from .. import _core +from .._validation import strict_integer_array, strict_offsets from ._shared import ( _as_1d_array, _as_coordinate_array, @@ -110,19 +111,12 @@ def deserialize(cls, serialization): def _normalize_projection_offsets(offsets, ndim: int) -> list[list[int]]: - normalized: list[list[int]] = [] - for index, offset in enumerate(offsets): - values = [int(v) for v in offset] - if len(values) != ndim: - raise ValueError( - f"offsets[{index}] must have length {ndim}, got length={len(values)}" - ) - normalized.append(values) - return normalized + return [list(offset) for offset in strict_offsets(offsets, ndim)] def _normalize_projection_strides(strides, ndim: int) -> list[int]: - values = [int(v) for v in strides] + values_array = strict_integer_array(strides, "strides", dtype=np.uint64, ndim=1) + values = [int(v) for v in values_array] if len(values) != ndim: raise ValueError( f"strides must have length {ndim}, got length={len(values)}" diff --git a/src/bioimage_cpp/graph/_shared.py b/src/bioimage_cpp/graph/_shared.py index 7c5ac39..2ece8e7 100644 --- a/src/bioimage_cpp/graph/_shared.py +++ b/src/bioimage_cpp/graph/_shared.py @@ -11,6 +11,7 @@ import numpy as np from .. import _core +from .._validation import strict_index, strict_integer_array, strict_offsets _REGION_ADJACENCY_GRAPH_BY_DTYPE = { @@ -25,46 +26,48 @@ def _as_shape(shape, ndim: int, name: str = "shape") -> list[int]: - array = np.asarray(shape) + array = strict_integer_array(shape, name, dtype=np.uint64, ndim=1, non_negative=True) 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") if np.any(array <= 0): raise ValueError(f"{name} dimensions must be greater than zero") return [int(axis_size) for axis_size in array] def _as_coordinate_array(coordinate, ndim: int, name: str) -> np.ndarray: - array = np.asarray(coordinate, dtype=np.uint64) + array = strict_integer_array( + coordinate, name, dtype=np.uint64, ndim=1, non_negative=True + ) if array.ndim != 1 or array.shape[0] != ndim: raise ValueError(f"{name} must be a 1D sequence of length {ndim}") return np.ascontiguousarray(array) def _as_offset_array(offset, ndim: int, name: str) -> np.ndarray: - array = np.asarray(offset, dtype=np.int64) + array = strict_integer_array(offset, name, dtype=np.int64, ndim=1) if array.ndim != 1 or array.shape[0] != ndim: raise ValueError(f"{name} must be a 1D sequence of length {ndim}") return np.ascontiguousarray(array) def _as_uv_array(uvs, name: str) -> np.ndarray: - array = np.asarray(uvs, dtype=np.uint64) + array = strict_integer_array(uvs, name, dtype=np.uint64, non_negative=True) if array.ndim != 2 or array.shape[1] != 2: raise ValueError(f"{name} must have shape (n_edges, 2)") return np.ascontiguousarray(array) def _as_node_array(nodes, name: str) -> np.ndarray: - array = np.asarray(nodes, dtype=np.uint64) + array = strict_integer_array(nodes, name, dtype=np.uint64, ndim=1, non_negative=True) if array.ndim != 1: raise ValueError(f"{name} must be a 1D array") return np.ascontiguousarray(array) def _as_serialization_array(serialization) -> np.ndarray: - array = np.asarray(serialization, dtype=np.uint64) + array = strict_integer_array( + serialization, "serialization", dtype=np.uint64, ndim=1, non_negative=True + ) if array.ndim != 1: raise ValueError("serialization must be a 1D array") if array.size < 2: @@ -76,16 +79,11 @@ def _as_serialization_array(serialization) -> np.ndarray: def _copy_graph(graph) -> _core.UndirectedGraph: - # `uv_ids()` always returns a unique list (graphs deduplicate on insert), - # so we can use the bulk constructor that skips per-edge hash dedup — - # significantly faster than `insert_edges` for large graphs. The result - # is a ``_core.UndirectedGraph``; downstream code (objectives, solvers, - # validators) uses base-class methods that work identically. - if graph.number_of_edges == 0: - return _core.UndirectedGraph(int(graph.number_of_nodes)) - return _core.UndirectedGraph.from_unique_edges( - int(graph.number_of_nodes), graph.uv_ids() - ) + # `clone()` preserves edge ids and accepts any valid insertion order, and + # defers CSR construction. The public `from_unique_edges` fast path is not + # usable here: it intentionally requires sorted, canonical, duplicate-free + # input, which `uv_ids()` (insertion order) does not guarantee. + return graph.clone() def _as_edge_costs(edge_costs, graph) -> np.ndarray: @@ -98,7 +96,9 @@ def _as_edge_costs(edge_costs, graph) -> np.ndarray: def _as_node_labels(labels, graph) -> np.ndarray: - array = np.asarray(labels, dtype=np.uint64) + array = strict_integer_array( + labels, "labels", dtype=np.uint64, ndim=1, non_negative=True + ) if array.ndim != 1: raise ValueError("labels must be a 1D array") if array.shape[0] != graph.number_of_nodes: @@ -118,7 +118,9 @@ def _as_1d_array(values, dtype, name: str, expected_size: int) -> np.ndarray: def _dense_labels(labels) -> np.ndarray: - labels = np.asarray(labels, dtype=np.uint64) + labels = strict_integer_array( + labels, "labels", dtype=np.uint64, non_negative=True + ) _, dense = np.unique(labels, return_inverse=True) return np.ascontiguousarray(dense.astype(np.uint64, copy=False)) @@ -154,10 +156,7 @@ def _normalize_labels(labels: np.ndarray) -> np.ndarray: def _normalize_number_of_threads(number_of_threads: int) -> int: - number_of_threads = int(number_of_threads) - if number_of_threads < 0: - raise ValueError("number_of_threads must be non-negative") - return number_of_threads + return strict_index(number_of_threads, "number_of_threads", minimum=0) def _resolve_weight_dtype(array, name: str) -> np.ndarray: @@ -216,7 +215,7 @@ def _as_grid_data(values, graph, name: str, *, with_channels: bool) -> np.ndarra def _normalize_grid_offsets(offsets, ndim: int, n_channels: int) -> list[tuple[int, ...]]: - normalized = [tuple(int(value) for value in offset) for offset in offsets] + normalized = strict_offsets(offsets, ndim) if len(normalized) != n_channels: raise ValueError( "offsets length must match affinities channel count, got " diff --git a/src/bioimage_cpp/graph/features/__init__.py b/src/bioimage_cpp/graph/features/__init__.py index 9e4e0b6..99bc7c4 100644 --- a/src/bioimage_cpp/graph/features/__init__.py +++ b/src/bioimage_cpp/graph/features/__init__.py @@ -22,6 +22,7 @@ import numpy as np from .. import _core +from ..._validation import strict_offsets from .._shared import ( _as_grid_data, _as_uv_array, @@ -244,7 +245,9 @@ def lifted_edges_from_affinities( f"rag shape={tuple(rag.shape)}, labels shape={label_array.shape}" ) - normalized_offsets = [tuple(int(value) for value in offset) for offset in offsets] + if np.asarray(offsets).size == 0: + return np.empty((0, 2), dtype=np.uint64) + normalized_offsets = strict_offsets(offsets, label_array.ndim) if any(len(offset) != label_array.ndim for offset in normalized_offsets): raise ValueError("each offset must have length matching labels ndim") @@ -516,7 +519,7 @@ def _accumulate_affinity_features( f"affinities shape={affinity_array.shape}, labels shape={label_array.shape}" ) - normalized_offsets = [tuple(int(value) for value in offset) for offset in offsets] + normalized_offsets = strict_offsets(offsets, label_array.ndim) if len(normalized_offsets) != affinity_array.shape[0]: raise ValueError( "offsets length must match affinities channel count, got " @@ -555,7 +558,7 @@ def _accumulate_lifted_affinity_features( f"affinities shape={affinity_array.shape}, labels shape={label_array.shape}" ) - normalized_offsets = [tuple(int(value) for value in offset) for offset in offsets] + normalized_offsets = strict_offsets(offsets, label_array.ndim) if len(normalized_offsets) != affinity_array.shape[0]: raise ValueError( "offsets length must match affinities channel count, got " diff --git a/src/bioimage_cpp/graph/lifted_multicut/__init__.py b/src/bioimage_cpp/graph/lifted_multicut/__init__.py index e7d30df..c5a4dcb 100644 --- a/src/bioimage_cpp/graph/lifted_multicut/__init__.py +++ b/src/bioimage_cpp/graph/lifted_multicut/__init__.py @@ -186,14 +186,10 @@ def __init__( base_graph = graph base_costs = _as_edge_costs(edge_costs, base_graph) - # Use the bulk constructor for the lifted graph's base portion to - # bypass the per-edge hash dedup that ``insert_edges`` performs. - if int(base_graph.number_of_edges) > 0: - lifted_graph = _core.UndirectedGraph.from_unique_edges( - int(base_graph.number_of_nodes), base_graph.uv_ids() - ) - else: - lifted_graph = _core.UndirectedGraph(int(base_graph.number_of_nodes)) + # Preserve the base graph's edge IDs: base costs are indexed in that + # order, whereas the sorted bulk constructor is only valid when its + # input already has lexicographic edge order. + lifted_graph = base_graph.clone() weights_list = [base_costs.copy()] diff --git a/src/bioimage_cpp/graph/mutex_watershed.py b/src/bioimage_cpp/graph/mutex_watershed.py index 69d490a..61138a5 100644 --- a/src/bioimage_cpp/graph/mutex_watershed.py +++ b/src/bioimage_cpp/graph/mutex_watershed.py @@ -183,6 +183,10 @@ def semantic_mutex_watershed_clustering( int(mutex_uv_array.shape[0]), ) semantic_uv_array = _as_uv_array(semantic_node_classes, "semantic_node_classes") + if semantic_uv_array.size and np.any( + semantic_uv_array[:, 1] > np.iinfo(np.int64).max + ): + raise ValueError("semantic class ids must be <= np.iinfo(np.int64).max") semantic_cost_array = _as_1d_array( semantic_cost_array, semantic_cost_array.dtype, diff --git a/src/bioimage_cpp/segmentation/mutex_watershed.py b/src/bioimage_cpp/segmentation/mutex_watershed.py index 1927a84..b8a29fa 100644 --- a/src/bioimage_cpp/segmentation/mutex_watershed.py +++ b/src/bioimage_cpp/segmentation/mutex_watershed.py @@ -7,6 +7,7 @@ import numpy as np from .. import _core +from .._validation import strict_index, strict_offsets _MUTEX_WATERSHED_BY_DTYPE = { np.dtype("float32"): _core._mutex_watershed_grid_float32, @@ -80,8 +81,8 @@ def mutex_watershed( f"affinities must have one of dtypes ({supported}), got dtype={dtype}" ) from error - normalized_offsets = [tuple(int(value) for value in offset) for offset in offsets] spatial_ndim = array.ndim - 1 + normalized_offsets = strict_offsets(offsets, spatial_ndim) if len(normalized_offsets) != array.shape[0]: raise ValueError( "offsets length must match affinities channel count, got " @@ -92,8 +93,9 @@ def mutex_watershed( "each offset must have length matching the spatial ndim, got " f"spatial ndim={spatial_ndim}" ) - if number_of_attractive_channels < 0: - raise ValueError("number_of_attractive_channels must be non-negative") + number_of_attractive_channels = strict_index( + number_of_attractive_channels, "number_of_attractive_channels", minimum=0 + ) if number_of_attractive_channels > array.shape[0]: raise ValueError("number_of_attractive_channels must be <= number of channels") @@ -199,10 +201,10 @@ def semantic_mutex_watershed( f"affinities must have one of dtypes ({supported}), got dtype={dtype}" ) from error - normalized_offsets = [tuple(int(value) for value in offset) for offset in offsets] - number_of_offsets = len(normalized_offsets) number_of_channels = int(array.shape[0]) spatial_ndim = array.ndim - 1 + normalized_offsets = strict_offsets(offsets, spatial_ndim) + number_of_offsets = len(normalized_offsets) if number_of_offsets == 0: raise ValueError("offsets must not be empty") if number_of_channels <= number_of_offsets: @@ -216,8 +218,9 @@ def semantic_mutex_watershed( "each offset must have length matching the spatial ndim, got " f"spatial ndim={spatial_ndim}" ) - if number_of_attractive_channels < 0: - raise ValueError("number_of_attractive_channels must be non-negative") + number_of_attractive_channels = strict_index( + number_of_attractive_channels, "number_of_attractive_channels", minimum=0 + ) if number_of_attractive_channels > number_of_offsets: raise ValueError( "number_of_attractive_channels must be <= len(offsets)" @@ -259,7 +262,9 @@ def _normalize_strides( raise ValueError("randomized_strides requires strides") return None - normalized = tuple(int(stride) for stride in strides) + normalized = tuple( + strict_index(stride, "stride", minimum=1) for stride in strides + ) if len(normalized) != spatial_ndim: raise ValueError( "strides length must match the spatial ndim, got " diff --git a/src/bioimage_cpp/segmentation/watershed.py b/src/bioimage_cpp/segmentation/watershed.py index 08d3a6b..80582dc 100644 --- a/src/bioimage_cpp/segmentation/watershed.py +++ b/src/bioimage_cpp/segmentation/watershed.py @@ -7,6 +7,7 @@ import numpy as np from .. import _core +from .._validation import strict_offsets _WATERSHED_BY_DTYPE: dict[np.dtype, dict[np.dtype, object]] = { np.dtype("float32"): { @@ -47,10 +48,11 @@ def watershed( ) -> np.ndarray: """Run a marker-controlled watershed on a 2D or 3D heightmap. - Pixels with non-zero ``markers`` values are treated as seeds. Their - label is propagated to neighbouring pixels in order of increasing - ``image`` value, using axis-aligned connectivity (4-neighbours in 2D, - 6-neighbours in 3D). Tie-breaking on equal heights is unspecified. + Pixels with non-zero ``markers`` values are treated as seeds. Heights are + quantized to 65,536 levels and propagated with Meyer-style monotone + flooding using axis-aligned connectivity. Ordering within one quantization + level is unspecified, including for distinct input values that quantize to + the same level. Parameters ---------- @@ -201,7 +203,7 @@ def watershed_from_affinities( f"markers must have one of dtypes ({supported}), got dtype={markers_array.dtype}" ) from error - normalized_offsets = [tuple(int(value) for value in offset) for offset in offsets] + normalized_offsets = strict_offsets(offsets, spatial_ndim) if len(normalized_offsets) != n_channels: raise ValueError( "offsets count must equal affinities channel count, got " diff --git a/src/bioimage_cpp/transformation/_transformation.py b/src/bioimage_cpp/transformation/_transformation.py index 205f77c..22b7bf4 100644 --- a/src/bioimage_cpp/transformation/_transformation.py +++ b/src/bioimage_cpp/transformation/_transformation.py @@ -364,6 +364,14 @@ def affine_transform( contiguous = np.ascontiguousarray(array) typed_fill = _normalize_fill_value(fill_value, dtype) output = _validate_or_allocate_out(out, output_shape, dtype) + # Guard against `out` aliasing an input: the sampling loop reads input + # samples that later output pixels still need. + if np.shares_memory(contiguous, output): + contiguous = contiguous.copy() + if np.shares_memory(normalized_matrix, output): + normalized_matrix = normalized_matrix.copy() + if np.shares_memory(starts, output): + starts = starts.copy() run(contiguous, normalized_matrix, starts, output, order, typed_fill) return output @@ -444,5 +452,10 @@ def map_coordinates( contiguous = np.ascontiguousarray(array) typed_fill = _normalize_fill_value(fill_value, dtype) output = _validate_or_allocate_out(out, output_shape, dtype) + # Guard against `out` aliasing an input (see affine_transform). + if np.shares_memory(contiguous, output): + contiguous = contiguous.copy() + if np.shares_memory(coords, output): + coords = coords.copy() run(contiguous, coords, output, order, typed_fill) return output diff --git a/tests/affinities/test_compute_affinities.py b/tests/affinities/test_compute_affinities.py index 6e25db5..7123091 100644 --- a/tests/affinities/test_compute_affinities.py +++ b/tests/affinities/test_compute_affinities.py @@ -4,6 +4,22 @@ import bioimage_cpp as bic +def test_minimum_int64_offset_is_safely_out_of_bounds(): + labels = np.ones((2, 2), dtype=np.uint32) + affinities, valid = bic.affinities.compute_affinities( + labels, [[np.iinfo(np.int64).min, 0]], number_of_threads=1 + ) + assert not affinities.any() + assert not valid.any() + + +def test_fractional_offsets_are_rejected(): + with pytest.raises(TypeError, match="integers"): + bic.affinities.compute_affinities( + np.ones((2, 2), dtype=np.uint32), [[0.5, 0]] + ) + + def _numpy_reference(labels, offsets, ignore_label=None): """Slow but obvious reference: nested Python loops over voxels.""" labels = np.asarray(labels) diff --git a/tests/affinities/test_compute_embedding_distances.py b/tests/affinities/test_compute_embedding_distances.py index 9427ae0..7b6beda 100644 --- a/tests/affinities/test_compute_embedding_distances.py +++ b/tests/affinities/test_compute_embedding_distances.py @@ -4,6 +4,14 @@ import bioimage_cpp as bic +def test_minimum_int64_offset_is_safely_out_of_bounds(): + values = np.ones((1, 2, 2), dtype=np.float32) + distances = bic.affinities.compute_embedding_distances( + values, [[np.iinfo(np.int64).min, 0]], number_of_threads=1 + ) + assert not distances.any() + + def _numpy_reference(values, offsets, norm): """Slow but obvious reference: nested Python loops over voxels.""" values = np.asarray(values, dtype=np.float64) diff --git a/tests/distance/test_geodesic_distance.py b/tests/distance/test_geodesic_distance.py index 3441212..f3a0daa 100644 --- a/tests/distance/test_geodesic_distance.py +++ b/tests/distance/test_geodesic_distance.py @@ -206,6 +206,29 @@ def test_mask_invalid_arguments(): bic.distance.geodesic_distance_field(mask, np.zeros((1, 3), np.int64)) +def test_mask_rejects_background_sources_and_invalid_speed(): + mask = np.ones((4, 4), np.uint8) + mask[1, 1] = 0 + with pytest.raises(ValueError, match="foreground"): + bic.distance.geodesic_distance_field(mask, np.array([[1, 1]], np.int64)) + for value in (0.0, -1.0, np.nan, np.inf): + speed = np.ones(mask.shape, np.float64) + speed[0, 0] = value + with pytest.raises(ValueError, match="finite and strictly positive"): + bic.distance.geodesic_distance_field( + mask, np.array([[0, 0]], np.int64), speed=speed + ) + + +def test_mesh_rejects_non_finite_vertices(): + verts, faces, _ = flat_grid_mesh(3) + verts[0, 0] = np.nan + with pytest.raises(ValueError, match="finite"): + bic.distance.geodesic_distance_field_mesh( + verts, faces, np.array([0], np.int64) + ) + + # --------------------------------------------------------------------------- # # mask: gradient of the field # --------------------------------------------------------------------------- # diff --git a/tests/distance/test_non_maximum_distance_suppression.py b/tests/distance/test_non_maximum_distance_suppression.py index 194e6d3..00390cd 100644 --- a/tests/distance/test_non_maximum_distance_suppression.py +++ b/tests/distance/test_non_maximum_distance_suppression.py @@ -128,6 +128,15 @@ def test_deterministic(): assert np.array_equal(a, b) +def test_thread_counts_are_identical(): + rng = np.random.default_rng(17) + dm = rng.uniform(0, 5, size=(40, 40)).astype(np.float32) + points = rng.integers(0, 40, size=(300, 2), dtype=np.int64) + single = nms(dm, points, number_of_threads=1) + multi = nms(dm, points, number_of_threads=4) + np.testing.assert_array_equal(single, multi) + + def test_invalid_points_ndim(): dm = np.ones((10, 10), dtype=np.float32) with pytest.raises(ValueError): diff --git a/tests/graph/test_rag_coordinates.py b/tests/graph/test_rag_coordinates.py index 5128c88..fe61a2b 100644 --- a/tests/graph/test_rag_coordinates.py +++ b/tests/graph/test_rag_coordinates.py @@ -4,6 +4,15 @@ import bioimage_cpp as bic +@pytest.mark.parametrize("number_of_threads", [1, 4]) +def test_negative_labels_raise_without_terminating(number_of_threads): + labels = np.zeros((8, 8), dtype=np.int32) + labels[3, 3] = -1 + rag = bic.graph.region_adjacency_graph(np.zeros((8, 8), dtype=np.uint32)) + with pytest.raises(ValueError, match="negative"): + bic.graph.rag_coordinates(rag, labels, number_of_threads=number_of_threads) + + def _labels_2d(): return np.array( [ diff --git a/tests/graph/test_rag_features.py b/tests/graph/test_rag_features.py index 15aeb08..77a2bd9 100644 --- a/tests/graph/test_rag_features.py +++ b/tests/graph/test_rag_features.py @@ -202,3 +202,14 @@ def test_rag_feature_validation(): bic.graph.features.affinity_features( rag, labels, np.ones((2, 1, 2)), offsets=[[0, 1]] ) + + +@pytest.mark.parametrize("number_of_threads", [1, 4]) +def test_negative_labels_raise_without_terminating(number_of_threads): + labels = np.zeros((8, 8), dtype=np.int32) + labels[3, 3] = -1 + rag = bic.graph.region_adjacency_graph(np.zeros((8, 8), dtype=np.uint32)) + with pytest.raises(ValueError, match="negative"): + bic.graph.features.edge_map_features( + rag, labels, np.zeros((8, 8)), number_of_threads=number_of_threads + ) diff --git a/tests/graph/test_semantic_mutex_watershed_graph.py b/tests/graph/test_semantic_mutex_watershed_graph.py index 8afdd19..c2a7896 100644 --- a/tests/graph/test_semantic_mutex_watershed_graph.py +++ b/tests/graph/test_semantic_mutex_watershed_graph.py @@ -274,3 +274,16 @@ def test_out_of_range_semantic_node_raises(): graph, edge_costs, mutex_uvs, mutex_costs, semantic_node_classes, semantic_costs, ) + + +def test_semantic_class_above_int64_range_raises(): + graph = bic.graph.UndirectedGraph(1) + with pytest.raises(ValueError, match="semantic class ids"): + bic.graph.mutex_watershed.semantic_mutex_watershed_clustering( + graph, + np.empty(0, np.float64), + np.empty((0, 2), np.uint64), + np.empty(0, np.float64), + np.array([[0, 2**63]], np.uint64), + np.array([1.0], np.float64), + ) diff --git a/tests/graph/test_undirected_graph.py b/tests/graph/test_undirected_graph.py index f9198ac..bdd489a 100644 --- a/tests/graph/test_undirected_graph.py +++ b/tests/graph/test_undirected_graph.py @@ -155,3 +155,21 @@ def test_undirected_graph_freeze_is_callable_after_inserts_and_after_reads(): np.testing.assert_array_equal( np.sort(graph.node_adjacency(0)[:, 0]), np.array([1, 2], 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) + + +@pytest.mark.parametrize( + "uvs", + [ + np.array([[1, 0]], np.uint64), + np.array([[0, 1], [0, 1]], np.uint64), + np.array([[1, 2], [0, 1]], np.uint64), + ], +) +def test_from_unique_edges_validates_fast_path_preconditions(uvs): + with pytest.raises(ValueError): + bic.graph.UndirectedGraph.from_unique_edges(3, uvs) diff --git a/tests/segmentation/test_relabel_sequential.py b/tests/segmentation/test_relabel_sequential.py index 9311179..a9f887b 100644 --- a/tests/segmentation/test_relabel_sequential.py +++ b/tests/segmentation/test_relabel_sequential.py @@ -219,3 +219,40 @@ def test_rejects_non_int_offset(): def test_rejects_bool_offset(): with pytest.raises(ValueError, match="offset must be a positive integer"): bic.segmentation.relabel_sequential(np.array([1, 2], dtype=np.uint32), offset=True) + + +def test_maximal_label_raises_instead_of_segfaulting(): + # max_value + 1 for the dense forward map would wrap size_t to zero; the + # kernel must reject this rather than write out of bounds. + labels = np.array([0, np.iinfo(np.uint64).max], dtype=np.uint64) + with pytest.raises(OverflowError): + bic.segmentation.relabel_sequential(labels) + + +def test_offset_plus_label_count_overflowing_dtype_raises(): + # offset == UINT32_MAX leaves room for exactly one new label; two distinct + # labels overflow the label dtype and must raise before allocating the + # multi-gigabyte inverse map. + labels = np.array([1, 2], dtype=np.uint32) + with pytest.raises(OverflowError): + bic.segmentation.relabel_sequential(labels, offset=int(np.iinfo(np.uint32).max)) + + +def test_offset_plus_label_count_overflowing_uint64_raises(): + labels = np.array([1, 2], dtype=np.uint64) + with pytest.raises(OverflowError): + bic.segmentation.relabel_sequential(labels, offset=int(np.iinfo(np.uint64).max)) + + +@pytest.mark.parametrize("dtype", SUPPORTED_DTYPES) +def test_maps_remain_plain_ndarrays(dtype): + # The minimal safe fix keeps the dense 3-tuple contract: forward_map and + # inverse_map are plain ndarrays of the input dtype. + labels = np.array([0, 5, 10, 5, 0, 200], dtype=dtype) + relabeled, forward_map, inverse_map = bic.segmentation.relabel_sequential(labels) + assert isinstance(forward_map, np.ndarray) + assert isinstance(inverse_map, np.ndarray) + assert forward_map.dtype == np.dtype(dtype) + assert inverse_map.dtype == np.dtype(dtype) + assert relabeled.dtype == np.dtype(dtype) + np.testing.assert_array_equal(inverse_map[relabeled], labels) diff --git a/tests/test_filters.py b/tests/test_filters.py index 13d6dc8..bce5a6d 100644 --- a/tests/test_filters.py +++ b/tests/test_filters.py @@ -13,6 +13,19 @@ import bioimage_cpp.filters as bf +@pytest.mark.parametrize("value", [np.nan, np.inf]) +def test_filters_reject_non_finite_inputs_and_parameters(value): + image = np.ones((5, 5), dtype=np.float32) + bad_image = image.copy() + bad_image[0, 0] = value + with pytest.raises(ValueError, match="finite"): + bf.gaussian_smoothing(bad_image, 1.0) + with pytest.raises(ValueError, match="finite"): + bf.gaussian_smoothing(image, value) + with pytest.raises(ValueError, match="finite"): + bf.gaussian_smoothing(image, 1.0, window_size=value) + + SHAPES_2D = [(7, 11), (32, 32), (48, 64)] SHAPES_3D = [(5, 7, 11), (8, 12, 16)] SIGMAS = [0.7, 1.5, 3.0] diff --git a/tests/test_flow.py b/tests/test_flow.py index adfb7a4..7b73cc1 100644 --- a/tests/test_flow.py +++ b/tests/test_flow.py @@ -153,6 +153,28 @@ def test_restrict_to_mask_keeps_density_inside_mask(): assert (density[~mask] == 0).all() +@pytest.mark.parametrize("method", ["euler", "rk2"]) +def test_restrict_to_mask_freezes_at_last_valid_position(method): + flow = np.zeros((2, 3, 3), dtype=np.float32) + flow[1] = 1.0 + mask = np.zeros((3, 3), dtype=bool) + mask[1, 1] = True + density = bic.flow.compute_flow_density( + flow, mask, n_iter=3, dt=1.0, tol=0.0, + method=method, restrict_to_mask=True, + ) + expected = np.zeros((3, 3), dtype=np.float32) + expected[1, 1] = 1.0 + np.testing.assert_array_equal(density, expected) + + +def test_flow_rejects_non_finite_values(): + flow = np.zeros((2, 3, 3), dtype=np.float32) + flow[0, 1, 1] = np.nan + with pytest.raises(ValueError, match="finite"): + bic.flow.compute_flow_density(flow, np.ones((3, 3), bool)) + + def test_rk2_runs(): rng = np.random.default_rng(0) flow = rng.normal(scale=0.3, size=(2, 12, 12)).astype(np.float32) diff --git a/tests/test_transformation.py b/tests/test_transformation.py index 5300220..d1f5c84 100644 --- a/tests/test_transformation.py +++ b/tests/test_transformation.py @@ -529,3 +529,38 @@ def test_map_coordinates_invalid_inputs_raise(): bic.transformation.map_coordinates(data, good, out=np.empty((3, 3), dtype=np.float32)) with pytest.raises(TypeError, match="dtype"): bic.transformation.map_coordinates(data, good, out=np.empty((4, 5), dtype=np.float64)) + + +@pytest.mark.parametrize("function", ["affine", "coordinates"]) +def test_resampling_supports_output_aliasing_input(function): + data = np.arange(9, dtype=np.int32).reshape(3, 3) + matrix = _matrix(2, translation=[-1, 0]) + if function == "affine": + expected = bic.transformation.affine_transform( + data, matrix, order=0, fill_value=-1 + ) + returned = bic.transformation.affine_transform( + data, matrix, order=0, fill_value=-1, out=data + ) + else: + coords = _affine_coords(matrix, data.shape) + expected = bic.transformation.map_coordinates( + data, coords, order=0, fill_value=-1 + ) + returned = bic.transformation.map_coordinates( + data, coords, order=0, fill_value=-1, out=data + ) + assert returned is data + np.testing.assert_array_equal(data, expected) + + +def test_non_finite_coordinates_use_fill_and_matrix_is_rejected(): + data = np.arange(9, dtype=np.float32).reshape(3, 3) + coords = np.indices(data.shape, dtype=np.float64) + coords[0, 1, 1] = np.nan + result = bic.transformation.map_coordinates(data, coords, order=0, fill_value=-9) + assert result[1, 1] == -9 + matrix = _matrix(2) + matrix[0, 2] = np.nan + with pytest.raises(ValueError, match="finite"): + bic.transformation.affine_transform(data, matrix)