diff --git a/.gitignore b/.gitignore index 1ce03dc..d495756 100644 --- a/.gitignore +++ b/.gitignore @@ -27,3 +27,5 @@ venv/ # Data *.h5 + +CLAUDE.md diff --git a/AGENTS.md b/AGENTS.md index 5279916..ed1d2c7 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -435,6 +435,11 @@ Public functions should have concise docstrings documenting: - background-label behavior, if relevant - axis and coordinate conventions +When adding or changing public functionality, update `MIGRATION_GUIDE.md` so +users migrating from `nifty`, `affogato`, or related packages know the +corresponding `bioimage-cpp` API, important behavioral differences, and any +intentional improvements. + ## What an agent should do when modifying this repository When adding or changing code: diff --git a/CMakeLists.txt b/CMakeLists.txt index dbc3141..e12d235 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,6 +11,7 @@ find_package(nanobind CONFIG REQUIRED) nanobind_add_module(_core NB_STATIC + src/bindings/blocking.cxx src/bindings/module.cxx src/bindings/graph.cxx src/bindings/segmentation.cxx diff --git a/MIGRATION_GUIDE.md b/MIGRATION_GUIDE.md new file mode 100644 index 0000000..22a4591 --- /dev/null +++ b/MIGRATION_GUIDE.md @@ -0,0 +1,293 @@ +# Migration Guide + +This guide explains how to migrate code that used selected `nifty` or +`affogato` functionality to `bioimage-cpp`. + +`bioimage-cpp` is not a drop-in compatibility layer. The package keeps a +smaller, NumPy-first API with Python-style method names, explicit dtype support, +and no I/O dependencies in the C++ core. + +## Imports + +Use: + +```python +import bioimage_cpp as bic +``` + +Graph functionality is under `bic.graph`, segmentation functionality is under +`bic.segmentation`, and small utility functions are under `bic.utils`. + +## Blocking + +Nifty: + +```python +import nifty.tools as nt + +blocking = nt.blocking([0, 0], [100, 80], [32, 32]) +block = blocking.getBlock(0) +block_with_halo = blocking.getBlockWithHalo(0, [8, 8]) +``` + +bioimage-cpp: + +```python +import bioimage_cpp as bic + +blocking = bic.Blocking([0, 0], [100, 80], [32, 32]) +block = blocking.get_block(0) +block_with_halo = blocking.get_block_with_halo(0, [8, 8]) +``` + +Name changes: + +| nifty-style name | bioimage-cpp name | +| --- | --- | +| `roiBegin` | `roi_begin` | +| `roiEnd` | `roi_end` | +| `blockShape` | `block_shape` | +| `blockShift` | `block_shift` | +| `blocksPerAxis` | `blocks_per_axis` | +| `numberOfBlocks` | `number_of_blocks` | +| `blockGridPosition` | `block_grid_position` | +| `getNeighborId` | `get_neighbor_id` | +| `getBlock` | `get_block` | +| `getBlockWithHalo` | `get_block_with_halo` | +| `addHalo` | `add_halo` | +| `coordinatesToBlockId` | `coordinates_to_block_id` | +| `getBlockIdsInBoundingBox` | `get_block_ids_in_bounding_box` | +| `getBlockIdsOverlappingBoundingBox` | `get_block_ids_overlapping_bounding_box` | +| `getLocalOverlaps` | `get_local_overlaps` | +| `getBlockIdsInSlice` | `get_block_ids_in_slice` | + +Intentional improvements over nifty: + +- `coordinates_to_block_id` accounts for both `roi_begin` and `block_shift`. +- `get_block_ids_overlapping_bounding_box` works for any dimensionality, not + only 3D. +- Bounding boxes use NumPy-style half-open intervals: `[begin, end)`. +- `get_local_overlaps` returns `None` if blocks do not overlap; otherwise it + returns `(begin_a, end_a, begin_b, end_b)` in local coordinates. + +## Undirected Graphs + +Nifty: + +```python +import nifty.graph as ng + +graph = ng.undirectedGraph(4) +edge_id = graph.insertEdge(0, 1) +uvs = graph.uvIds() +``` + +bioimage-cpp: + +```python +import bioimage_cpp as bic + +graph = bic.graph.UndirectedGraph(4) +edge_id = graph.insert_edge(0, 1) +uvs = graph.uv_ids() +``` + +The convenience constructor is: + +```python +graph = bic.graph.undirected_graph(4) +graph = bic.graph.UndirectedGraph.from_edges(4, [[0, 1], [1, 2]]) +``` + +Important differences: + +- Nodes are fixed at construction and have ids `0 .. number_of_nodes - 1`. +- Re-inserting an existing undirected edge returns the existing edge id. +- Bulk methods accept array-like inputs and return NumPy arrays. +- Python-style names are preferred. A few nifty-style aliases are still present + on `UndirectedGraph` for convenience, but new code should use snake_case. + +Common method/property mapping: + +| nifty-style name | bioimage-cpp name | +| --- | --- | +| `numberOfNodes` | `number_of_nodes` | +| `numberOfEdges` | `number_of_edges` | +| `nodeIdUpperBound` | `node_id_upper_bound` | +| `edgeIdUpperBound` | `edge_id_upper_bound` | +| `insertEdge` | `insert_edge` | +| `insertEdges` | `insert_edges` | +| `findEdge` | `find_edge` | +| `findEdges` | `find_edges` | +| `uvIds` | `uv_ids` | +| `nodeAdjacency` | `node_adjacency` | +| `serializationSize` | `serialization_size` | +| `extractSubgraphFromNodes` | `extract_subgraph_from_nodes` | +| `edgesFromNodeList` | `edges_from_node_list` | + +## Region Adjacency Graphs + +Nifty: + +```python +import nifty.graph.rag as nrag + +rag = nrag.gridRag(labels) +uvs = rag.uvIds() +``` + +bioimage-cpp: + +```python +import bioimage_cpp as bic + +rag = bic.graph.region_adjacency_graph(labels) +uvs = rag.uv_ids() +``` + +Notes: + +- Supported label dtypes are `uint32`, `uint64`, `int32`, and `int64`. +- Labels must be 2D or 3D. +- Negative signed labels are rejected. +- Nodes correspond to label ids from `0` to `labels.max()`. +- Edge ids are deterministic; RAG edges are sorted lexicographically by + endpoint ids. +- Non-contiguous labels are copied to contiguous memory before entering C++. + +## RAG Boundary and Affinity Features + +Nifty has RAG feature helpers such as `accumulateEdgeMeanAndLength`, +`accumulateEdgeStandartFeatures`, and affinity feature accumulation helpers. +In `bioimage-cpp`, these are exposed as explicit NumPy-returning functions. + +Simple edge-map features: + +```python +rag = bic.graph.region_adjacency_graph(labels) +features = bic.graph.edge_map_features(rag, labels, edge_map) +``` + +The columns are: + +```python +bic.graph.SIMPLE_EDGE_FEATURE_NAMES +# ("mean", "size") +``` + +Complex edge-map features: + +```python +features = bic.graph.edge_map_features_complex(rag, labels, edge_map) +``` + +The columns are: + +```python +bic.graph.COMPLEX_EDGE_FEATURE_NAMES +# ("mean", "median", "std", "min", "max", "p5", "p10", +# "p25", "p75", "p90", "p95", "size") +``` + +Affinity features: + +```python +features = bic.graph.affinity_features( + rag, + labels, + affinities, + offsets=[[0, 1], [1, 0]], +) +``` + +Complex affinity features: + +```python +features = bic.graph.affinity_features_complex( + rag, + labels, + affinities, + offsets=[[0, 1], [1, 0]], +) +``` + +Notes: + +- `edge_map` must have the same shape as `labels`. +- `affinities` must have shape `(channels, *labels.shape)`. +- `offsets` must have one offset per channel in NumPy axis order. +- Feature arrays use `float64` output. +- `number_of_threads=0` uses the library default; pass a positive integer for a + fixed thread count. + +## Mutex Watershed + +Affogato: + +```python +from affogato.segmentation import compute_mws_segmentation + +seg = compute_mws_segmentation( + weights, + offsets, + number_of_attractive_channels=3, + strides=[1, 1, 1], +) +``` + +bioimage-cpp: + +```python +import bioimage_cpp as bic + +seg = bic.segmentation.mutex_watershed( + weights, + offsets, + number_of_attractive_channels=3, + strides=[1, 1, 1], +) +``` + +Important migration notes: + +- `bioimage-cpp` expects the first `number_of_attractive_channels` channels to + be attractive merge-edge weights and the remaining channels to be mutex-edge + weights. +- Supported affinity dtypes are `float32` and `float64`. +- Inputs must represent 2D or 3D grids with shapes `(channels, y, x)` or + `(channels, z, y, x)`. +- Non-contiguous affinity arrays are copied to contiguous memory. +- `strides` sub-sample mutex edges only; attractive edges are always kept. +- `randomized_strides=True` uses NumPy's global random state, so existing + `np.random.seed(...)` workflows remain deterministic. +- A boolean `mask` may be passed. Edges touching `False` pixels are ignored and + masked pixels are set to label `0`. +- Output labels are `uint64`, consecutive, and 1-based for foreground pixels. + +## Dictionary-Based Relabeling + +If you used a small helper to apply a dictionary to an integer label array, use +`take_dict`: + +```python +labels = np.array([1, 3, 2, 1], dtype=np.uint64) +relabeling = {1: 10, 2: 20, 3: 30} + +out = bic.utils.take_dict(relabeling, labels) +``` + +Notes: + +- Supported input dtypes are `uint32`, `uint64`, `int32`, and `int64`. +- Output has the same shape and dtype as the input array. +- Every value in the input must be present in the mapping. +- Non-contiguous inputs are copied before entering C++. + +## I/O and Build Dependencies + +`bioimage-cpp` intentionally does not replace nifty or affogato I/O helpers. +Load TIFF, HDF5, zarr, N5, OME-NGFF, and related formats with existing Python +libraries, then pass NumPy arrays to `bioimage-cpp`. + +The package is designed for small PyPI wheels and does not depend on nifty, +vigra, HDF5, z5, xtensor, pybind11, or other large C++ libraries. diff --git a/include/bioimage_cpp/blocking.hxx b/include/bioimage_cpp/blocking.hxx new file mode 100644 index 0000000..a2ac2e3 --- /dev/null +++ b/include/bioimage_cpp/blocking.hxx @@ -0,0 +1,487 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace bioimage_cpp { + +using Coordinate = std::int64_t; +using CoordinateVector = std::vector; + +class Block { +public: + Block() = default; + + Block(CoordinateVector begin, CoordinateVector end) + : begin_(std::move(begin)), end_(std::move(end)) { + if (begin_.size() != end_.size()) { + throw std::invalid_argument("block begin and end must have the same length"); + } + for (std::size_t axis = 0; axis < begin_.size(); ++axis) { + if (end_[axis] < begin_[axis]) { + throw std::invalid_argument("block end must be >= begin for every axis"); + } + } + } + + const CoordinateVector &begin() const { + return begin_; + } + + const CoordinateVector &end() const { + return end_; + } + + CoordinateVector shape() const { + CoordinateVector result(begin_.size()); + for (std::size_t axis = 0; axis < begin_.size(); ++axis) { + result[axis] = end_[axis] - begin_[axis]; + } + return result; + } + + std::size_t ndim() const { + return begin_.size(); + } + +private: + CoordinateVector begin_; + CoordinateVector end_; +}; + +class BlockWithHalo { +public: + BlockWithHalo() = default; + + BlockWithHalo(Block outer_block, Block inner_block) + : outer_block_(std::move(outer_block)), + inner_block_(std::move(inner_block)), + inner_block_local_() { + if (outer_block_.ndim() != inner_block_.ndim()) { + throw std::invalid_argument("outer_block and inner_block must have the same ndim"); + } + + CoordinateVector local_begin(inner_block_.ndim()); + CoordinateVector local_end(inner_block_.ndim()); + const auto inner_shape = inner_block_.shape(); + for (std::size_t axis = 0; axis < inner_block_.ndim(); ++axis) { + local_begin[axis] = inner_block_.begin()[axis] - outer_block_.begin()[axis]; + local_end[axis] = local_begin[axis] + inner_shape[axis]; + } + inner_block_local_ = Block(std::move(local_begin), std::move(local_end)); + } + + const Block &outer_block() const { + return outer_block_; + } + + const Block &inner_block() const { + return inner_block_; + } + + const Block &inner_block_local() const { + return inner_block_local_; + } + +private: + Block outer_block_; + Block inner_block_; + Block inner_block_local_; +}; + +struct LocalOverlaps { + CoordinateVector overlap_begin_a; + CoordinateVector overlap_end_a; + CoordinateVector overlap_begin_b; + CoordinateVector overlap_end_b; +}; + +class Blocking { +public: + Blocking( + CoordinateVector roi_begin, + CoordinateVector roi_end, + CoordinateVector block_shape, + CoordinateVector block_shift = {} + ) + : roi_begin_(std::move(roi_begin)), + roi_end_(std::move(roi_end)), + block_shape_(std::move(block_shape)), + block_shift_(std::move(block_shift)), + blocks_per_axis_(), + blocks_per_axis_strides_(), + number_of_blocks_(1) { + validate_constructor_arguments(); + + const auto ndim = roi_begin_.size(); + blocks_per_axis_.assign(ndim, 0); + blocks_per_axis_strides_.assign(ndim, 1); + + for (std::size_t axis = 0; axis < ndim; ++axis) { + const auto dim_size = roi_end_[axis] - (roi_begin_[axis] - block_shift_[axis]); + blocks_per_axis_[axis] = ceil_div(dim_size, block_shape_[axis]); + number_of_blocks_ *= static_cast(blocks_per_axis_[axis]); + } + + for (std::ptrdiff_t axis = static_cast(ndim) - 2; axis >= 0; --axis) { + const auto index = static_cast(axis); + blocks_per_axis_strides_[index] = + blocks_per_axis_strides_[index + 1] * blocks_per_axis_[index + 1]; + } + } + + const CoordinateVector &roi_begin() const { + return roi_begin_; + } + + const CoordinateVector &roi_end() const { + return roi_end_; + } + + const CoordinateVector &block_shape() const { + return block_shape_; + } + + const CoordinateVector &block_shift() const { + return block_shift_; + } + + const CoordinateVector &blocks_per_axis() const { + return blocks_per_axis_; + } + + std::uint64_t number_of_blocks() const { + return number_of_blocks_; + } + + CoordinateVector block_grid_position(const std::uint64_t block_id) const { + require_valid_block_id(block_id); + CoordinateVector result(ndim()); + for (std::size_t axis = 0; axis < ndim(); ++axis) { + result[axis] = block_axis_position(block_id, axis); + } + return result; + } + + std::int64_t get_neighbor_id( + const std::uint64_t block_id, + const std::size_t axis, + const bool lower + ) const { + require_valid_block_id(block_id); + require_valid_axis(axis); + + const auto position = block_axis_position(block_id, axis); + if ((lower && position == 0) || (!lower && position == blocks_per_axis_[axis] - 1)) { + return -1; + } + + const auto stride = blocks_per_axis_strides_[axis]; + return static_cast(block_id) + (lower ? -stride : stride); + } + + Block get_block(const std::uint64_t block_id) const { + require_valid_block_id(block_id); + + CoordinateVector begin_coord(ndim()); + CoordinateVector end_coord(ndim()); + auto index = block_id; + for (std::size_t axis = 0; axis < ndim(); ++axis) { + const auto block_coord = static_cast(index / blocks_per_axis_strides_[axis]); + index -= static_cast(block_coord) * blocks_per_axis_strides_[axis]; + + const auto begin_at_axis = + (roi_begin_[axis] - block_shift_[axis]) + block_coord * block_shape_[axis]; + begin_coord[axis] = std::max(begin_at_axis, roi_begin_[axis]); + end_coord[axis] = std::min(begin_at_axis + block_shape_[axis], roi_end_[axis]); + } + return Block(std::move(begin_coord), std::move(end_coord)); + } + + BlockWithHalo get_block_with_halo( + const std::uint64_t block_id, + const CoordinateVector &halo_begin, + const CoordinateVector &halo_end + ) const { + const auto inner_block = get_block(block_id); + return add_halo(inner_block, halo_begin, halo_end); + } + + BlockWithHalo get_block_with_halo( + const std::uint64_t block_id, + const CoordinateVector &halo + ) const { + return get_block_with_halo(block_id, halo, halo); + } + + BlockWithHalo add_halo( + const Block &inner_block, + const CoordinateVector &halo_begin, + const CoordinateVector &halo_end + ) const { + require_vector_length(halo_begin, "halo_begin"); + require_vector_length(halo_end, "halo_end"); + if (inner_block.ndim() != ndim()) { + throw std::invalid_argument("inner_block ndim must match blocking ndim"); + } + + CoordinateVector outer_begin(ndim()); + CoordinateVector outer_end(ndim()); + for (std::size_t axis = 0; axis < ndim(); ++axis) { + if (halo_begin[axis] < 0 || halo_end[axis] < 0) { + throw std::invalid_argument("halo values must be non-negative"); + } + outer_begin[axis] = std::max(inner_block.begin()[axis] - halo_begin[axis], roi_begin_[axis]); + outer_end[axis] = std::min(inner_block.end()[axis] + halo_end[axis], roi_end_[axis]); + } + return BlockWithHalo(Block(std::move(outer_begin), std::move(outer_end)), inner_block); + } + + BlockWithHalo add_halo(const Block &inner_block, const CoordinateVector &halo) const { + return add_halo(inner_block, halo, halo); + } + + std::uint64_t coordinates_to_block_id(const CoordinateVector &coordinates) const { + require_vector_length(coordinates, "coordinates"); + + std::uint64_t block_id = 0; + for (std::size_t axis = 0; axis < ndim(); ++axis) { + if (coordinates[axis] < roi_begin_[axis] || coordinates[axis] >= roi_end_[axis]) { + throw std::out_of_range("coordinates must lie inside the half-open ROI"); + } + const auto shifted = coordinates[axis] - (roi_begin_[axis] - block_shift_[axis]); + const auto block_coord = shifted / block_shape_[axis]; + block_id += static_cast(block_coord) * blocks_per_axis_strides_[axis]; + } + return block_id; + } + + std::vector get_block_ids_in_bounding_box( + const CoordinateVector &box_begin, + const CoordinateVector &box_end + ) const { + require_box(box_begin, box_end); + + std::vector result; + for (std::uint64_t block_id = 0; block_id < number_of_blocks_; ++block_id) { + const auto block = get_block(block_id); + bool enclosed = true; + for (std::size_t axis = 0; axis < ndim(); ++axis) { + if (block.begin()[axis] < box_begin[axis] || block.end()[axis] > box_end[axis]) { + enclosed = false; + break; + } + } + if (enclosed) { + result.push_back(block_id); + } + } + return result; + } + + std::vector get_block_ids_overlapping_bounding_box( + const CoordinateVector &box_begin, + const CoordinateVector &box_end + ) const { + require_box(box_begin, box_end); + + CoordinateVector first(ndim()); + CoordinateVector last(ndim()); + for (std::size_t axis = 0; axis < ndim(); ++axis) { + const auto overlap_begin = std::max(box_begin[axis], roi_begin_[axis]); + const auto overlap_end = std::min(box_end[axis], roi_end_[axis]); + if (overlap_begin >= overlap_end) { + return {}; + } + first[axis] = block_axis_position_for_coordinate(overlap_begin, axis); + last[axis] = block_axis_position_for_coordinate(overlap_end - 1, axis); + } + + std::vector result; + CoordinateVector position = first; + while (true) { + std::uint64_t block_id = 0; + for (std::size_t axis = 0; axis < ndim(); ++axis) { + block_id += static_cast(position[axis]) * blocks_per_axis_strides_[axis]; + } + result.push_back(block_id); + + std::ptrdiff_t axis = static_cast(ndim()) - 1; + for (; axis >= 0; --axis) { + const auto index = static_cast(axis); + if (position[index] < last[index]) { + ++position[index]; + for (std::size_t reset_axis = index + 1; reset_axis < ndim(); ++reset_axis) { + position[reset_axis] = first[reset_axis]; + } + break; + } + } + if (axis < 0) { + break; + } + } + return result; + } + + std::optional get_local_overlaps( + const std::uint64_t block_a_id, + const std::uint64_t block_b_id, + const CoordinateVector &block_halo + ) const { + require_vector_length(block_halo, "block_halo"); + + const auto block_a = get_block_with_halo(block_a_id, block_halo).outer_block(); + const auto block_b = get_block_with_halo(block_b_id, block_halo).outer_block(); + + CoordinateVector global_begin(ndim()); + CoordinateVector global_end(ndim()); + for (std::size_t axis = 0; axis < ndim(); ++axis) { + global_begin[axis] = std::max(block_a.begin()[axis], block_b.begin()[axis]); + global_end[axis] = std::min(block_a.end()[axis], block_b.end()[axis]); + if (global_begin[axis] >= global_end[axis]) { + return std::nullopt; + } + } + + LocalOverlaps overlaps; + overlaps.overlap_begin_a.resize(ndim()); + overlaps.overlap_end_a.resize(ndim()); + overlaps.overlap_begin_b.resize(ndim()); + overlaps.overlap_end_b.resize(ndim()); + for (std::size_t axis = 0; axis < ndim(); ++axis) { + overlaps.overlap_begin_a[axis] = global_begin[axis] - block_a.begin()[axis]; + overlaps.overlap_end_a[axis] = global_end[axis] - block_a.begin()[axis]; + overlaps.overlap_begin_b[axis] = global_begin[axis] - block_b.begin()[axis]; + overlaps.overlap_end_b[axis] = global_end[axis] - block_b.begin()[axis]; + } + return overlaps; + } + + std::vector get_block_ids_in_slice( + const Coordinate z, + const CoordinateVector &block_halo + ) const { + require_vector_length(block_halo, "block_halo"); + if (ndim() == 0) { + return {}; + } + + std::vector result; + for (std::uint64_t block_id = 0; block_id < number_of_blocks_; ++block_id) { + const auto block = get_block_with_halo(block_id, block_halo).outer_block(); + if (z >= block.begin()[0] && z < block.end()[0]) { + result.push_back(block_id); + } + } + return result; + } + + std::size_t ndim() const { + return roi_begin_.size(); + } + +private: + static Coordinate ceil_div(const Coordinate numerator, const Coordinate denominator) { + return numerator == 0 ? 0 : (numerator - 1) / denominator + 1; + } + + void validate_constructor_arguments() { + if (roi_begin_.empty()) { + throw std::invalid_argument("roi_begin must not be empty"); + } + if (roi_end_.size() != roi_begin_.size()) { + throw std::invalid_argument("roi_end length must match roi_begin length"); + } + if (block_shape_.size() != roi_begin_.size()) { + throw std::invalid_argument("block_shape length must match roi_begin length"); + } + if (block_shift_.empty()) { + block_shift_.assign(roi_begin_.size(), 0); + } + if (block_shift_.size() != roi_begin_.size()) { + throw std::invalid_argument("block_shift length must match roi_begin length"); + } + + for (std::size_t axis = 0; axis < roi_begin_.size(); ++axis) { + if (roi_end_[axis] < roi_begin_[axis]) { + throw std::invalid_argument("roi_end must be >= roi_begin for every axis"); + } + if (block_shape_[axis] <= 0) { + throw std::invalid_argument("block_shape values must be positive"); + } + if (block_shift_[axis] < 0) { + throw std::invalid_argument("block_shift values must be non-negative"); + } + } + } + + void require_vector_length(const CoordinateVector &coordinates, const char *name) const { + if (coordinates.size() != ndim()) { + throw std::invalid_argument( + std::string(name) + " length must match blocking ndim" + ); + } + } + + void require_box(const CoordinateVector &box_begin, const CoordinateVector &box_end) const { + require_vector_length(box_begin, "box_begin"); + require_vector_length(box_end, "box_end"); + for (std::size_t axis = 0; axis < ndim(); ++axis) { + if (box_end[axis] < box_begin[axis]) { + throw std::invalid_argument("box_end must be >= box_begin for every axis"); + } + } + } + + void require_valid_axis(const std::size_t axis) const { + if (axis >= ndim()) { + throw std::out_of_range("axis is out of range"); + } + } + + void require_valid_block_id(const std::uint64_t block_id) const { + if (block_id >= number_of_blocks_) { + throw std::out_of_range("block_id is out of range"); + } + } + + Coordinate block_axis_position(const std::uint64_t block_id, const std::size_t axis) const { + require_valid_axis(axis); + auto index = block_id; + Coordinate position = 0; + for (std::size_t current_axis = 0; current_axis <= axis; ++current_axis) { + position = static_cast(index / blocks_per_axis_strides_[current_axis]); + index -= static_cast(position) * blocks_per_axis_strides_[current_axis]; + } + return position; + } + + Coordinate block_axis_position_for_coordinate( + const Coordinate coordinate, + const std::size_t axis + ) const { + const auto shifted = coordinate - (roi_begin_[axis] - block_shift_[axis]); + auto position = shifted / block_shape_[axis]; + position = std::max(0, position); + position = std::min(position, blocks_per_axis_[axis] - 1); + return position; + } + + CoordinateVector roi_begin_; + CoordinateVector roi_end_; + CoordinateVector block_shape_; + CoordinateVector block_shift_; + CoordinateVector blocks_per_axis_; + std::vector blocks_per_axis_strides_; + std::uint64_t number_of_blocks_; +}; + +} // namespace bioimage_cpp diff --git a/src/bindings/blocking.cxx b/src/bindings/blocking.cxx new file mode 100644 index 0000000..5090193 --- /dev/null +++ b/src/bindings/blocking.cxx @@ -0,0 +1,175 @@ +#include "blocking.hxx" + +#include "bioimage_cpp/blocking.hxx" + +#include +#include +#include + +#include +#include + +namespace nb = nanobind; + +namespace bioimage_cpp::bindings { +namespace { + +BlockWithHalo get_block_with_symmetric_halo( + const Blocking &blocking, + const std::uint64_t block_id, + const CoordinateVector &halo +) { + return blocking.get_block_with_halo(block_id, halo); +} + +BlockWithHalo get_block_with_asymmetric_halo( + const Blocking &blocking, + const std::uint64_t block_id, + const CoordinateVector &halo_begin, + const CoordinateVector &halo_end +) { + return blocking.get_block_with_halo(block_id, halo_begin, halo_end); +} + +BlockWithHalo add_symmetric_halo( + const Blocking &blocking, + const Block &inner_block, + const CoordinateVector &halo +) { + return blocking.add_halo(inner_block, halo); +} + +BlockWithHalo add_asymmetric_halo( + const Blocking &blocking, + const Block &inner_block, + const CoordinateVector &halo_begin, + const CoordinateVector &halo_end +) { + return blocking.add_halo(inner_block, halo_begin, halo_end); +} + +nb::object get_local_overlaps( + const Blocking &blocking, + const std::uint64_t block_a_id, + const std::uint64_t block_b_id, + const CoordinateVector &block_halo +) { + const auto overlaps = blocking.get_local_overlaps(block_a_id, block_b_id, block_halo); + if (!overlaps.has_value()) { + return nb::none(); + } + return nb::cast(std::make_tuple( + overlaps->overlap_begin_a, + overlaps->overlap_end_a, + overlaps->overlap_begin_b, + overlaps->overlap_end_b + )); +} + +} // namespace + +void bind_blocking(nb::module_ &m) { + nb::class_(m, "Block") + .def( + nb::init(), + nb::arg("begin"), + nb::arg("end") + ) + .def_prop_ro("begin", &Block::begin) + .def_prop_ro("end", &Block::end) + .def_prop_ro("shape", &Block::shape) + .def_prop_ro("ndim", &Block::ndim); + + nb::class_(m, "BlockWithHalo") + .def( + nb::init(), + nb::arg("outer_block"), + nb::arg("inner_block") + ) + .def_prop_ro("outer_block", &BlockWithHalo::outer_block) + .def_prop_ro("inner_block", &BlockWithHalo::inner_block) + .def_prop_ro("inner_block_local", &BlockWithHalo::inner_block_local); + + nb::class_(m, "Blocking") + .def( + nb::init(), + nb::arg("roi_begin"), + nb::arg("roi_end"), + nb::arg("block_shape"), + nb::arg("block_shift") = CoordinateVector{} + ) + .def_prop_ro("roi_begin", &Blocking::roi_begin) + .def_prop_ro("roi_end", &Blocking::roi_end) + .def_prop_ro("block_shape", &Blocking::block_shape) + .def_prop_ro("block_shift", &Blocking::block_shift) + .def_prop_ro("blocks_per_axis", &Blocking::blocks_per_axis) + .def_prop_ro("number_of_blocks", &Blocking::number_of_blocks) + .def_prop_ro("ndim", &Blocking::ndim) + .def("block_grid_position", &Blocking::block_grid_position, nb::arg("block_id")) + .def( + "get_neighbor_id", + &Blocking::get_neighbor_id, + nb::arg("block_id"), + nb::arg("axis"), + nb::arg("lower") + ) + .def("get_block", &Blocking::get_block, nb::arg("block_id")) + .def( + "get_block_with_halo", + &get_block_with_symmetric_halo, + nb::arg("block_id"), + nb::arg("halo") + ) + .def( + "get_block_with_halo", + &get_block_with_asymmetric_halo, + nb::arg("block_id"), + nb::arg("halo_begin"), + nb::arg("halo_end") + ) + .def( + "add_halo", + &add_symmetric_halo, + nb::arg("inner_block"), + nb::arg("halo") + ) + .def( + "add_halo", + &add_asymmetric_halo, + nb::arg("inner_block"), + nb::arg("halo_begin"), + nb::arg("halo_end") + ) + .def( + "coordinates_to_block_id", + &Blocking::coordinates_to_block_id, + nb::arg("coordinates") + ) + .def( + "get_block_ids_in_bounding_box", + &Blocking::get_block_ids_in_bounding_box, + nb::arg("box_begin"), + nb::arg("box_end") + ) + .def( + "get_block_ids_overlapping_bounding_box", + &Blocking::get_block_ids_overlapping_bounding_box, + nb::arg("box_begin"), + nb::arg("box_end") + ) + .def( + "get_local_overlaps", + &get_local_overlaps, + nb::arg("block_a_id"), + nb::arg("block_b_id"), + nb::arg("block_halo") + ) + .def( + "get_block_ids_in_slice", + &Blocking::get_block_ids_in_slice, + nb::arg("z"), + nb::arg("block_halo") + ); +} + +} // namespace bioimage_cpp::bindings diff --git a/src/bindings/blocking.hxx b/src/bindings/blocking.hxx new file mode 100644 index 0000000..1da3704 --- /dev/null +++ b/src/bindings/blocking.hxx @@ -0,0 +1,9 @@ +#pragma once + +#include + +namespace bioimage_cpp::bindings { + +void bind_blocking(nanobind::module_ &m); + +} // namespace bioimage_cpp::bindings diff --git a/src/bindings/module.cxx b/src/bindings/module.cxx index 0db2821..8b9206a 100644 --- a/src/bindings/module.cxx +++ b/src/bindings/module.cxx @@ -1,3 +1,4 @@ +#include "blocking.hxx" #include "graph.hxx" #include "segmentation.hxx" #include "utils.hxx" @@ -8,6 +9,7 @@ namespace nb = nanobind; NB_MODULE(_core, m) { m.doc() = "C++ extension module for bioimage_cpp."; + bioimage_cpp::bindings::bind_blocking(m); bioimage_cpp::bindings::bind_graph(m); bioimage_cpp::bindings::bind_segmentation(m); bioimage_cpp::bindings::bind_utils(m); diff --git a/src/bioimage_cpp/__init__.py b/src/bioimage_cpp/__init__.py index 2eba9e0..4879fb3 100644 --- a/src/bioimage_cpp/__init__.py +++ b/src/bioimage_cpp/__init__.py @@ -1,8 +1,17 @@ """Dependency-light bioimage analysis algorithms backed by C++.""" from ._version import __version__ +from ._core import Block, Blocking, BlockWithHalo from . import graph from . import segmentation from . import utils -__all__ = ["__version__", "graph", "segmentation", "utils"] +__all__ = [ + "__version__", + "Block", + "Blocking", + "BlockWithHalo", + "graph", + "segmentation", + "utils", +] diff --git a/tests/test_blocking.py b/tests/test_blocking.py new file mode 100644 index 0000000..c7e670f --- /dev/null +++ b/tests/test_blocking.py @@ -0,0 +1,137 @@ +import pytest + +import bioimage_cpp as bic + + +def test_block_shape_and_halo_local_coordinates(): + block = bic.Block([2, 3], [7, 11]) + + assert block.begin == [2, 3] + assert block.end == [7, 11] + assert block.shape == [5, 8] + assert block.ndim == 2 + + with_halo = bic.BlockWithHalo(bic.Block([0, 1], [9, 12]), block) + assert with_halo.outer_block.begin == [0, 1] + assert with_halo.outer_block.end == [9, 12] + assert with_halo.inner_block.begin == [2, 3] + assert with_halo.inner_block_local.begin == [2, 2] + assert with_halo.inner_block_local.end == [7, 10] + + +def test_blocking_grid_and_blocks_match_nifty_layout(): + blocking = bic.Blocking([0, 0], [10, 7], [4, 3]) + + assert blocking.roi_begin == [0, 0] + assert blocking.roi_end == [10, 7] + assert blocking.block_shape == [4, 3] + assert blocking.block_shift == [0, 0] + assert blocking.blocks_per_axis == [3, 3] + assert blocking.number_of_blocks == 9 + assert blocking.ndim == 2 + + assert blocking.block_grid_position(0) == [0, 0] + assert blocking.block_grid_position(5) == [1, 2] + assert blocking.get_block(0).begin == [0, 0] + assert blocking.get_block(0).end == [4, 3] + assert blocking.get_block(5).begin == [4, 6] + assert blocking.get_block(5).end == [8, 7] + assert blocking.get_block(8).begin == [8, 6] + assert blocking.get_block(8).end == [10, 7] + + +def test_neighbor_ids(): + blocking = bic.Blocking([0, 0], [10, 7], [4, 3]) + + assert blocking.get_neighbor_id(4, axis=0, lower=True) == 1 + assert blocking.get_neighbor_id(4, axis=0, lower=False) == 7 + assert blocking.get_neighbor_id(4, axis=1, lower=True) == 3 + assert blocking.get_neighbor_id(4, axis=1, lower=False) == 5 + assert blocking.get_neighbor_id(1, axis=0, lower=True) == -1 + assert blocking.get_neighbor_id(7, axis=0, lower=False) == -1 + assert blocking.get_neighbor_id(3, axis=1, lower=True) == -1 + assert blocking.get_neighbor_id(5, axis=1, lower=False) == -1 + + +def test_shifted_roi_coordinates_to_block_id_and_overlapping_box(): + blocking = bic.Blocking([10, 20], [20, 31], [4, 5], [2, 1]) + + assert blocking.blocks_per_axis == [3, 3] + assert blocking.get_block(0).begin == [10, 20] + assert blocking.get_block(0).end == [12, 24] + assert blocking.get_block(4).begin == [12, 24] + assert blocking.get_block(4).end == [16, 29] + + assert blocking.coordinates_to_block_id([10, 20]) == 0 + assert blocking.coordinates_to_block_id([12, 24]) == 4 + assert blocking.coordinates_to_block_id([19, 30]) == 8 + + assert blocking.get_block_ids_overlapping_bounding_box([11, 23], [17, 30]) == list( + range(9) + ) + + +def test_block_ids_in_bounding_box_requires_full_enclosure(): + blocking = bic.Blocking([0, 0], [10, 10], [4, 4]) + + assert blocking.get_block_ids_in_bounding_box([4, 4], [10, 10]) == [4, 5, 7, 8] + assert blocking.get_block_ids_in_bounding_box([5, 5], [10, 10]) == [8] + + +def test_overlapping_bounding_box_is_dimension_independent(): + blocking = bic.Blocking([0, 0, 0, 0], [6, 6, 6, 6], [3, 3, 3, 3]) + + assert blocking.blocks_per_axis == [2, 2, 2, 2] + assert blocking.get_block_ids_overlapping_bounding_box( + [2, 2, 2, 2], [4, 4, 4, 4] + ) == list(range(16)) + assert blocking.get_block_ids_overlapping_bounding_box( + [0, 0, 0, 0], [3, 3, 3, 3] + ) == [0] + + +def test_block_with_halo_and_add_halo_clip_to_roi(): + blocking = bic.Blocking([0, 0], [10, 10], [4, 4]) + + block = blocking.get_block_with_halo(0, [2, 1]) + assert block.inner_block.begin == [0, 0] + assert block.inner_block.end == [4, 4] + assert block.outer_block.begin == [0, 0] + assert block.outer_block.end == [6, 5] + assert block.inner_block_local.begin == [0, 0] + assert block.inner_block_local.end == [4, 4] + + asymmetric = blocking.add_halo(blocking.get_block(4), [1, 2], [3, 4]) + assert asymmetric.outer_block.begin == [3, 2] + assert asymmetric.outer_block.end == [10, 10] + assert asymmetric.inner_block_local.begin == [1, 2] + assert asymmetric.inner_block_local.end == [5, 6] + + +def test_local_overlaps_return_local_coordinates_or_none(): + blocking = bic.Blocking([0, 0], [10, 10], [5, 5]) + + overlaps = blocking.get_local_overlaps(0, 1, [1, 1]) + assert overlaps == ([0, 4], [6, 6], [0, 0], [6, 2]) + + assert blocking.get_local_overlaps(0, 3, [0, 0]) is None + + +def test_block_ids_in_slice_uses_first_axis_and_halo(): + blocking = bic.Blocking([0, 0, 0], [6, 4, 4], [3, 2, 2]) + + assert blocking.get_block_ids_in_slice(3, [0, 0, 0]) == list(range(4, 8)) + assert blocking.get_block_ids_in_slice(3, [1, 0, 0]) == list(range(8)) + + +def test_invalid_inputs_raise_clear_errors(): + with pytest.raises(ValueError, match="block_shape values must be positive"): + bic.Blocking([0, 0], [10, 10], [4, 0]) + + blocking = bic.Blocking([0, 0], [10, 10], [4, 4]) + with pytest.raises(IndexError, match="block_id is out of range"): + blocking.get_block(9) + with pytest.raises(IndexError, match="coordinates must lie inside"): + blocking.coordinates_to_block_id([10, 0]) + with pytest.raises(ValueError, match="box_end must be >= box_begin"): + blocking.get_block_ids_overlapping_bounding_box([5, 0], [4, 1])