From 22c75cca7dd7c99de3792650e6139a129f0b0992 Mon Sep 17 00:00:00 2001 From: Valentin Boussot Date: Thu, 6 Aug 2026 15:10:52 +0200 Subject: [PATCH 1/7] feat(py): add itk_transform_resample_bounding_box for out-of-core resampling Compute the region of a moving image needed to resample a fixed image grid, from geometry alone: the pixel buffers are never read and the Dask graphs are never computed. That is the point -- describe two images and a transform with a few numbers, learn exactly which block a resample will touch, and only then move pixels. The transform is an ITK one, including the CompositeTransform an Elastix registration returns, and it maps fixed points into moving space. The image geometry is built the way ngff_image_to_itk_image builds it, RFC-4 direction included, so the transform is applied in the space it was produced in. ResampleBoundingBox keys everything by dimension name -- the pipeline reports arrays fastest-axis-first, the reverse of the Zarr order -- clamps negative start indices rather than letting them wrap, and crop() returns a lazily sliced NgffImage with a corrected translation. --- py/ngff_zarr/__init__.py | 7 + .../itk_transform_resample_bounding_box.py | 448 +++++++++++++++ ...est_itk_transform_resample_bounding_box.py | 516 ++++++++++++++++++ 3 files changed, 971 insertions(+) create mode 100644 py/ngff_zarr/itk_transform_resample_bounding_box.py create mode 100644 py/test/test_itk_transform_resample_bounding_box.py diff --git a/py/ngff_zarr/__init__.py b/py/ngff_zarr/__init__.py index f759d124..72914843 100644 --- a/py/ngff_zarr/__init__.py +++ b/py/ngff_zarr/__init__.py @@ -24,6 +24,10 @@ write_hcs_well_image, ) from .itk_image_to_ngff_image import itk_image_to_ngff_image +from .itk_transform_resample_bounding_box import ( + ResampleBoundingBox, + itk_transform_resample_bounding_box, +) from .lif_to_ngff_image import ( has_mosaic_dimension, lif_file_to_ngff_images, @@ -118,6 +122,9 @@ "nibabel_image_to_ngff_image", "extract_omero_metadata_from_nibabel", "ngff_image_to_itk_image", + # Out-of-core resampling + "itk_transform_resample_bounding_box", + "ResampleBoundingBox", "memory_usage", "task_count", "to_multiscales", diff --git a/py/ngff_zarr/itk_transform_resample_bounding_box.py b/py/ngff_zarr/itk_transform_resample_bounding_box.py new file mode 100644 index 00000000..f7ac9066 --- /dev/null +++ b/py/ngff_zarr/itk_transform_resample_bounding_box.py @@ -0,0 +1,448 @@ +# SPDX-FileCopyrightText: Copyright (c) Fideus Labs LLC +# SPDX-License-Identifier: MIT +"""Find the region of a moving image needed to resample a fixed image grid.""" + +import math +from collections.abc import Sequence +from dataclasses import dataclass, field + +import numpy as np + +from .ngff_image import NgffImage +from .rfc4 import anatomical_orientation_to_itk_direction + +_SPATIAL_DIMS = ("x", "y", "z") + + +@dataclass +class ResampleBoundingBox: + """The region of a moving image needed to resample a fixed image grid. + + Every mapping is keyed by dimension name, so callers never have to track + whether a given array is in Zarr or ITK axis order. + """ + + #: Spatial dimension names, in the moving image's (Zarr) order. + dims: tuple[str, ...] + #: Start of the region in the moving image's index space. May be negative + #: when the transformed grid extends past the moving image origin. + start_index: dict[str, int] + #: Size of the region in moving-image pixels, before clamping. + size: dict[str, int] + #: Tight extent of the transformed fixed grid, in moving-image physical + #: coordinates. Independent of ``padding``. + corners_min: dict[str, float] + corners_max: dict[str, float] + #: Physical coordinates of the padded region's first and last pixel. + padded_corners_min: dict[str, float] + padded_corners_max: dict[str, float] + #: Shape of the moving image, used to clamp the region into bounds. + moving_shape: dict[str, int] = field(default_factory=dict) + + def clamped(self) -> dict[str, tuple[int, int]]: + """Return the region intersected with the moving image bounds. + + :return: ``{dim: (start, stop)}`` with ``0 <= start <= stop <= shape``. + :rtype: dict[str, tuple[int, int]] + """ + bounds = {} + for dim in self.dims: + extent = self.moving_shape[dim] + start = max(0, min(self.start_index[dim], extent)) + stop = max(start, min(self.start_index[dim] + self.size[dim], extent)) + bounds[dim] = (start, stop) + return bounds + + @property + def is_empty(self) -> bool: + """Whether the region does not overlap the moving image at all.""" + return any(stop <= start for start, stop in self.clamped().values()) + + def slices(self, dims: Sequence[str] | None = None) -> tuple[slice, ...]: + """Slices selecting the region, clamped to the moving image bounds. + + Negative start indices are clamped rather than passed through, since a + negative bound would silently wrap around instead of raising. + + :param dims: Dimension order of the array to be sliced. Defaults to the + spatial dims. Dimensions outside the region (``t``, ``c``) get a + full slice. + :type dims: Sequence[str] | None + + :return: One slice per entry in ``dims``. + :rtype: tuple[slice, ...] + """ + order = tuple(dims) if dims is not None else self.dims + bounds = self.clamped() + return tuple( + slice(*bounds[dim]) if dim in bounds else slice(None) for dim in order + ) + + def crop(self, moving: NgffImage) -> NgffImage | None: + """Lazily crop ``moving`` to this region, correcting its translation. + + The returned image indexes the same Dask graph -- nothing is computed -- + so only this region's chunks are read when it is finally materialized. + + :param moving: The moving image the region was computed against. + :type moving: NgffImage + + :return: The cropped image, or ``None`` when the region does not + overlap the moving image. + :rtype: NgffImage | None + """ + if self.is_empty: + return None + bounds = self.clamped() + translation = dict(moving.translation) + for dim, (start, _) in bounds.items(): + translation[dim] = moving.translation[dim] + start * moving.scale[dim] + return NgffImage( + data=moving.data[self.slices(moving.dims)], + dims=moving.dims, + scale=dict(moving.scale), + translation=translation, + name=moving.name, + axes_units=moving.axes_units, + axes_orientations=moving.axes_orientations, + axes_types=moving.axes_types, + channel_names=moving.channel_names, + channel_colors=moving.channel_colors, + ) + + +def _spatial_dims(ngff_image: NgffImage) -> list[str]: + return [dim for dim in ngff_image.dims if dim in _SPATIAL_DIMS] + + +def _check_geometry(label: str, image: NgffImage, spatial: Sequence[str]) -> None: + """Require a finite scale and translation for every spatial axis. + + A missing or non-finite entry would otherwise reach the pipeline and come + back as a plausible-looking region computed from the wrong geometry. + """ + for name, values in (("scale", image.scale), ("translation", image.translation)): + for dim in spatial: + if dim not in values: + msg = f"{label} image {name} has no entry for dimension '{dim}'" + raise ValueError(msg) + if not np.isfinite(values[dim]): + msg = ( + f"{label} image {name} for dimension '{dim}' is " + f"{values[dim]}; it must be finite" + ) + raise ValueError(msg) + if name == "scale": + for dim in spatial: + if values[dim] == 0: + msg = f"{label} image scale for dimension '{dim}' is zero" + raise ValueError(msg) + + +def _check_index_arrays(result: dict, itk_dims: Sequence[str]) -> None: + """Reject a pipeline result that cannot fill one bound per dimension. + + The returned mappings are built from these two arrays alone, so a short or + non-finite result would otherwise surface as an opaque failure further + down, or as an undefined bound a caller slices with. + """ + for name in ("paddedStartIndex", "paddedSize"): + values = result[name] + if len(values) != len(itk_dims): + msg = ( + f"the pipeline reported {len(values)} {name} values for " + f"{len(itk_dims)} dimensions" + ) + raise ValueError(msg) + for dim, value in zip(itk_dims, values): + if not math.isfinite(value): + msg = ( + f"the {name} for dimension '{dim}' is not finite; check " + "the transform and the image scale and translation" + ) + raise ValueError(msg) + + +def _check_region_contains_corners(result: dict, itk_dims: Sequence[str]) -> None: + """Reject a region that does not contain the points it was derived from. + + The pipeline computes the integer region in 32-bit index space while + reporting the corners as doubles, so a grid whose extent exceeds that range + comes back as a wrapped -- and typically empty -- region. Padding is + non-negative, so a correct region always spans its own tight corners; + checking that catches the wrap without assuming an axis direction. + """ + corners_min = result["corners"]["min"] + corners_max = result["corners"]["max"] + padded_min = result["paddedCorners"]["min"] + padded_max = result["paddedCorners"]["max"] + + for axis, dim in enumerate(itk_dims): + values = ( + corners_min[axis], + corners_max[axis], + padded_min[axis], + padded_max[axis], + ) + if any(value is None or not np.isfinite(value) for value in values): + msg = ( + f"the resample region for dimension '{dim}' is not finite; " + "check the transform and the image scale and translation" + ) + raise ValueError(msg) + # paddedCorners hold the start and end index positions, which invert + # when an axis direction is negative, so compare against the span. + low = min(padded_min[axis], padded_max[axis]) + high = max(padded_min[axis], padded_max[axis]) + tight_low = min(corners_min[axis], corners_max[axis]) + tight_high = max(corners_min[axis], corners_max[axis]) + if tight_low < low or tight_high > high: + msg = ( + f"the resample region for dimension '{dim}' does not contain " + f"the transformed grid ({tight_low} to {tight_high} lies " + f"outside {low} to {high}). The region spans more than the " + "index range the pipeline can represent; check for a scale " + "mismatch between the fixed and moving images." + ) + raise ValueError(msg) + + +def _itk_direction(ngff_image: NgffImage, itk_dims: Sequence[str]) -> np.ndarray: + """Direction matrix from RFC-4 orientation, matching ngff_image_to_itk_image. + + All-or-nothing: unless every spatial axis carries an orientation that maps + onto an LPS axis, the direction falls back to identity. + """ + direction = np.eye(len(itk_dims)) + orientations = ngff_image.axes_orientations + if not orientations: + return direction + + columns = [] + for dim in itk_dims: + orientation = orientations.get(dim) + if orientation is None: + return direction + column = anatomical_orientation_to_itk_direction(orientation.value) + if column is None: + return direction + columns.append(column) + + for col_index, column in enumerate(columns): + for row in range(len(itk_dims)): + direction[row, col_index] = column[row] + return direction + + +def _metadata_only_itk_image( + ngff_image: NgffImage, + itk_dims: Sequence[str], + direction: np.ndarray, +): + """Build an ITK-Wasm image carrying geometry only, with an empty buffer. + + ``ngff_image.data`` is never converted to an array -- only its ``shape`` and + ``dtype`` are read, which are Dask metadata. This is what lets the region be + computed against images whose pixels are remote or enormous. + """ + from itkwasm import Image, ImageType, IntTypes, PixelTypes + + dims = tuple(ngff_image.dims) + image_type = ImageType( + dimension=len(itk_dims), + componentType=IntTypes.UInt8, + pixelType=PixelTypes.Scalar, + components=1, + ) + return Image( + imageType=image_type, + name=ngff_image.name, + origin=[float(ngff_image.translation[dim]) for dim in itk_dims], + spacing=[float(ngff_image.scale[dim]) for dim in itk_dims], + direction=direction, + size=[int(ngff_image.data.shape[dims.index(dim)]) for dim in itk_dims], + metadata={}, + data=np.empty((0,), dtype=np.uint8), + ) + + +def _as_itk_transform_list(transform) -> list: + """Normalize a supported transform input into an ITK-Wasm transform list. + + ITK applies the *last* entry of a transform list first, and so does an + ``itk.CompositeTransform``, so the order carries over unchanged. + """ + from itkwasm import Transform as ItkTransform + + if isinstance(transform, ItkTransform): + return [transform] + + if isinstance(transform, dict): + return [_transform_from_dict(transform)] + + if isinstance(transform, (list, tuple)) and transform: + entries = list(transform) + if all(isinstance(entry, ItkTransform) for entry in entries): + return entries + if all(isinstance(entry, dict) for entry in entries): + return [_transform_from_dict(entry) for entry in entries] + + # An itk.Transform, including the CompositeTransform Elastix returns. + if hasattr(transform, "GetTransformTypeAsString"): + import itk + + as_dict = itk.dict_from_transform(transform) + if isinstance(as_dict, dict): + as_dict = [as_dict] + return [_transform_from_dict(entry) for entry in as_dict] + + msg = ( + f"unsupported transform type {type(transform).__name__}. Expected an " + "RFC-5 coordinate transformation, an itk.Transform, or an ITK-Wasm " + "Transform / TransformList." + ) + raise TypeError(msg) + + +def _transform_from_dict(entry: dict): + from itkwasm import Transform as ItkTransform + from itkwasm import TransformType + + entry = dict(entry) + transform_type = entry.get("transformType") + if isinstance(transform_type, dict): + entry["transformType"] = TransformType(**transform_type) + return ItkTransform(**entry) + + +def itk_transform_resample_bounding_box( + transform, + fixed: NgffImage, + moving: NgffImage, + padding: int = 1, +) -> ResampleBoundingBox: + """Compute the moving-image region needed to resample a fixed image grid. + + The region is derived from image geometry alone -- the pixel buffers of + ``fixed`` and ``moving`` are never read, and their Dask graphs are never + computed. That is the point: describe two images and a transform with a few + numbers, learn exactly which block of the moving image a resample will + touch, and only then move pixels. + + The transform acts on ITK physical space, so the image geometry is built + the way :func:`ngff_zarr.ngff_image_to_itk_image` builds it, including the + direction matrix derived from RFC-4 anatomical orientation. It maps *fixed* + points into *moving* space, matching the direction registration libraries + return. + + :param transform: An ``itk.Transform`` (including the ``CompositeTransform`` + an Elastix registration returns), or an ITK-Wasm ``Transform`` / + ``TransformList``. + + :param fixed: The image whose grid is resampled. Geometry only. + :type fixed: NgffImage + + :param moving: The image to be sampled. Geometry only. + :type moving: NgffImage + + :param padding: Pixels of padding added per side. The default of 1 covers + linear interpolation, which reads one neighbor beyond the continuous + index bound. Use 0 for the tight region, or more for wider kernels. + :type padding: int + + :return: The region, keyed by dimension name in Zarr order. + :rtype: ResampleBoundingBox + + :raises ValueError: If ``fixed`` and ``moving`` do not share the same + spatial dimensions, if there are not 2, 3 or 4 of them, if a scale or + translation entry is missing, zero or non-finite, if ``padding`` is + negative, if the transform couples spatial and non-spatial axes, or if + the region spans more than the index range the pipeline can represent. + """ + from itkwasm_downsample import resample_bounding_box + + if not isinstance(padding, int) or isinstance(padding, bool) or padding < 0: + msg = f"padding must be a non-negative integer, got {padding!r}" + raise ValueError(msg) + + fixed_spatial = _spatial_dims(fixed) + moving_spatial = _spatial_dims(moving) + if fixed_spatial != moving_spatial: + msg = ( + f"fixed image has spatial dims {tuple(fixed_spatial)} but moving " + f"image has {tuple(moving_spatial)}; they must match" + ) + raise ValueError(msg) + if len(fixed_spatial) < 2: + msg = ( + f"images have {len(fixed_spatial)} spatial dims " + f"{tuple(fixed_spatial)}; only 2 and 3 are supported" + ) + raise ValueError(msg) + + for label, image in (("fixed", fixed), ("moving", moving)): + _check_geometry(label, image, fixed_spatial) + + # ITK orders points fastest-axis-first, the reverse of the Zarr order. + itk_dims = list(reversed(fixed_spatial)) + + fixed_dims = tuple(fixed.dims) + fixed_extent = { + dim: int(fixed.data.shape[fixed_dims.index(dim)]) for dim in fixed_spatial + } + moving_dims = tuple(moving.dims) + moving_shape = { + dim: int(moving.data.shape[moving_dims.index(dim)]) for dim in moving_spatial + } + if any(extent == 0 for extent in fixed_extent.values()): + # A grid with a zero-length axis has no samples to resample, so there + # is nothing to fetch. The pipeline reports an all-zero region here; + # returning it directly keeps the corners meaningful rather than the + # sentinels a degenerate boundary walk produces. + zeros = dict.fromkeys(fixed_spatial, 0) + return ResampleBoundingBox( + dims=tuple(fixed_spatial), + start_index=dict(zeros), + size=dict(zeros), + corners_min=dict.fromkeys(fixed_spatial, 0.0), + corners_max=dict.fromkeys(fixed_spatial, 0.0), + padded_corners_min=dict.fromkeys(fixed_spatial, 0.0), + padded_corners_max=dict.fromkeys(fixed_spatial, 0.0), + moving_shape=moving_shape, + ) + + transform_list = _as_itk_transform_list(transform) + fixed_direction = _itk_direction(fixed, itk_dims) + moving_direction = _itk_direction(moving, itk_dims) + + result = resample_bounding_box( + transform_list, + _metadata_only_itk_image(fixed, itk_dims, fixed_direction), + _metadata_only_itk_image(moving, itk_dims, moving_direction), + padding=padding, + ) + + _check_index_arrays(result, itk_dims) + _check_region_contains_corners(result, itk_dims) + + # The pipeline reports arrays fastest-axis-first; key them by dimension + # name so no caller has to remember that, and order the mappings like + # ``dims`` so a printed result reads in Zarr order. + def by_dim(values) -> dict: + indexed = dict(zip(itk_dims, values)) + return {dim: indexed[dim] for dim in fixed_spatial} + + return ResampleBoundingBox( + dims=tuple(fixed_spatial), + start_index={k: int(v) for k, v in by_dim(result["paddedStartIndex"]).items()}, + size={k: int(v) for k, v in by_dim(result["paddedSize"]).items()}, + corners_min={k: float(v) for k, v in by_dim(result["corners"]["min"]).items()}, + corners_max={k: float(v) for k, v in by_dim(result["corners"]["max"]).items()}, + padded_corners_min={ + k: float(v) for k, v in by_dim(result["paddedCorners"]["min"]).items() + }, + padded_corners_max={ + k: float(v) for k, v in by_dim(result["paddedCorners"]["max"]).items() + }, + moving_shape=moving_shape, + ) diff --git a/py/test/test_itk_transform_resample_bounding_box.py b/py/test/test_itk_transform_resample_bounding_box.py new file mode 100644 index 00000000..754ca96b --- /dev/null +++ b/py/test/test_itk_transform_resample_bounding_box.py @@ -0,0 +1,516 @@ +# SPDX-FileCopyrightText: Copyright (c) Fideus Labs LLC +# SPDX-License-Identifier: MIT +"""Tests for the RFC-5 to ITK transform bridge and the resample bounding box. + +The expected regions here are not taken from the pipeline itself. They come +either from the worked examples in the ITK-Wasm ``resample-bounding-box`` +documentation, or from ``_oracle_region`` below, which recomputes the region +from first principles in NGFF axis order. That keeps the tests honest about +the two conventions that differ between RFC-5 and ITK: axis order and +sequence composition order. +""" + +import itertools + +import dask.array as da +import numpy as np +import pytest +from ngff_zarr import ( + RAS, + NgffImage, + itk_transform_resample_bounding_box, + ngff_image_to_itk_image, +) +from ngff_zarr.itk_transform_resample_bounding_box import ( + _itk_direction, + _metadata_only_itk_image, +) + + +def _translation(offset): + """An ITK-Wasm translation, in ITK (fastest-axis-first) order.""" + from itkwasm import ( + FloatTypes, + Transform, + TransformParameterizations, + TransformType, + ) + + dimension = len(offset) + return [ + Transform( + transformType=TransformType( + transformParameterization=TransformParameterizations.Translation, + parametersValueType=FloatTypes.Float64, + inputDimension=dimension, + outputDimension=dimension, + ), + numberOfFixedParameters=0, + numberOfParameters=dimension, + fixedParameters=np.empty((0,), dtype=np.float64), + parameters=np.asarray(offset, dtype=np.float64), + ) + ] + + +def _identity(dimension): + return _translation([0.0] * dimension) + + +def _image(dims, shape, scale, translation, orientations=None): + """A geometry-only NgffImage; the data is never meant to be computed.""" + return NgffImage( + data=da.zeros(tuple(shape[d] for d in dims), dtype=np.uint8, chunks=8), + dims=list(dims), + scale={d: float(scale[d]) for d in dims}, + translation={d: float(translation[d]) for d in dims}, + axes_orientations=orientations, + ) + + +def _affine(matrix, offset): + """An ITK-Wasm affine: row-major matrix then translation, centre at origin.""" + from itkwasm import ( + FloatTypes, + Transform, + TransformParameterizations, + TransformType, + ) + + dimension = len(offset) + parameters = np.concatenate([np.asarray(matrix, dtype=float).ravel(), offset]) + return [ + Transform( + transformType=TransformType( + transformParameterization=TransformParameterizations.Affine, + parametersValueType=FloatTypes.Float64, + inputDimension=dimension, + outputDimension=dimension, + ), + numberOfFixedParameters=dimension, + numberOfParameters=len(parameters), + fixedParameters=np.zeros(dimension, dtype=np.float64), + parameters=parameters, + ) + ] + + +def _oracle_region(matrix, offset, fixed, moving, spatial, padding): + """Recompute the region in NGFF order, independently of the pipeline. + + A linear map sends the fixed rectangle to a convex region, so sampling the + corners is exact. + """ + shape = [fixed.data.shape[list(fixed.dims).index(d)] for d in spatial] + fixed_scale = np.array([fixed.scale[d] for d in spatial]) + fixed_translation = np.array([fixed.translation[d] for d in spatial]) + moving_scale = np.array([moving.scale[d] for d in spatial]) + moving_translation = np.array([moving.translation[d] for d in spatial]) + + index_min = np.full(len(spatial), np.inf) + index_max = np.full(len(spatial), -np.inf) + for corner in itertools.product(*[(0, n - 1) for n in shape]): + point = fixed_translation + fixed_scale * np.array(corner, dtype=float) + moved = matrix @ point + offset + continuous_index = (moved - moving_translation) / moving_scale + index_min = np.minimum(index_min, continuous_index) + index_max = np.maximum(index_max, continuous_index) + + start = np.floor(index_min).astype(np.int64) - padding + end = np.ceil(index_max).astype(np.int64) + padding + return start, np.maximum(end - start + 1, 0) + + +def test_mismatched_spatial_dims_are_rejected(): + fixed = _image( + "zyx", + {"z": 4, "y": 4, "x": 4}, + {"z": 1, "y": 1, "x": 1}, + {"z": 0, "y": 0, "x": 0}, + ) + moving = _image("yx", {"y": 4, "x": 4}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + with pytest.raises(ValueError, match="they must match"): + itk_transform_resample_bounding_box(_identity(2), fixed, moving) + + +def test_pixel_buffers_are_never_computed(): + """The whole point: geometry in, region out, no pixels touched.""" + computed = [] + + def explode(block): + computed.append(True) + raise AssertionError("pixel data was materialized") + + poison = da.zeros((32, 32), dtype=np.uint8, chunks=8).map_blocks( + explode, dtype=np.uint8 + ) + fixed = NgffImage( + data=poison, + dims=["y", "x"], + scale={"y": 1.0, "x": 1.0}, + translation={"y": 0.0, "x": 0.0}, + ) + + # dask probes the block function once while building the graph above, so + # only what happens from here on counts. + computed.clear() + bounding_box = itk_transform_resample_bounding_box( + _identity(2), fixed, fixed, padding=1 + ) + + assert not computed + assert bounding_box.start_index == {"y": -1, "x": -1} + assert bounding_box.size == {"y": 34, "x": 34} + + +def test_non_spatial_axes_are_passed_through(): + shape = {"t": 3, "c": 2, "z": 4, "y": 8, "x": 8} + unit = dict.fromkeys(shape, 1) + zero = dict.fromkeys(shape, 0) + fixed = _image("tczyx", shape, unit, zero) + moving = _image("tczyx", {"t": 3, "c": 2, "z": 16, "y": 32, "x": 32}, unit, zero) + + bounding_box = itk_transform_resample_bounding_box( + _translation([3.0, 2.0, 1.0]), fixed, moving, padding=0 + ) + + assert bounding_box.dims == ("z", "y", "x") + assert bounding_box.start_index == {"z": 1, "y": 2, "x": 3} + slices = bounding_box.slices(moving.dims) + assert slices[0] == slice(None) # t + assert slices[1] == slice(None) # c + assert slices[2:] == (slice(1, 5), slice(2, 10), slice(3, 11)) + + +def test_region_outside_the_moving_image_is_empty(): + fixed = _image("yx", {"y": 4, "x": 4}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + moving = _image("yx", {"y": 32, "x": 32}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + + bounding_box = itk_transform_resample_bounding_box( + _translation([1000.0, 1000.0]), fixed, moving, padding=1 + ) + + assert bounding_box.is_empty + assert bounding_box.crop(moving) is None + + +def test_negative_start_index_is_clamped_not_wrapped(): + """A negative slice bound would wrap silently instead of raising.""" + fixed = _image("yx", {"y": 4, "x": 4}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + moving = _image("yx", {"y": 32, "x": 32}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + + bounding_box = itk_transform_resample_bounding_box( + _identity(2), fixed, moving, padding=2 + ) + + assert bounding_box.start_index == {"y": -2, "x": -2} + assert bounding_box.clamped() == {"y": (0, 6), "x": (0, 6)} + assert bounding_box.slices() == (slice(0, 6), slice(0, 6)) + assert not bounding_box.is_empty + + +def test_crop_is_lazy_and_shifts_the_translation(): + fixed = _image("yx", {"y": 64, "x": 64}, {"y": 1, "x": 1}, {"y": 512, "x": 1024}) + moving = _image("yx", {"y": 4096, "x": 4096}, {"y": 2, "x": 2}, {"y": 0, "x": 0}) + + bounding_box = itk_transform_resample_bounding_box( + _translation([0.0, 0.0]), fixed, moving, padding=1 + ) + cropped = bounding_box.crop(moving) + + assert isinstance(cropped.data, da.Array) + bounds = bounding_box.clamped() + assert cropped.data.shape == tuple( + stop - start for start, stop in (bounds["y"], bounds["x"]) + ) + for dim in ("y", "x"): + expected = moving.translation[dim] + bounds[dim][0] * moving.scale[dim] + assert np.isclose(cropped.translation[dim], expected) + # Reads a small corner rather than the whole moving image. + assert np.prod(cropped.data.shape) < 0.01 * np.prod(moving.data.shape) + + +def test_crop_preserves_orientation_and_scale(): + moving = _image( + "zyx", + {"z": 32, "y": 32, "x": 32}, + {"z": 2.0, "y": 1.0, "x": 0.5}, + {"z": 1.0, "y": 2.0, "x": 3.0}, + RAS, + ) + fixed = _image( + "zyx", + {"z": 4, "y": 4, "x": 4}, + {"z": 2.0, "y": 1.0, "x": 0.5}, + {"z": 1.0, "y": 2.0, "x": 3.0}, + RAS, + ) + bounding_box = itk_transform_resample_bounding_box(_identity(3), fixed, moving) + cropped = bounding_box.crop(moving) + + assert cropped.scale == moving.scale + assert cropped.axes_orientations == RAS + assert list(cropped.dims) == list(moving.dims) + + +@pytest.mark.parametrize("orientations", [None, RAS]) +def test_metadata_only_geometry_matches_ngff_image_to_itk_image(orientations): + """The ITK path must land in the same physical space Elastix registered in.""" + image = _image( + "zyx", + {"z": 4, "y": 8, "x": 16}, + {"z": 3.0, "y": 2.0, "x": 1.0}, + {"z": 30.0, "y": 20.0, "x": 10.0}, + orientations, + ) + reference = ngff_image_to_itk_image(image, wasm=True) + + itk_dims = ["x", "y", "z"] + geometry = _metadata_only_itk_image( + image, itk_dims, _itk_direction(image, itk_dims) + ) + + assert list(geometry.origin) == list(reference.origin) + assert list(geometry.spacing) == list(reference.spacing) + assert list(geometry.size) == list(reference.size) + assert np.allclose(np.asarray(geometry.direction), np.asarray(reference.direction)) + assert geometry.data.size == 0 + + +def test_ras_orientation_yields_a_non_identity_direction(): + image = _image( + "zyx", + {"z": 4, "y": 8, "x": 16}, + {"z": 1, "y": 1, "x": 1}, + {"z": 0, "y": 0, "x": 0}, + RAS, + ) + direction = _itk_direction(image, ["x", "y", "z"]) + assert np.allclose(direction, np.diag([-1.0, -1.0, 1.0])) + + +def test_itk_translation_transform_is_accepted(): + itk = pytest.importorskip("itk") + + fixed = _image("yx", {"y": 16, "x": 16}, {"y": 2, "x": 2}, {"y": 20, "x": 10}) + moving = _image("yx", {"y": 64, "x": 64}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + transform = itk.TranslationTransform[itk.D, 2].New() + transform.SetOffset([10.0, 5.0]) # ITK order (x, y) + + bounding_box = itk_transform_resample_bounding_box( + transform, fixed, moving, padding=1 + ) + + assert bounding_box.start_index == {"y": 24, "x": 19} + assert bounding_box.size == {"y": 33, "x": 33} + + +def test_itk_composite_transform_is_accepted(): + """An ITK composite applies its last-added transform first.""" + itk = pytest.importorskip("itk") + + translation = itk.TranslationTransform[itk.D, 2].New() + translation.SetOffset([10.0, 0.0]) + scaling = itk.AffineTransform[itk.D, 2].New() + scaling.SetMatrix(itk.matrix_from_array(np.array([[2.0, 0.0], [0.0, 2.0]]))) + scaling.SetTranslation([0.0, 0.0]) + scaling.SetCenter([0.0, 0.0]) + composite = itk.CompositeTransform[itk.D, 2].New() + composite.AddTransform(translation) + composite.AddTransform(scaling) + + fixed = _image("yx", {"y": 16, "x": 16}, {"y": 2, "x": 2}, {"y": 20, "x": 10}) + moving = _image("yx", {"y": 512, "x": 512}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + + bounding_box = itk_transform_resample_bounding_box( + composite, fixed, moving, padding=0 + ) + + # Cross-checked against the composite's own point mapping. + low = composite.TransformPoint([10.0, 20.0]) + high = composite.TransformPoint([40.0, 50.0]) + assert np.isclose(bounding_box.corners_min["x"], low[0]) + assert np.isclose(bounding_box.corners_min["y"], low[1]) + assert np.isclose(bounding_box.corners_max["x"], high[0]) + assert np.isclose(bounding_box.corners_max["y"], high[1]) + + +def test_index_range_overflow_is_reported_not_silently_empty(): + """A grid too large for the pipeline's index space must not read as empty. + + The region is computed in 32-bit index space while the corners come back as + doubles, so a fixed/moving scale mismatch of this size wraps the integer + region. Reporting "no overlap" for a grid that covers the whole moving + image would be a silent, wrong answer. + """ + fixed = _image("yx", {"y": 16, "x": 16}, {"y": 1e9, "x": 1e9}, {"y": 0, "x": 0}) + moving = _image("yx", {"y": 64, "x": 64}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + + with pytest.raises(ValueError, match="does not contain the transformed grid"): + itk_transform_resample_bounding_box(_identity(2), fixed, moving, padding=1) + + +@pytest.mark.parametrize("bad", [-1, 1.5, float("nan"), float("inf")]) +def test_padding_that_is_not_a_non_negative_integer_is_rejected(bad): + fixed = _image("yx", {"y": 4, "x": 4}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + with pytest.raises(ValueError, match="padding must be a non-negative integer"): + itk_transform_resample_bounding_box(_identity(2), fixed, fixed, padding=bad) + + +def test_unsupported_spatial_dimensionality_is_rejected(): + fixed = _image("x", {"x": 8}, {"x": 1}, {"x": 0}) + with pytest.raises(ValueError, match="only 2 and 3 are supported"): + itk_transform_resample_bounding_box(_identity(2), fixed, fixed) + + +@pytest.mark.parametrize("bad", [float("nan"), float("inf")]) +def test_non_finite_geometry_is_rejected(bad): + fixed = _image("yx", {"y": 4, "x": 4}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + fixed.scale["x"] = bad + with pytest.raises(ValueError, match="must be finite"): + itk_transform_resample_bounding_box(_identity(2), fixed, fixed) + + +def test_missing_scale_entry_is_rejected(): + fixed = _image("yx", {"y": 4, "x": 4}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + del fixed.scale["x"] + with pytest.raises(ValueError, match="no entry for dimension 'x'"): + itk_transform_resample_bounding_box(_identity(2), fixed, fixed) + + +def test_zero_scale_is_rejected(): + fixed = _image("yx", {"y": 4, "x": 4}, {"y": 1, "x": 0}, {"y": 0, "x": 0}) + with pytest.raises(ValueError, match="scale for dimension 'x' is zero"): + itk_transform_resample_bounding_box(_identity(2), fixed, fixed) + + +def test_degenerate_fixed_grid_yields_an_empty_region(): + fixed = _image("yx", {"y": 0, "x": 4}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + moving = _image("yx", {"y": 32, "x": 32}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + + bounding_box = itk_transform_resample_bounding_box(_identity(2), fixed, moving) + + assert bounding_box.is_empty + assert bounding_box.crop(moving) is None + + +def _constant_displacement_field(itk, shift, size=8, spacing=8.0): + """A displacement field that shifts every point by ``shift`` (ITK order).""" + field = itk.Image[itk.Vector[itk.D, 2], 2].New() + region = itk.ImageRegion[2]() + extent = itk.Size[2]() + extent[0], extent[1] = size, size + region.SetSize(extent) + field.SetRegions(region) + field.SetSpacing([spacing, spacing]) + field.SetOrigin([0.0, 0.0]) + field.Allocate() + offset = itk.Vector[itk.D, 2]() + offset[0], offset[1] = shift + field.FillBuffer(offset) + + transform = itk.DisplacementFieldTransform[itk.D, 2].New() + transform.SetDisplacementField(field) + return transform + + +def test_non_linear_displacement_field_is_supported(): + """Non-linear transforms are handled: the pipeline walks the grid boundary. + + This is what deformable registration needs, so it must not be refused + along with the transforms that genuinely cannot work. + """ + itk = pytest.importorskip("itk") + + transform = _constant_displacement_field(itk, (5.0, 3.0)) + assert not transform.IsLinear() + + fixed = _image("yx", {"y": 32, "x": 32}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + moving = _image("yx", {"y": 256, "x": 256}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + + bounding_box = itk_transform_resample_bounding_box( + transform, fixed, moving, padding=1 + ) + + # A constant field shifts ITK (x, y) by (5, 3), so NGFF (y, x) by (3, 5). + # The fixed grid spans 0..31 on both axes. + assert bounding_box.corners_min == {"y": 3.0, "x": 5.0} + assert bounding_box.corners_max == {"y": 34.0, "x": 36.0} + assert bounding_box.start_index == {"y": 2, "x": 4} + assert bounding_box.size == {"y": 34, "x": 34} + + +def test_unsupported_transform_type_is_rejected(): + fixed = _image("yx", {"y": 4, "x": 4}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + with pytest.raises(TypeError, match="unsupported transform type"): + itk_transform_resample_bounding_box("not a transform", fixed, fixed) + + +def test_asymmetric_three_dimensional_affine_matches_oracle(): + """An asymmetric sheared affine makes any axis-order slip visible.""" + spatial = ("z", "y", "x") + fixed = _image( + spatial, + {"z": 4, "y": 8, "x": 16}, + {"z": 3.0, "y": 2.0, "x": 1.0}, + {"z": 30.0, "y": 20.0, "x": 10.0}, + ) + moving = _image( + spatial, + {"z": 64, "y": 128, "x": 256}, + {"z": 1.5, "y": 0.5, "x": 0.25}, + {"z": -5.0, "y": 7.0, "x": 3.0}, + ) + # Stated in NGFF order for the oracle; the transform is built in ITK order, + # so both the rows and the columns reverse. + matrix = np.array([[1.0, 0.2, 0.0], [0.0, 2.0, 0.3], [0.5, 0.0, 1.0]]) + offset = np.array([4.0, -6.0, 11.0]) + reversal = np.eye(3)[::-1] + + bounding_box = itk_transform_resample_bounding_box( + _affine(reversal @ matrix @ reversal, reversal @ offset), + fixed, + moving, + padding=2, + ) + + expected_start, expected_size = _oracle_region( + matrix, offset, fixed, moving, spatial, 2 + ) + assert [bounding_box.start_index[d] for d in spatial] == expected_start.tolist() + assert [bounding_box.size[d] for d in spatial] == expected_size.tolist() + + +def _stub_result(padded_start_index, padded_size): + corners = {"min": [0.0, 0.0], "max": [1.0, 1.0]} + return { + "paddedStartIndex": padded_start_index, + "paddedSize": padded_size, + "corners": corners, + "paddedCorners": {"min": [0.0, 0.0], "max": [1.0, 1.0]}, + } + + +@pytest.mark.parametrize( + ("padded_start_index", "padded_size", "match"), + [ + ([0], [4, 4], "1 paddedStartIndex values for 2 dimensions"), + ([0, 0], [4], "1 paddedSize values for 2 dimensions"), + ([0, float("nan")], [4, 4], "paddedStartIndex for dimension 'y'"), + ([0, 0], [4, float("inf")], "paddedSize for dimension 'y'"), + ], +) +def test_unusable_pipeline_index_arrays_are_rejected( + monkeypatch, padded_start_index, padded_size, match +): + import itkwasm_downsample + + monkeypatch.setattr( + itkwasm_downsample, + "resample_bounding_box", + lambda *args, **kwargs: _stub_result(padded_start_index, padded_size), + ) + + fixed = _image("yx", {"y": 4, "x": 4}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + with pytest.raises(ValueError, match=match): + itk_transform_resample_bounding_box(_identity(2), fixed, fixed) From 1a89a63fc13f402025f887151131cf7b1b417f78 Mon Sep 17 00:00:00 2001 From: Valentin Boussot Date: Thu, 6 Aug 2026 15:10:52 +0200 Subject: [PATCH 2/7] feat(ts): add itkTransformResampleBoundingBox for out-of-core resampling TypeScript counterpart of the Python implementation, following the existing node/browser split. ResampleBoundingBox exposes selection() for a zarrita selection and croppedTranslation() for the shifted origin. --- ts/scripts/build_npm.ts | 2 +- ts/src/browser-mod.ts | 5 + ...transform_resample_bounding_box-browser.ts | 34 ++ ...tk_transform_resample_bounding_box-node.ts | 47 ++ ..._transform_resample_bounding_box-shared.ts | 434 ++++++++++++++ .../io/itk_transform_resample_bounding_box.ts | 21 + ts/src/mod.ts | 1 + ...tk_transform_resample_bounding_box_test.ts | 555 ++++++++++++++++++ 8 files changed, 1098 insertions(+), 1 deletion(-) create mode 100644 ts/src/io/itk_transform_resample_bounding_box-browser.ts create mode 100644 ts/src/io/itk_transform_resample_bounding_box-node.ts create mode 100644 ts/src/io/itk_transform_resample_bounding_box-shared.ts create mode 100644 ts/src/io/itk_transform_resample_bounding_box.ts create mode 100644 ts/test/itk_transform_resample_bounding_box_test.ts diff --git a/ts/scripts/build_npm.ts b/ts/scripts/build_npm.ts index 62f0e1ff..de84e249 100644 --- a/ts/scripts/build_npm.ts +++ b/ts/scripts/build_npm.ts @@ -197,7 +197,7 @@ async function createPackageJson(): Promise { dependencies: { "@fideus-labs/fizarrita": "^1.3.0", "@fideus-labs/worker-pool": "^1.0.0", - "@itk-wasm/downsample": "^1.8.1", + "@itk-wasm/downsample": "^2.0.0", "itk-wasm": "^1.0.0-b.196", "@zarrita/storage": "^0.1.4", zod: "^4.0.2", diff --git a/ts/src/browser-mod.ts b/ts/src/browser-mod.ts index f96dacb7..f5fb49c0 100644 --- a/ts/src/browser-mod.ts +++ b/ts/src/browser-mod.ts @@ -23,6 +23,11 @@ export { itkImageToNgffImage, type ItkImageToNgffImageOptions, } from "./io/itk_image_to_ngff_image.ts"; +export { itkTransformResampleBoundingBox } from "./io/itk_transform_resample_bounding_box-browser.ts"; +export { + type ItkTransformResampleBoundingBoxOptions, + ResampleBoundingBox, +} from "./io/itk_transform_resample_bounding_box-shared.ts"; export { dataTypeToComponentType, ngffImageToItkImage, diff --git a/ts/src/io/itk_transform_resample_bounding_box-browser.ts b/ts/src/io/itk_transform_resample_bounding_box-browser.ts new file mode 100644 index 00000000..e4744db7 --- /dev/null +++ b/ts/src/io/itk_transform_resample_bounding_box-browser.ts @@ -0,0 +1,34 @@ +// SPDX-FileCopyrightText: Copyright (c) Fideus Labs LLC +// SPDX-License-Identifier: MIT + +/** Browser implementation: dispatches through the package's web worker pool. */ + +import { resampleBoundingBox } from "@itk-wasm/downsample"; +import type { TransformList } from "itk-wasm"; +import type { NgffImage } from "../types/ngff_image.ts"; +import { + type ItkTransformResampleBoundingBoxOptions, + type ResampleBoundingBox, + resampleBoundingBoxShared, +} from "./itk_transform_resample_bounding_box-shared.ts"; + +/** + * Compute the moving-image region needed to resample a fixed image grid. + * + * Browser counterpart of the Node implementation; see + * `itk_transform_resample_bounding_box-node.ts` for the full description. + */ +export function itkTransformResampleBoundingBox( + transform: TransformList, + fixed: NgffImage, + moving: NgffImage, + options: ItkTransformResampleBoundingBoxOptions = {}, +): Promise { + return resampleBoundingBoxShared( + resampleBoundingBox, + transform, + fixed, + moving, + options, + ); +} diff --git a/ts/src/io/itk_transform_resample_bounding_box-node.ts b/ts/src/io/itk_transform_resample_bounding_box-node.ts new file mode 100644 index 00000000..6065fac0 --- /dev/null +++ b/ts/src/io/itk_transform_resample_bounding_box-node.ts @@ -0,0 +1,47 @@ +// SPDX-FileCopyrightText: Copyright (c) Fideus Labs LLC +// SPDX-License-Identifier: MIT + +/** Node/Deno implementation: native WASM, no web worker. */ + +import { resampleBoundingBoxNode } from "@itk-wasm/downsample"; +import type { TransformList } from "itk-wasm"; +import type { NgffImage } from "../types/ngff_image.ts"; +import { + type ItkTransformResampleBoundingBoxOptions, + type ResampleBoundingBox, + resampleBoundingBoxShared, +} from "./itk_transform_resample_bounding_box-shared.ts"; + +/** + * Compute the moving-image region needed to resample a fixed image grid. + * + * The region is derived from image geometry alone -- the pixel buffers of + * `fixed` and `moving` are never read. That is the point: describe two images + * and a transform with a few numbers, learn exactly which block of the moving + * image a resample will touch, and only then move pixels. + * + * The transform acts on ITK physical space, so the geometry is built the way + * {@link ngffImageToItkImage} builds it, including the direction matrix derived + * from RFC-4 anatomical orientation. It maps *fixed* points into *moving* + * space. + * + * @param transform An ITK-Wasm `TransformList`. + * @param fixed The image whose grid is resampled. Geometry only. + * @param moving The image to be sampled. Geometry only. + * @param options Padding options. + * @returns The region, keyed by dimension name in Zarr order. + */ +export function itkTransformResampleBoundingBox( + transform: TransformList, + fixed: NgffImage, + moving: NgffImage, + options: ItkTransformResampleBoundingBoxOptions = {}, +): Promise { + return resampleBoundingBoxShared( + resampleBoundingBoxNode, + transform, + fixed, + moving, + options, + ); +} diff --git a/ts/src/io/itk_transform_resample_bounding_box-shared.ts b/ts/src/io/itk_transform_resample_bounding_box-shared.ts new file mode 100644 index 00000000..c1371ed4 --- /dev/null +++ b/ts/src/io/itk_transform_resample_bounding_box-shared.ts @@ -0,0 +1,434 @@ +// SPDX-FileCopyrightText: Copyright (c) Fideus Labs LLC +// SPDX-License-Identifier: MIT + +/** + * Find the region of a moving image needed to resample a fixed image grid. + * + * Environment-agnostic core. The node and browser modules supply the + * `resampleBoundingBox` pipeline; everything else lives here. + */ + +import * as zarr from "zarrita"; +import type { Image, TransformList } from "itk-wasm"; +import { NgffImage } from "../types/ngff_image.ts"; +import { anatomicalOrientationToItkDirection } from "../types/rfc4.ts"; + +const SPATIAL_DIMS = ["x", "y", "z"]; + +/** The shape of the JSON the `resample-bounding-box` pipeline emits. */ +interface RawBoundingBox { + paddedStartIndex: number[]; + paddedSize: number[]; + paddedCorners: { min: number[]; max: number[] }; + corners: { min: number[]; max: number[] }; +} + +/** Options for {@link itkTransformResampleBoundingBox}. */ +export interface ItkTransformResampleBoundingBoxOptions { + /** + * Pixels of padding added per side. The default of 1 covers linear + * interpolation, which reads one neighbor beyond the continuous index + * bound. Use 0 for the tight region, or more for wider kernels. + */ + padding?: number; +} + +/** + * The region of a moving image needed to resample a fixed image grid. + * + * Every record is keyed by dimension name, so callers never have to track + * whether a given array is in Zarr or ITK axis order. + */ +export class ResampleBoundingBox { + /** Spatial dimension names, in the moving image's (Zarr) order. */ + readonly dims: string[]; + /** + * Start of the region in the moving image's index space. May be negative + * when the transformed grid extends past the moving image origin. + */ + readonly startIndex: Record; + /** Size of the region in moving-image pixels, before clamping. */ + readonly size: Record; + /** + * Tight extent of the transformed fixed grid, in moving-image physical + * coordinates. Independent of `padding`. + */ + readonly cornersMin: Record; + readonly cornersMax: Record; + /** Physical coordinates of the padded region's first and last pixel. */ + readonly paddedCornersMin: Record; + readonly paddedCornersMax: Record; + /** Shape of the moving image, used to clamp the region into bounds. */ + readonly movingShape: Record; + + constructor(options: { + dims: string[]; + startIndex: Record; + size: Record; + cornersMin: Record; + cornersMax: Record; + paddedCornersMin: Record; + paddedCornersMax: Record; + movingShape: Record; + }) { + this.dims = [...options.dims]; + this.startIndex = { ...options.startIndex }; + this.size = { ...options.size }; + this.cornersMin = { ...options.cornersMin }; + this.cornersMax = { ...options.cornersMax }; + this.paddedCornersMin = { ...options.paddedCornersMin }; + this.paddedCornersMax = { ...options.paddedCornersMax }; + this.movingShape = { ...options.movingShape }; + } + + /** + * The region intersected with the moving image bounds. + * + * @returns `{dim: [start, stop]}` with `0 <= start <= stop <= shape`. + */ + clamped(): Record { + const bounds: Record = {}; + for (const dim of this.dims) { + const extent = this.movingShape[dim]; + const start = Math.max(0, Math.min(this.startIndex[dim], extent)); + const stop = Math.max( + start, + Math.min(this.startIndex[dim] + this.size[dim], extent), + ); + bounds[dim] = [start, stop]; + } + return bounds; + } + + /** Whether the region does not overlap the moving image at all. */ + get isEmpty(): boolean { + return Object.values(this.clamped()).some(([start, stop]) => stop <= start); + } + + /** + * A zarrita selection for the region, clamped to the moving image bounds. + * + * Negative start indices are clamped rather than passed through, since a + * negative bound would silently wrap around instead of raising. + * + * @param dims Dimension order of the array to be read. Defaults to the + * spatial dims. Dimensions outside the region (`t`, `c`) select everything. + */ + selection(dims?: string[]): (zarr.Slice | null)[] { + const order = dims ?? this.dims; + const bounds = this.clamped(); + return order.map((dim) => + dim in bounds ? zarr.slice(bounds[dim][0], bounds[dim][1]) : null + ); + } + + /** + * The translation the cropped region would have, per dimension. + * + * Reading `selection()` yields an array whose origin has moved; this is the + * matching `translation` for an {@link NgffImage} built from it. + */ + croppedTranslation(moving: NgffImage): Record { + const translation = { ...moving.translation }; + for (const [dim, [start]] of Object.entries(this.clamped())) { + translation[dim] = moving.translation[dim] + start * moving.scale[dim]; + } + return translation; + } +} + +function spatialDims(image: NgffImage): string[] { + return image.dims.filter((dim) => SPATIAL_DIMS.includes(dim)); +} + +/** + * Require a finite scale and translation for every spatial axis. + * + * A missing or non-finite entry would otherwise reach the pipeline and come + * back as a plausible-looking region computed from the wrong geometry. + */ +function checkGeometry( + label: string, + image: NgffImage, + spatial: string[], +): void { + for ( + const [name, values] of [ + ["scale", image.scale], + ["translation", image.translation], + ] as const + ) { + for (const dim of spatial) { + const value = values[dim]; + if (value === undefined) { + throw new Error( + `${label} image ${name} has no entry for dimension '${dim}'`, + ); + } + if (!Number.isFinite(value)) { + throw new Error( + `${label} image ${name} for dimension '${dim}' is ${value}; ` + + `it must be finite`, + ); + } + if (name === "scale" && value === 0) { + throw new Error( + `${label} image scale for dimension '${dim}' is zero`, + ); + } + } + } +} + +/** + * Reject a region that does not contain the points it was derived from. + * + * The pipeline computes the integer region in 32-bit index space while + * reporting the corners as doubles, so a grid whose extent exceeds that range + * comes back as a wrapped -- and typically empty -- region. Padding is + * non-negative, so a correct region always spans its own tight corners; + * checking that catches the wrap without assuming an axis direction. + */ +// The returned records are built from these two arrays alone, so a short or +// non-finite pipeline result would otherwise become an undefined bound that +// only surfaces once a caller slices with it. +function checkIndexArrays(raw: RawBoundingBox, itkDims: string[]): void { + const arrays: [string, number[]][] = [ + ["paddedStartIndex", raw.paddedStartIndex], + ["paddedSize", raw.paddedSize], + ]; + for (const [name, values] of arrays) { + if (values?.length !== itkDims.length) { + throw new Error( + `the pipeline reported ${values?.length ?? 0} ${name} values for ` + + `${itkDims.length} dimensions`, + ); + } + itkDims.forEach((dim, axis) => { + if (!Number.isFinite(values[axis])) { + throw new Error( + `the ${name} for dimension '${dim}' is not finite; check the ` + + `transform and the image scale and translation`, + ); + } + }); + } +} + +function checkRegionContainsCorners( + raw: RawBoundingBox, + itkDims: string[], +): void { + itkDims.forEach((dim, axis) => { + const values = [ + raw.corners.min[axis], + raw.corners.max[axis], + raw.paddedCorners.min[axis], + raw.paddedCorners.max[axis], + ]; + if (values.some((value) => value === null || !Number.isFinite(value))) { + throw new Error( + `the resample region for dimension '${dim}' is not finite; check ` + + `the transform and the image scale and translation`, + ); + } + // paddedCorners hold the start and end index positions, which invert when + // an axis direction is negative, so compare against the span. + const low = Math.min( + raw.paddedCorners.min[axis], + raw.paddedCorners.max[axis], + ); + const high = Math.max( + raw.paddedCorners.min[axis], + raw.paddedCorners.max[axis], + ); + const tightLow = Math.min(raw.corners.min[axis], raw.corners.max[axis]); + const tightHigh = Math.max(raw.corners.min[axis], raw.corners.max[axis]); + if (tightLow < low || tightHigh > high) { + throw new Error( + `the resample region for dimension '${dim}' does not contain the ` + + `transformed grid (${tightLow} to ${tightHigh} lies outside ` + + `${low} to ${high}). The region spans more than the index range ` + + `the pipeline can represent; check for a scale mismatch between ` + + `the fixed and moving images.`, + ); + } + }); +} + +function identityDirection(dimension: number): Float64Array { + const direction = new Float64Array(dimension * dimension); + for (let i = 0; i < dimension; i++) direction[i * dimension + i] = 1.0; + return direction; +} + +/** + * Direction matrix from RFC-4 orientation, matching `ngffImageToItkImage`. + * + * All-or-nothing: unless every spatial axis carries an orientation that maps + * onto an LPS axis, the direction falls back to identity. + */ +export function itkDirection( + image: NgffImage, + itkDims: string[], +): Float64Array { + const dimension = itkDims.length; + const direction = identityDirection(dimension); + + const orientations = image.axesOrientations; + if (!orientations) return direction; + + const columns: number[][] = []; + for (const dim of itkDims) { + const orientation = orientations[dim]; + if (orientation === undefined) return direction; + const column = anatomicalOrientationToItkDirection(orientation.value); + if (column === undefined) return direction; + columns.push(column); + } + + for (let col = 0; col < dimension; col++) { + for (let row = 0; row < dimension; row++) { + direction[row * dimension + col] = columns[col][row]; + } + } + return direction; +} + +/** + * Build an ITK-Wasm image carrying geometry only, with an empty buffer. + * + * The zarr array is never read -- only its `shape` is consulted. This is what + * lets the region be computed against images whose pixels are remote or + * enormous. + */ +export function metadataOnlyItkImage( + image: NgffImage, + itkDims: string[], + direction: Float64Array, +): Image { + const itkImage: Image = { + imageType: { + dimension: itkDims.length, + componentType: "uint8", + pixelType: "Scalar", + components: 1, + }, + name: image.name, + origin: itkDims.map((dim) => image.translation[dim]), + spacing: itkDims.map((dim) => image.scale[dim]), + direction, + size: itkDims.map((dim) => image.data.shape[image.dims.indexOf(dim)]), + metadata: new Map(), + data: new Uint8Array(0), + }; + return itkImage; +} + +/** + * Compute the moving-image region needed to resample a fixed image grid. + * + * @param pipeline The environment's `resampleBoundingBox` implementation. + */ +export async function resampleBoundingBoxShared( + pipeline: ( + transform: TransformList, + fixed: Image, + moving: Image, + options: { padding?: number }, + ) => Promise<{ boundingBox: unknown }>, + transform: TransformList, + fixed: NgffImage, + moving: NgffImage, + options: ItkTransformResampleBoundingBoxOptions = {}, +): Promise { + const padding = options.padding ?? 1; + if (!Number.isInteger(padding) || padding < 0) { + throw new Error(`padding must be a non-negative integer, got ${padding}`); + } + + const fixedSpatial = spatialDims(fixed); + const movingSpatial = spatialDims(moving); + if (fixedSpatial.join(",") !== movingSpatial.join(",")) { + throw new Error( + `fixed image has spatial dims [${fixedSpatial.join(", ")}] but moving ` + + `image has [${movingSpatial.join(", ")}]; they must match`, + ); + } + if (fixedSpatial.length < 2) { + throw new Error( + `images have ${fixedSpatial.length} spatial dims ` + + `[${fixedSpatial.join(", ")}]; only 2 and 3 are supported`, + ); + } + + checkGeometry("fixed", fixed, fixedSpatial); + checkGeometry("moving", moving, fixedSpatial); + + // ITK orders points fastest-axis-first, the reverse of the Zarr order. + const itkDims = [...fixedSpatial].reverse(); + + const movingShape: Record = {}; + for (const dim of movingSpatial) { + movingShape[dim] = moving.data.shape[moving.dims.indexOf(dim)]; + } + const fixedDegenerate = fixedSpatial.some( + (dim) => fixed.data.shape[fixed.dims.indexOf(dim)] === 0, + ); + if (fixedDegenerate) { + // A grid with a zero-length axis has no samples to resample, so there is + // nothing to fetch. + const zeros: Record = {}; + for (const dim of fixedSpatial) zeros[dim] = 0; + return new ResampleBoundingBox({ + dims: fixedSpatial, + startIndex: { ...zeros }, + size: { ...zeros }, + cornersMin: { ...zeros }, + cornersMax: { ...zeros }, + paddedCornersMin: { ...zeros }, + paddedCornersMax: { ...zeros }, + movingShape, + }); + } + + // An ITK transform list acts on ITK physical space, so the geometry is built + // the way ngffImageToItkImage builds it, direction included. + const transformList = transform; + const fixedDirection = itkDirection(fixed, itkDims); + const movingDirection = itkDirection(moving, itkDims); + + const { boundingBox } = await pipeline( + transformList, + metadataOnlyItkImage(fixed, itkDims, fixedDirection), + metadataOnlyItkImage(moving, itkDims, movingDirection), + { padding }, + ); + const raw = boundingBox as RawBoundingBox; + checkIndexArrays(raw, itkDims); + checkRegionContainsCorners(raw, itkDims); + + // The pipeline reports arrays fastest-axis-first; key them by dimension + // name so no caller has to remember that, and order the records like `dims` + // so a serialized result reads in Zarr order. + const byDim = (values: number[]): Record => { + const indexed: Record = {}; + itkDims.forEach((dim, index) => { + indexed[dim] = values[index]; + }); + const record: Record = {}; + for (const dim of fixedSpatial) record[dim] = indexed[dim]; + return record; + }; + + return new ResampleBoundingBox({ + dims: fixedSpatial, + startIndex: byDim(raw.paddedStartIndex), + size: byDim(raw.paddedSize), + cornersMin: byDim(raw.corners.min), + cornersMax: byDim(raw.corners.max), + paddedCornersMin: byDim(raw.paddedCorners.min), + paddedCornersMax: byDim(raw.paddedCorners.max), + movingShape, + }); +} diff --git a/ts/src/io/itk_transform_resample_bounding_box.ts b/ts/src/io/itk_transform_resample_bounding_box.ts new file mode 100644 index 00000000..4d3b1c1d --- /dev/null +++ b/ts/src/io/itk_transform_resample_bounding_box.ts @@ -0,0 +1,21 @@ +// SPDX-FileCopyrightText: Copyright (c) Fideus Labs LLC +// SPDX-License-Identifier: MIT + +/** + * Resample bounding box support. + * + * This module provides conditional exports for browser and Node environments. + * The actual implementation is delegated to environment-specific modules: + * - itk_transform_resample_bounding_box-browser.ts: WebWorker-based functions + * - itk_transform_resample_bounding_box-node.ts: native WASM for Node/Deno + * + * For Deno runtime, we default to the node implementation. + * For browser bundlers, they should use conditional exports in package.json + * to resolve to the browser implementation. + */ + +export { itkTransformResampleBoundingBox } from "./itk_transform_resample_bounding_box-node.ts"; +export { + type ItkTransformResampleBoundingBoxOptions, + ResampleBoundingBox, +} from "./itk_transform_resample_bounding_box-shared.ts"; diff --git a/ts/src/mod.ts b/ts/src/mod.ts index 0cbcba25..c9310d26 100644 --- a/ts/src/mod.ts +++ b/ts/src/mod.ts @@ -5,6 +5,7 @@ export { config, setWorkerPoolSize } from "./config.ts"; export * from "./io/from_ngff_zarr.ts"; export * from "./io/hcs.ts"; export * from "./io/itk_image_to_ngff_image.ts"; +export * from "./io/itk_transform_resample_bounding_box.ts"; export * from "./io/ngff_image_to_itk_image.ts"; export type { MemoryStoreToZipOptions } from "./io/rfc9_zip.ts"; // RFC-9 exports diff --git a/ts/test/itk_transform_resample_bounding_box_test.ts b/ts/test/itk_transform_resample_bounding_box_test.ts new file mode 100644 index 00000000..460d7c38 --- /dev/null +++ b/ts/test/itk_transform_resample_bounding_box_test.ts @@ -0,0 +1,555 @@ +// SPDX-FileCopyrightText: Copyright (c) Fideus Labs LLC +// SPDX-License-Identifier: MIT + +/** + * RFC-5 to ITK transform bridge and resample bounding box tests. + * + * Mirrors `py/test/test_itk_transform_resample_bounding_box.py`. The expected + * regions are not taken from the pipeline itself: they come either from the + * worked examples in the ITK-Wasm `resample-bounding-box` documentation, or + * from `oracleRegion` below, which recomputes the region from first principles + * in NGFF axis order. That keeps the tests honest about the two conventions + * that differ between RFC-5 and ITK: axis order and sequence composition + * order. + */ + +import { assertAlmostEquals, assertEquals, assertRejects } from "@std/assert"; +import * as zarr from "zarrita"; +import { itkTransformResampleBoundingBox, NgffImage } from "../src/mod.ts"; +import { resampleBoundingBoxShared } from "../src/io/itk_transform_resample_bounding_box-shared.ts"; +import { RAS } from "../src/types/rfc4.ts"; + +/** An ITK-Wasm translation, in ITK (fastest-axis-first) order. */ +// deno-lint-ignore no-explicit-any +function itkTranslation(offset: number[]): any { + return [{ + transformType: { + transformParameterization: "Translation", + parametersValueType: "float64", + inputDimension: offset.length, + outputDimension: offset.length, + }, + name: "TranslationTransform", + inputSpaceName: "", + outputSpaceName: "", + numberOfFixedParameters: 0, + numberOfParameters: offset.length, + fixedParameters: new Float64Array(0), + parameters: new Float64Array(offset), + metadata: new Map(), + }]; +} + +// deno-lint-ignore no-explicit-any +const identity = (dimension: number): any => + itkTranslation(new Array(dimension).fill(0)); + +/** An ITK-Wasm affine, row-major matrix then translation, centre at the origin. */ +// deno-lint-ignore no-explicit-any +function itkAffine(matrix: number[][], offset: number[]): any { + const dimension = offset.length; + return [{ + transformType: { + transformParameterization: "Affine", + parametersValueType: "float64", + inputDimension: dimension, + outputDimension: dimension, + }, + name: "AffineTransform", + inputSpaceName: "", + outputSpaceName: "", + numberOfFixedParameters: dimension, + numberOfParameters: dimension * dimension + dimension, + fixedParameters: new Float64Array(dimension), + parameters: new Float64Array([...matrix.flat(), ...offset]), + metadata: new Map(), + }]; +} + +/** Reverse a square matrix's rows and columns: NGFF order <-> ITK order. */ +function reversed(matrix: number[][]): number[][] { + const order = matrix.map((_, i) => matrix.length - 1 - i); + return order.map((r) => order.map((c) => matrix[r][c])); +} + +/** A geometry-only NgffImage; the data is never meant to be read. */ +async function geometryImage( + dims: string[], + shape: Record, + scale: Record, + translation: Record, + axesOrientations?: Record, +): Promise { + const store = new Map(); + const root = zarr.root(store); + const arrayShape = dims.map((dim) => shape[dim]); + const data = await zarr.create(root.resolve("data"), { + shape: arrayShape, + chunk_shape: arrayShape.map((n) => Math.min(n, 32)), + data_type: "uint8", + fill_value: 0, + }); + return new NgffImage({ + data, + dims, + scale, + translation, + name: "image", + axesUnits: undefined, + axesOrientations, + computedCallbacks: undefined, + }); +} + +/** Recompute the region in NGFF order, independently of the pipeline. */ +function oracleRegion( + matrix: number[][], + offset: number[], + fixedShape: number[], + fixedScale: number[], + fixedTranslation: number[], + movingScale: number[], + movingTranslation: number[], + padding: number, +): { start: number[]; size: number[] } { + const ndim = fixedShape.length; + const indexMin = new Array(ndim).fill(Infinity); + const indexMax = new Array(ndim).fill(-Infinity); + + // A linear map sends the fixed rectangle to a convex region, so sampling + // the corners is exact. + for (let mask = 0; mask < 1 << ndim; mask++) { + const corner = Array.from( + { length: ndim }, + (_, i) => (mask >> i) & 1 ? fixedShape[i] - 1 : 0, + ); + const point = corner.map((c, i) => fixedTranslation[i] + fixedScale[i] * c); + for (let row = 0; row < ndim; row++) { + let moved = offset[row]; + for (let col = 0; col < ndim; col++) { + moved += matrix[row][col] * point[col]; + } + const continuousIndex = (moved - movingTranslation[row]) / + movingScale[row]; + indexMin[row] = Math.min(indexMin[row], continuousIndex); + indexMax[row] = Math.max(indexMax[row], continuousIndex); + } + } + + const start = indexMin.map((v) => Math.floor(v) - padding); + const end = indexMax.map((v) => Math.ceil(v) + padding); + return { start, size: end.map((e, i) => Math.max(e - start[i] + 1, 0)) }; +} + +Deno.test("mismatched spatial dims are rejected", async () => { + const fixed = await geometryImage(["z", "y", "x"], { z: 4, y: 4, x: 4 }, { + z: 1, + y: 1, + x: 1, + }, { z: 0, y: 0, x: 0 }); + const moving = await geometryImage(["y", "x"], { y: 4, x: 4 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + + await assertRejects( + () => itkTransformResampleBoundingBox(identity(2), fixed, moving), + Error, + "they must match", + ); +}); + +Deno.test("pixel data is never read", async () => { + // A store that throws on any chunk read: if the pipeline touched pixels, + // this would reject rather than return a region. + const backing = new Map(); + const poison = { + get(key: string): Promise { + if (key.endsWith("zarr.json")) { + return Promise.resolve(backing.get(key)); + } + throw new Error("pixel data was read"); + }, + set(key: string, value: Uint8Array): Promise { + backing.set(key, value); + return Promise.resolve(); + }, + delete(key: string): Promise { + return Promise.resolve(backing.delete(key)); + }, + }; + const root = zarr.root(poison as never); + const data = await zarr.create(root.resolve("data"), { + shape: [32, 32], + chunk_shape: [8, 8], + data_type: "uint8", + fill_value: 0, + }); + const image = new NgffImage({ + data, + dims: ["y", "x"], + scale: { y: 1, x: 1 }, + translation: { y: 0, x: 0 }, + name: "image", + axesUnits: undefined, + computedCallbacks: undefined, + }); + + const boundingBox = await itkTransformResampleBoundingBox( + identity(2), + image, + image, + { padding: 1 }, + ); + + assertEquals(boundingBox.startIndex, { y: -1, x: -1 }); + assertEquals(boundingBox.size, { y: 34, x: 34 }); +}); + +Deno.test("non-spatial axes are passed through", async () => { + const dims = ["t", "c", "z", "y", "x"]; + const fixed = await geometryImage(dims, { t: 3, c: 2, z: 4, y: 8, x: 8 }, { + t: 1, + c: 1, + z: 1, + y: 1, + x: 1, + }, { t: 0, c: 0, z: 0, y: 0, x: 0 }); + const moving = await geometryImage( + dims, + { t: 3, c: 2, z: 16, y: 32, x: 32 }, + { + t: 1, + c: 1, + z: 1, + y: 1, + x: 1, + }, + { t: 0, c: 0, z: 0, y: 0, x: 0 }, + ); + + const boundingBox = await itkTransformResampleBoundingBox( + itkTranslation([3, 2, 1]), + fixed, + moving, + { padding: 0 }, + ); + + assertEquals(boundingBox.dims, ["z", "y", "x"]); + assertEquals(boundingBox.startIndex, { z: 1, y: 2, x: 3 }); + + const selection = boundingBox.selection(moving.dims); + assertEquals(selection[0], null); // t + assertEquals(selection[1], null); // c +}); + +Deno.test("a region outside the moving image is empty", async () => { + const fixed = await geometryImage(["y", "x"], { y: 4, x: 4 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + const moving = await geometryImage(["y", "x"], { y: 32, x: 32 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + + const boundingBox = await itkTransformResampleBoundingBox( + itkTranslation([1000, 1000]), + fixed, + moving, + { padding: 1 }, + ); + + assertEquals(boundingBox.isEmpty, true); +}); + +Deno.test("a negative start index is clamped, not wrapped", async () => { + const fixed = await geometryImage(["y", "x"], { y: 4, x: 4 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + const moving = await geometryImage(["y", "x"], { y: 32, x: 32 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + + const boundingBox = await itkTransformResampleBoundingBox( + identity(2), + fixed, + moving, + { padding: 2 }, + ); + + assertEquals(boundingBox.startIndex, { y: -2, x: -2 }); + assertEquals(boundingBox.clamped(), { y: [0, 6], x: [0, 6] }); + assertEquals(boundingBox.isEmpty, false); +}); + +Deno.test("the cropped translation shifts by start * scale", async () => { + const fixed = await geometryImage(["y", "x"], { y: 64, x: 64 }, { + y: 1, + x: 1, + }, { y: 512, x: 1024 }); + const moving = await geometryImage(["y", "x"], { y: 1024, x: 1024 }, { + y: 2, + x: 2, + }, { y: 0, x: 0 }); + + const boundingBox = await itkTransformResampleBoundingBox( + itkTranslation([0, 0]), + fixed, + moving, + { padding: 1 }, + ); + const translation = boundingBox.croppedTranslation(moving); + const bounds = boundingBox.clamped(); + + for (const dim of ["y", "x"]) { + assertAlmostEquals( + translation[dim], + moving.translation[dim] + bounds[dim][0] * moving.scale[dim], + ); + } +}); + +Deno.test("RAS orientation yields a non-identity direction", async () => { + const { itkDirection } = await import( + "../src/io/itk_transform_resample_bounding_box-shared.ts" + ); + const image = await geometryImage( + ["z", "y", "x"], + { z: 4, y: 8, x: 16 }, + { z: 1, y: 1, x: 1 }, + { z: 0, y: 0, x: 0 }, + RAS as never, + ); + const direction = itkDirection(image, ["x", "y", "z"]); + assertEquals(Array.from(direction), [-1, 0, 0, 0, -1, 0, 0, 0, 1]); +}); + +Deno.test("an index range overflow is reported, not silently empty", async () => { + // The region is computed in 32-bit index space while the corners come back + // as doubles, so a fixed/moving scale mismatch this large wraps the integer + // region. Reporting "no overlap" for a grid covering the whole moving image + // would be a silent, wrong answer. + const fixed = await geometryImage(["y", "x"], { y: 16, x: 16 }, { + y: 1e9, + x: 1e9, + }, { y: 0, x: 0 }); + const moving = await geometryImage(["y", "x"], { y: 64, x: 64 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + + await assertRejects( + () => + itkTransformResampleBoundingBox(identity(2), fixed, moving, { + padding: 1, + }), + Error, + "does not contain the transformed grid", + ); +}); + +Deno.test("padding that is not a non-negative integer is rejected", async () => { + const fixed = await geometryImage(["y", "x"], { y: 4, x: 4 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + for (const padding of [-1, 1.5, NaN, Infinity]) { + await assertRejects( + () => + itkTransformResampleBoundingBox(identity(2), fixed, fixed, { + padding, + }), + Error, + "padding must be a non-negative integer", + ); + } +}); + +Deno.test("unsupported spatial dimensionality is rejected", async () => { + const fixed = await geometryImage(["x"], { x: 8 }, { x: 1 }, { x: 0 }); + await assertRejects( + () => itkTransformResampleBoundingBox(identity(2), fixed, fixed), + Error, + "only 2 and 3 are supported", + ); +}); + +Deno.test("a missing scale entry is rejected rather than defaulted", async () => { + // Python raises here; TypeScript used to substitute spacing 1 and return a + // plausible but wrong region. + const fixed = await geometryImage(["y", "x"], { y: 4, x: 4 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + const partial = new NgffImage({ + data: fixed.data, + dims: ["y", "x"], + scale: { y: 2 } as Record, + translation: { y: 0, x: 0 }, + name: "image", + axesUnits: undefined, + computedCallbacks: undefined, + }); + + await assertRejects( + () => itkTransformResampleBoundingBox(identity(2), partial, fixed), + Error, + "no entry for dimension 'x'", + ); +}); + +Deno.test("a non-finite scale is rejected", async () => { + const fixed = await geometryImage(["y", "x"], { y: 4, x: 4 }, { + y: 1, + x: NaN, + }, { y: 0, x: 0 }); + await assertRejects( + () => itkTransformResampleBoundingBox(identity(2), fixed, fixed), + Error, + "must be finite", + ); +}); + +Deno.test("a zero scale is rejected", async () => { + const fixed = await geometryImage(["y", "x"], { y: 4, x: 4 }, { + y: 1, + x: 0, + }, { y: 0, x: 0 }); + await assertRejects( + () => itkTransformResampleBoundingBox(identity(2), fixed, fixed), + Error, + "scale for dimension 'x' is zero", + ); +}); + +Deno.test("a degenerate fixed grid yields an empty region", async () => { + const fixed = await geometryImage(["y", "x"], { y: 0, x: 4 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + const moving = await geometryImage(["y", "x"], { y: 32, x: 32 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + + const boundingBox = await itkTransformResampleBoundingBox( + identity(2), + fixed, + moving, + ); + + assertEquals(boundingBox.isEmpty, true); +}); + +Deno.test("an ITK-Wasm transform list is accepted", async () => { + const fixed = await geometryImage(["y", "x"], { y: 16, x: 16 }, { + y: 2, + x: 2, + }, { y: 20, x: 10 }); + const moving = await geometryImage(["y", "x"], { y: 64, x: 64 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + + // Stated directly in ITK order (x, y). + const transformList = [{ + transformType: { + transformParameterization: "Translation", + parametersValueType: "float64", + inputDimension: 2, + outputDimension: 2, + }, + name: "TranslationTransform", + inputSpaceName: "", + outputSpaceName: "", + numberOfFixedParameters: 0, + numberOfParameters: 2, + fixedParameters: new Float64Array(0), + parameters: new Float64Array([10, 5]), + metadata: new Map(), + }]; + + const boundingBox = await itkTransformResampleBoundingBox( + transformList as never, + fixed, + moving, + { padding: 1 }, + ); + + assertEquals(boundingBox.startIndex, { y: 24, x: 19 }); + assertEquals(boundingBox.size, { y: 33, x: 33 }); +}); + +Deno.test("an asymmetric 3D affine matches the oracle", () => { + // Stated in NGFF (z, y, x) order for the oracle; the transform itself is + // built in ITK order, so both row and column ordering reverse. + const matrix = [[1, 0.2, 0], [0, 2, 0.3], [0.5, 0, 1]]; + const offset = [4, -6, 11]; + return (async () => { + const dims = ["z", "y", "x"]; + const fixed = await geometryImage(dims, { z: 4, y: 8, x: 16 }, { + z: 3, + y: 2, + x: 1, + }, { z: 30, y: 20, x: 10 }); + const moving = await geometryImage(dims, { z: 64, y: 128, x: 256 }, { + z: 1.5, + y: 0.5, + x: 0.25, + }, { z: -5, y: 7, x: 3 }); + + const boundingBox = await itkTransformResampleBoundingBox( + itkAffine(reversed(matrix), [...offset].reverse()), + fixed, + moving, + { padding: 2 }, + ); + + const expected = oracleRegion( + matrix, + offset, + [4, 8, 16], + [3, 2, 1], + [30, 20, 10], + [1.5, 0.5, 0.25], + [-5, 7, 3], + 2, + ); + assertEquals(dims.map((d) => boundingBox.startIndex[d]), expected.start); + assertEquals(dims.map((d) => boundingBox.size[d]), expected.size); + })(); +}); + +Deno.test("unusable pipeline index arrays are rejected", async () => { + const fixed = await geometryImage(["y", "x"], { y: 4, x: 4 }, { + y: 1, + x: 1, + }, { y: 0, x: 0 }); + const corners = { min: [0, 0], max: [1, 1] }; + const cases: [number[], number[], string][] = [ + [[0], [4, 4], "1 paddedStartIndex values for 2 dimensions"], + [[0, 0], [4], "1 paddedSize values for 2 dimensions"], + [[0, NaN], [4, 4], "paddedStartIndex for dimension 'y'"], + [[0, 0], [4, Infinity], "paddedSize for dimension 'y'"], + ]; + + for (const [paddedStartIndex, paddedSize, message] of cases) { + const pipeline = () => + Promise.resolve({ + boundingBox: { + paddedStartIndex, + paddedSize, + corners, + paddedCorners: corners, + }, + }); + await assertRejects( + () => resampleBoundingBoxShared(pipeline, identity(2), fixed, fixed), + Error, + message, + ); + } +}); From c71e5e06b9b4c824a71d00be5d782dec37aa2a0c Mon Sep 17 00:00:00 2001 From: Valentin Boussot Date: Thu, 6 Aug 2026 15:10:52 +0200 Subject: [PATCH 3/7] docs: describe out-of-core resampling of a moving image --- docs/itk.md | 86 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) diff --git a/docs/itk.md b/docs/itk.md index 3bb5bf8d..9aae8611 100644 --- a/docs/itk.md +++ b/docs/itk.md @@ -45,3 +45,89 @@ Python dataclass like `NgffImage`. >>> # Back again >>> itk_wasm_image = nz.ngff_image_to_itk_image(ngff_image, wasm=True) ``` + +## Out-of-core resampling + +Resampling a fixed grid through a transform only ever reads moving-image +samples inside the transformed footprint of that grid. When the moving image is +large, remote, or chunked, materializing all of it to resample a small +overlapping region is wasteful. + +`itk_transform_resample_bounding_box` answers *which moving-image indices will +the resample actually read?* from image geometry alone. The pixel buffers are +never touched and the Dask graphs are never computed: + +```python +>>> import itk +>>> import ngff_zarr as nz +>>> +>>> # Any linear or deformable ITK transform, including the CompositeTransform +>>> # an Elastix registration returns. It maps fixed points into moving space. +>>> transform = registration_method.GetCombinedTransform() # doctest: +SKIP +>>> region = nz.itk_transform_resample_bounding_box( # doctest: +SKIP +... transform, fixed_block, moving) +>>> region.start_index # doctest: +SKIP +{'y': 11, 'x': -5} +``` + +The result is keyed by dimension name, so there is no ambiguity about axis +order -- the underlying pipeline reports arrays fastest-axis-first, the reverse +of the Zarr order. + +`region.crop(moving)` returns a lazily sliced `NgffImage` whose `translation` +has been shifted to match, ready to hand to `ngff_image_to_itk_image`: + +```python +>>> block = region.crop(moving) # doctest: +SKIP +>>> moving_itk = nz.ngff_image_to_itk_image(block, wasm=False) # doctest: +SKIP +``` + +Only that block's chunks are read. `crop` returns `None` when the transformed +grid does not overlap the moving image at all, so a tiling loop can skip it +instead of resampling nothing. Start indices may be negative when the grid +extends past the moving origin; `crop`, `slices` and `clamped` clamp into +bounds rather than letting a negative index wrap around. + +The image geometry is built the way `ngff_image_to_itk_image` builds it, +including the direction matrix derived from [RFC-4](./rfc4.md) anatomical +orientation, so the transform is applied in the space a registration produced +it in. + +Use `padding` to cover the interpolator's support. The default of `1` covers +linear interpolation, which reads one neighbor beyond the continuous index +bound; pass `0` for the tight region or a larger value for wider kernels. + +### Non-linear transforms + +Deformable registration is supported. For a linear transform the region is +derived from the transformed grid corners, which is exact because a linear map +sends a rectangle to a convex region. For a non-linear one that would +*under*-bound the region -- an interior edge point can map outside the hull of +the transformed corners -- so the whole grid boundary is walked instead. Cost is +proportional to the boundary, not the pixel count, and per block that boundary +is small. + +```{note} +`itk.BSplineTransform` currently aborts inside the `itkwasm-downsample` +pipeline. This is an upstream defect in how that pipeline reconstructs a +transform, not a limitation of the approach; every other parameterization +tested -- rigid, similarity, affine, versor, and displacement fields -- works. +``` + +## TypeScript + +The TypeScript package provides `itkTransformResampleBoundingBox`. It is async, +takes options as an object, and returns a `ResampleBoundingBox` whose +`selection()` yields a zarrita selection instead of Python slices: + +```typescript +import { itkTransformResampleBoundingBox, zarrGet } from "@fideus-labs/ngff-zarr"; + +const region = await itkTransformResampleBoundingBox(transform, fixed, moving, { + padding: 1, +}); + +if (!region.isEmpty) { + const block = await zarrGet(moving.data, region.selection(moving.dims)); +} +``` From d526e9550e30b24767037d00974ce59fbbb4baf1 Mon Sep 17 00:00:00 2001 From: Valentin Boussot Date: Fri, 7 Aug 2026 20:53:53 +0200 Subject: [PATCH 4/7] test(py): rebuild the brain_two_components baseline in canonical axis order Merging main brought #623, which normalizes generated axes to the spec order (t, c, z, y, x). The brain_two_components DASK_IMAGE_GAUSSIAN baseline in the v0.21.0 testing-data archive was generated before that change, with the component axis last, so the key sets no longer match. Regenerate that baseline against the merged code and pin the updated archive. --- py/test/_data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/py/test/_data.py b/py/test/_data.py index a10efea1..d0a9aae4 100644 --- a/py/test/_data.py +++ b/py/test/_data.py @@ -17,7 +17,7 @@ from zarr.storage import MemoryStore test_data_ipfs_cid = "bafybeifqibhcomn4u42aqrgvttyfteysbspvzez5sbezcqj5yylzzafpma" -test_data_sha256 = "e24fc764f562b68724665a9d36fafb64661b9cb433653d12597dd66ed9f28439" +test_data_sha256 = "525dfae8fe52df4a18dc19de97f018e667e161bc83dc0584144c16d872349705" test_dir = Path(__file__).resolve().parent extract_dir = "data" From 68daf97434198a7e64ff3c822a8952d0169fea84 Mon Sep 17 00:00:00 2001 From: Valentin Boussot Date: Mon, 10 Aug 2026 22:49:05 +0200 Subject: [PATCH 5/7] fix(ts): fall back to identity for 3D-only orientations on 2D images A superior/inferior orientation points along LPS z; truncating its column into a 2D direction matrix produced a singular matrix instead of the identity fallback. --- ts/src/io/itk_transform_resample_bounding_box-shared.ts | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/ts/src/io/itk_transform_resample_bounding_box-shared.ts b/ts/src/io/itk_transform_resample_bounding_box-shared.ts index c1371ed4..eaf73a71 100644 --- a/ts/src/io/itk_transform_resample_bounding_box-shared.ts +++ b/ts/src/io/itk_transform_resample_bounding_box-shared.ts @@ -284,6 +284,12 @@ export function itkDirection( if (orientation === undefined) return direction; const column = anatomicalOrientationToItkDirection(orientation.value); if (column === undefined) return direction; + // A column pointing outside the matrix's dimension (e.g. a + // superior/inferior orientation on a 2D image) would truncate to a + // singular matrix; keep the identity fallback instead. + if (column.slice(dimension).some((component) => component !== 0)) { + return direction; + } columns.push(column); } From 94813d428179c23d2f6d14655aff304287c7d0a2 Mon Sep 17 00:00:00 2001 From: Valentin Boussot Date: Mon, 10 Aug 2026 22:49:05 +0200 Subject: [PATCH 6/7] test(ts): type the axes orientations and assert the spatial selection Use AnatomicalOrientation for geometryImage's axesOrientations instead of casting RAS through never, and check the z/y/x selection entries in the non-spatial pass-through test. --- ...tk_transform_resample_bounding_box_test.ts | 32 +++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/ts/test/itk_transform_resample_bounding_box_test.ts b/ts/test/itk_transform_resample_bounding_box_test.ts index 460d7c38..e3ef2229 100644 --- a/ts/test/itk_transform_resample_bounding_box_test.ts +++ b/ts/test/itk_transform_resample_bounding_box_test.ts @@ -18,6 +18,7 @@ import * as zarr from "zarrita"; import { itkTransformResampleBoundingBox, NgffImage } from "../src/mod.ts"; import { resampleBoundingBoxShared } from "../src/io/itk_transform_resample_bounding_box-shared.ts"; import { RAS } from "../src/types/rfc4.ts"; +import type { AnatomicalOrientation } from "../src/types/rfc4.ts"; /** An ITK-Wasm translation, in ITK (fastest-axis-first) order. */ // deno-lint-ignore no-explicit-any @@ -78,7 +79,7 @@ async function geometryImage( shape: Record, scale: Record, translation: Record, - axesOrientations?: Record, + axesOrientations?: Record, ): Promise { const store = new Map(); const root = zarr.root(store); @@ -241,6 +242,9 @@ Deno.test("non-spatial axes are passed through", async () => { const selection = boundingBox.selection(moving.dims); assertEquals(selection[0], null); // t assertEquals(selection[1], null); // c + assertEquals(selection[2], zarr.slice(1, 5)); // z + assertEquals(selection[3], zarr.slice(2, 10)); // y + assertEquals(selection[4], zarr.slice(3, 11)); // x }); Deno.test("a region outside the moving image is empty", async () => { @@ -321,12 +325,36 @@ Deno.test("RAS orientation yields a non-identity direction", async () => { { z: 4, y: 8, x: 16 }, { z: 1, y: 1, x: 1 }, { z: 0, y: 0, x: 0 }, - RAS as never, + RAS, ); const direction = itkDirection(image, ["x", "y", "z"]); assertEquals(Array.from(direction), [-1, 0, 0, 0, -1, 0, 0, 0, 1]); }); +Deno.test("a 3D-only orientation on a 2D image falls back to identity", async () => { + const { itkDirection } = await import( + "../src/io/itk_transform_resample_bounding_box-shared.ts" + ); + const { AnatomicalOrientationValues, createAnatomicalOrientation } = + await import("../src/types/rfc4.ts"); + // An inferior/superior axis points along LPS z; truncating its column to + // 2D would produce a singular matrix. + const image = await geometryImage( + ["y", "x"], + { y: 8, x: 16 }, + { y: 1, x: 1 }, + { y: 0, x: 0 }, + { + x: createAnatomicalOrientation(AnatomicalOrientationValues.LeftToRight), + y: createAnatomicalOrientation( + AnatomicalOrientationValues.InferiorToSuperior, + ), + }, + ); + const direction = itkDirection(image, ["x", "y"]); + assertEquals(Array.from(direction), [1, 0, 0, 1]); +}); + Deno.test("an index range overflow is reported, not silently empty", async () => { // The region is computed in 32-bit index space while the corners come back // as doubles, so a fixed/moving scale mismatch this large wraps the integer From 1c69e4114596c9801b8b93e20e35fee3826f1991 Mon Sep 17 00:00:00 2001 From: Valentin Boussot Date: Tue, 11 Aug 2026 14:02:25 +0200 Subject: [PATCH 7/7] fix(py): correct the declared value type of float32 ITK transforms itk.dict_from_transform materializes the parameters as float64 while keeping the transform's declared float32 value type, so the pipeline read the buffer as raw float32 and computed a region from garbage for e.g. DisplacementFieldTransform[itk.F, 2]. Declare the type the buffer actually has. --- .../itk_transform_resample_bounding_box.py | 11 ++++++ ...est_itk_transform_resample_bounding_box.py | 34 ++++++++++++++++--- 2 files changed, 41 insertions(+), 4 deletions(-) diff --git a/py/ngff_zarr/itk_transform_resample_bounding_box.py b/py/ngff_zarr/itk_transform_resample_bounding_box.py index f7ac9066..89a921c6 100644 --- a/py/ngff_zarr/itk_transform_resample_bounding_box.py +++ b/py/ngff_zarr/itk_transform_resample_bounding_box.py @@ -311,6 +311,17 @@ def _transform_from_dict(entry: dict): entry = dict(entry) transform_type = entry.get("transformType") if isinstance(transform_type, dict): + # itk.dict_from_transform always materializes the parameters as + # float64 but keeps the transform's declared value type, so a + # float-parameterized transform (e.g. DisplacementFieldTransform + # [itk.F, 2]) arrives declared float32 with a float64 buffer. The + # pipeline then reads the buffer as float32 and computes a region + # from garbage. Declare the type the buffer actually has. + parameters = entry.get("parameters") + if parameters is not None: + actual = str(np.asarray(parameters).dtype) + if actual in ("float32", "float64"): + transform_type = {**transform_type, "parametersValueType": actual} entry["transformType"] = TransformType(**transform_type) return ItkTransform(**entry) diff --git a/py/test/test_itk_transform_resample_bounding_box.py b/py/test/test_itk_transform_resample_bounding_box.py index 754ca96b..12bc12b9 100644 --- a/py/test/test_itk_transform_resample_bounding_box.py +++ b/py/test/test_itk_transform_resample_bounding_box.py @@ -394,9 +394,11 @@ def test_degenerate_fixed_grid_yields_an_empty_region(): assert bounding_box.crop(moving) is None -def _constant_displacement_field(itk, shift, size=8, spacing=8.0): +def _constant_displacement_field(itk, shift, size=8, spacing=8.0, ctype=None): """A displacement field that shifts every point by ``shift`` (ITK order).""" - field = itk.Image[itk.Vector[itk.D, 2], 2].New() + if ctype is None: + ctype = itk.D + field = itk.Image[itk.Vector[ctype, 2], 2].New() region = itk.ImageRegion[2]() extent = itk.Size[2]() extent[0], extent[1] = size, size @@ -405,11 +407,11 @@ def _constant_displacement_field(itk, shift, size=8, spacing=8.0): field.SetSpacing([spacing, spacing]) field.SetOrigin([0.0, 0.0]) field.Allocate() - offset = itk.Vector[itk.D, 2]() + offset = itk.Vector[ctype, 2]() offset[0], offset[1] = shift field.FillBuffer(offset) - transform = itk.DisplacementFieldTransform[itk.D, 2].New() + transform = itk.DisplacementFieldTransform[ctype, 2].New() transform.SetDisplacementField(field) return transform @@ -440,6 +442,30 @@ def test_non_linear_displacement_field_is_supported(): assert bounding_box.size == {"y": 34, "x": 34} +def test_float_displacement_field_matches_double(): + """A float32-parameterized transform yields the double-precision region. + + itk.dict_from_transform materializes the parameters as float64 while + declaring the transform's own value type, so a float32 declaration over + a float64 buffer must be corrected, not read as raw float32. + """ + itk = pytest.importorskip("itk") + + fixed = _image("yx", {"y": 32, "x": 32}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + moving = _image("yx", {"y": 256, "x": 256}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) + + reference = itk_transform_resample_bounding_box( + _constant_displacement_field(itk, (5.0, 3.0)), fixed, moving, padding=1 + ) + float_transform = _constant_displacement_field(itk, (5.0, 3.0), ctype=itk.F) + bounding_box = itk_transform_resample_bounding_box( + float_transform, fixed, moving, padding=1 + ) + + assert bounding_box.start_index == reference.start_index == {"y": 2, "x": 4} + assert bounding_box.size == reference.size == {"y": 34, "x": 34} + + def test_unsupported_transform_type_is_rejected(): fixed = _image("yx", {"y": 4, "x": 4}, {"y": 1, "x": 1}, {"y": 0, "x": 0}) with pytest.raises(TypeError, match="unsupported transform type"):