Add itk_transform_resample_bounding_box

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));
+}
+```

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",