Convert coordinate transformations between RFC-5 and ITK

docs/rfc5.md

@@ -222,6 +222,39 @@ Use `Coordinates` instead of `Displacements` (and `axes_types={"c":
 "coordinate"}` on the field) for an absolute coordinate field; the store layout
 is identical.
 
+### Interoperating with ITK
+
+Every transformation in the table above converts to ITK with
+`ngff_transform_to_itk_transform`, and `itk_transform_to_ngff_transform`
+converts back. The second is how a registration result gets into the store:
+convert the `CompositeTransform` an Elastix registration returns and attach it
+to the multiscales metadata as shown above.
+
+`identity`, `scale`, `translation`, `rotation`, `affine`, `mapAxis`,
+`byDimension`, `bijection` and any `sequence` of them describe a linear mapping
+and are folded into the single affine ITK gets. A `mapAxis` becomes its
+permutation matrix, a `byDimension` writes each item into the rows its
+`outputAxes` name, and a `bijection` contributes its `forward` direction.
+
+Both directions reconcile the places where the conventions differ: RFC-5 orders
+parameters in Zarr axis order while ITK orders them fastest-axis-first, an
+RFC-5 `sequence` applies its first entry first while an ITK transform list
+applies its last entry first, and ITK's center of rotation is folded into the
+offset since an RFC-5 affine has none.
+
+A `displacements` or `coordinates` transformation converts too, with
+`itk_displacement_field_to_ngff_transform` and its inverse: the field is an
+array rather than a handful of numbers, so those return, and take, the field
+image alongside the transformation. Both reach ITK as a displacement field,
+since ITK has no absolute-coordinate transform; coming back, a field is
+`displacements`.
+
+Building on that, `resample_bounding_box` computes which region of a moving
+image a resample through the transformation would read, from geometry alone,
+and `resample` resamples the grid block by block through the same
+transformation. See [Out-of-core resampling](./itk.md#out-of-core-resampling)
+and [Converting transforms](./itk.md#converting-transforms).
+
 ## TypeScript
 
 The TypeScript package (`@fideus-labs/ngff-zarr`) mirrors the Python API. Field

py/examples/itk_elastix_transform_resample_s3.ipynb

@@ -16,8 +16,8 @@
     "   registers the coarsest level of the pyramids, which is small enough to fetch\n",
     "   in a few seconds. The transform lives in physical space, so it applies at\n",
     "   any level.\n",
-    "3. `itk_transform_resample` applies the transform at a fine level. For every\n",
-    "   output block it asks `itk_transform_resample_bounding_box` which region of\n",
+    "3. `resample` applies the transform at a fine level. For every\n",
+    "   output block it asks `resample_bounding_box` which region of\n",
     "   the moving image is needed, reads only the S3 chunks inside it, and\n",
     "   resamples through `itkwasm-downsample`. The blocks are tasks of one Dask\n",
     "   graph that reference the moving chunks directly, so a chunk that several\n",
@@ -65,9 +65,9 @@
     "from ngff_zarr import (\n",
     "    NgffImage,\n",
     "    from_ome_zarr,\n",
-    "    itk_transform_resample,\n",
-    "    itk_transform_resample_bounding_box,\n",
     "    ngff_image_to_itk_image,\n",
+    "    resample,\n",
+    "    resample_bounding_box,\n",
     "    to_multiscales,\n",
     "    to_ngff_image,\n",
     "    to_ome_zarr,\n",
@@ -285,7 +285,7 @@
    "source": [
     "## Resample a fine level, streaming from S3 to a local OME-Zarr\n",
     "\n",
-    "`itk_transform_resample` returns a lazy `NgffImage` on the fixed grid; nothing\n",
+    "`resample` returns a lazy `NgffImage` on the fixed grid; nothing\n",
     "has been read yet. `to_ome_zarr` then drives the computation: each output block\n",
     "reads the S3 chunks inside its own moving-image bounding box, resamples them,\n",
     "and is written to disk as soon as it is done, with sixteen blocks in flight at\n",
@@ -326,9 +326,7 @@
    "source": [
     "FINE = 2\n",
     "\n",
-    "resampled_brain = itk_transform_resample(\n",
-    "    brain_transform, fixed_levels[FINE], moving_levels[FINE]\n",
-    ")\n",
+    "resampled_brain = resample(brain_transform, fixed_levels[FINE], moving_levels[FINE])\n",
     "print(\n",
     "    f\"Output grid: shape={resampled_brain.data.shape}, blocks={resampled_brain.data.npartitions}\"\n",
     ")\n",
@@ -393,7 +391,7 @@
     "## The regions behind the blocks\n",
     "\n",
     "Each output block reads its own region of the moving image, and\n",
-    "`itk_transform_resample_bounding_box` computes that region for any sub-grid of\n",
+    "`resample_bounding_box` computes that region for any sub-grid of\n",
     "the fixed image from geometry alone. The views below use level 2 with 64³\n",
     "output blocks and one axial slice through the middle of the fixed brain; every\n",
     "block whose z-range contains that slice is drawn on it, and the moving-image\n",
@@ -453,7 +451,7 @@
     "        scale=fixed_view.scale,\n",
     "        translation=translation,\n",
     "    )\n",
-    "    return itk_transform_resample_bounding_box(\n",
+    "    return resample_bounding_box(\n",
     "        brain_transform, block_grid, moving_view, padding=padding\n",
     "    )\n",
     "\n",
@@ -623,7 +621,7 @@
     "## Benchmark\n",
     "\n",
     "Two knobs decide what a resample costs, and neither is an argument of\n",
-    "`itk_transform_resample`: the size of an output block is the chunking of the\n",
+    "`resample`: the size of an output block is the chunking of the\n",
     "fixed image, and the number of blocks in flight is Dask's `num_workers`.\n",
     "`benchmark` takes both, and reports the graph build, the wall clock of the\n",
     "write, the throughput, and the peak resident memory of the process while it\n",
@@ -702,7 +700,7 @@
     "        fixed = replace(fixed, data=fixed.data.rechunk((block,) * fixed.data.ndim))\n",
     "\n",
     "    start = time.perf_counter()\n",
-    "    image = itk_transform_resample(brain_transform, fixed, moving_levels[level])\n",
+    "    image = resample(brain_transform, fixed, moving_levels[level])\n",
     "    build = time.perf_counter() - start\n",
     "\n",
     "    store = f\"benchmark_level{level}.ome.zarr\"\n",

py/ngff_zarr/__init__.py

@@ -14,6 +14,10 @@
 )
 from .config import config
 from .detect_cli_io_backend import ConversionBackend, detect_cli_io_backend
+from .displacement_field_transform import (
+    itk_displacement_field_to_ngff_transform,
+    ngff_displacement_field_to_itk_transform,
+)
 from .from_ngff_zarr import from_ngff_zarr, from_ome_zarr
 from .hcs import (
     HCSPlate,
@@ -24,10 +28,9 @@
     write_hcs_well_image,
 )
 from .itk_image_to_ngff_image import itk_image_to_ngff_image
-from .itk_transform_resample import itk_transform_resample
-from .itk_transform_resample_bounding_box import (
-    ResampleBoundingBox,
-    itk_transform_resample_bounding_box,
+from .itk_transform_to_ngff_transform import (
+    itk_transform_to_ngff_matrix,
+    itk_transform_to_ngff_transform,
 )
 from .lif_to_ngff_image import (
     has_mosaic_dimension,
@@ -40,10 +43,16 @@
 from .multiscales import Multiscales, NgffMultiscales
 from .ngff_image import NgffImage
 from .ngff_image_to_itk_image import ngff_image_to_itk_image
+from .ngff_transform_to_itk_transform import ngff_transform_to_itk_transform
 from .nibabel_image_to_ngff_image import (
     extract_omero_metadata_from_nibabel,
     nibabel_image_to_ngff_image,
 )
+from .resample import resample
+from .resample_bounding_box import (
+    ResampleBoundingBox,
+    resample_bounding_box,
+)
 from .rfc4 import (
     LPS,
     RAS,
@@ -124,9 +133,15 @@
     "nibabel_image_to_ngff_image",
     "extract_omero_metadata_from_nibabel",
     "ngff_image_to_itk_image",
+    # RFC 5 - Coordinate transformations and ITK
+    "ngff_transform_to_itk_transform",
+    "itk_transform_to_ngff_matrix",
+    "itk_transform_to_ngff_transform",
+    "itk_displacement_field_to_ngff_transform",
+    "ngff_displacement_field_to_itk_transform",
     # Out-of-core resampling
-    "itk_transform_resample",
-    "itk_transform_resample_bounding_box",
+    "resample",
+    "resample_bounding_box",
     "ResampleBoundingBox",
     "memory_usage",
     "task_count",