Nipy / scverse / ITK spatial data convergence

[effigies]

At the numfocus summit, I've been talking with scverse and ITK folks, and there's definitely a sense that we're solving similar problems and probably should not be entirely duplicating effort. There's some work that already exists that overlaps with what we're trying to do here and is also intended to work with Dask, xarray and other structures beyond numpy ndarrays.

I and others worked on a doc where we've scoped out who has done what: https://demo.hedgedoc.org/VlVvlbOIS2qbDHIzXteKIw#

Current text follows:

---

Nipy / scverse / ITK spatial data convergence

This document is a product of the 2023 Numfocus summit. Its purpose is to identify points of convergence in spatial data and transformations.

Date: 13 September 2023
Authors:

• Chris Markiewicz (@effigies)
• Matt McCormick (@thewtex)

Existing structures, gaps and efforts to fill the gaps

Nipy (Neuroimaging in Python) structures

• nibabel.spatialimages.SpatialImage - Despite the implementations, this is mostly an abstract class. The concrete class people use would mostly be Nifti1Image, most of whose implementation is spread through the Nifti1Pair and AnalyzeImage superclasses.
  • SpatialImages are a (data, affine) pair, along with format-specific headers/metadata. The intrinsic affine is essential to all neuroimaging formats and use-cases.
• Meshes are currently heterogeneously defined, tightly tied to file formats. We have BIAP9 - The Coordinate Image API where we propose a generic structure for meshes/pointclouds and data sampled to them.
  • The goal is essentially to replicate the success of SpatialImage, tying sampled data to the manifold they are sampled on, while not requiring any more data be loaded into memory than is needed for the current operation.
  • Class hierarchy
    • Pointset - coordinates
      • TriangularMesh - Coordinates + adjacency data
      • ImageGrid - Coordinates of a SpatialImage (image affine applied to all data array indices)
• nitransforms is a project that grew out of a desire to unify the different transformation file formats and apply them in Python. It has turned out to have convergence with coordinate images, as transformations map points between spaces.
  • Class Hierarchy
    • SampledSpatialData is a representation of a point set as rows of coordinates. (Typically Nx3)
      • ImageGrid provides a SampledSpatialData view of a nibabel.spatialimages.SpatialImage
    • TransformBase - All transforms are mappings from coordinates in the target space to coordinates in the source space. Source and target spaces are RAS+ oriented, meaning an increasing x, y, or z value indicates movement to the right, anterior or superior of the origin.
      • Affine - 12 parameter relation of one space to another.
        • LinearTransformsMapping - series of affines that apply to images in a series. Canonically this would encode motion-correction affines to resample volumes to a reference.
      • DenseFieldTransform - ND array of displacements or deformations. Displacements are difference vectors to target locations, while deformations are target locations.
      • BSplineFieldTransform is a field transform represented as coefficients to a B-spline basis, the knot-spacing of which is encoded in the affine of the image. This produces a more compact field representation and precisely defines the interpolation.
      • TransformChain - List of transforms to apply from left to right

scverse structures

• spatialdata package
  • Transformations: spatialdata.transformations.transformations
    • BaseTransformation - Transform implementations, abstracted from serialization
    • NgffBaseTransformation - Serializations of transforms into NGFF (.ome.zarr)
      • Has input and output coordinate systems, which are (name, type[, unit]) structures stored in the metadata dictionary, e.g.
        {"type": "identity", "input": {"name": "CoordSysA", "axes": [{"type": "space", "name": "X"}, {"type": "space", "name": "Y"}]}, "input": {"name": "CoordSysB", "axes": [{"type": "space", "name": "X"}, {"type": "space", "name": "Y"}]}}
      • Subclasses:
        • Identity
        • Axis mapping: dict
        • Translation: JSON array of one float per dimension
        • Scale: JSON array of one float per dimension
        • Affine: JSON array of 6/12 parameters (for 2/3D cases)
        • Rotation: JSON array of unraveled rotation matrix
        • Sequence: list of transformations
        • NgffByDimension: ???
        • Unimplemented: displacements, inverse, bijection
      • Spec: https://github.com/bogovicj/ngff/blob/coord-transforms/latest/index.bs
        • Displacements at least are a zarr dataset
  • Data models
    • 2D/3D images and labels
    • Points (point clouds?)
    • Shapes (parametrically defined region?)
      • Polygons - sequence of world coordinates
    • Table (annotations on regions, which are labels or shapes)
      • AnnData
      • Never transformed/interpolated

ITK structures

• ITK supports a number of mature, well-defined transform types

  • They can be applied to N-dimensional images, label images, meshes, spatial objects aka geometries aka parametric shapes, and point sets aka point clouds.
  • They include methods to transform points, vectors, co-variant vectors, and second-rank tensors
  • They include jacobian implementations for registration and other use cases

• Definitions of the transforms are well-defined, including a human-readale description of the transform and its use cases, mathematical equations that define a transform, and a description of the transform parameters

  • This can be found in Book 2, Section 3.9 of the ITK Software Guide
  • Short Jupyter notebook tutorial on how transforms work with spatial images

• Example transforms

  • DisplacementFieldTransform - Displacements (deltas)
    • For anatomical images, these are displacements in LPS world space (like the image definition).
  • AffineTransform - dim*(dim-1)-parameter transform
  • BSplineTransform
  • CompositeTransform - a sequence of transforms
  • VersorTransform - rotations parameterized for manipulation and optimization
  • Euler2DTransform, Euler3DTransform, ScaleTransform, TranslationTransform, ...

• Neuroimaging details:

  • ITK's has support for many file formats used in neuroimaging and beyond
  • ITK loads images into Left-Posterior-Superior (LPS) orientation because this orientation is used by the clinical industry standard format for medical images, DICOM
  • NIfTI images are read into an LPS+ world space like other formats, so NIfTIs are reoriented on read/write. ITK needs to transform RAS-to-RAS transforms to LPS-to-LPS on ingestion, and neuroimaging tools that support ITK transforms convert to RAS-to-RAS.
    • For affines, this is as simple (in both directions) as. $\mathrm{diag}(-1, -1, 1, 1)\times A\times\mathrm{diag}(-1, -1, 1, 1)$.
    • For displacement fields, multiplying the first two volumes by $-1$ works.
    • For more detailed information and exercises and see this LPS-RAS-4x4 Jupyter Notebook tutorial

• File formats

  • ITK supports many transform files formats including:
    • An HDF5-based format (recommended)
    • A text-based format
    • MINC (Montreal Neuroimaging HDF5-based) format
    • .mat Matlab file format
    • DICOM
  • OME-Zarr support is in progress

spatial-image

• https://github.com/spatial-image/spatial-image provides a collection of classes that encode 1-5 dimensional arrays, where the dimensions are (channel, x, y, z, t). The goal is seamless interoperability in many contexts.
• Structures compatible with OME-Zarr

X5 (HDF5-based transform file format)

• Grew out of BIDS Extension Proposal 14
• Spec?