FIX: Ordering of ITK displacements fields when reading and writing

env.yml

@@ -24,6 +24,8 @@ dependencies:
   - nitime=0.10
   - scikit-image=0.22
   - scikit-learn=1.4
+  # SimpleITK, so build doesn't complain about building scikit from sources
+  - simpleitk=2.4
   # Utilities
   - graphviz=9.0
   - pandoc=3.1

nitransforms/base.py

@@ -202,30 +202,26 @@ def inverse(self):
     def ndindex(self):
         """List the indexes corresponding to the space grid."""
         if self._ndindex is None:
-            indexes = tuple([np.arange(s) for s in self._shape])
-            self._ndindex = np.array(np.meshgrid(*indexes, indexing="ij")).reshape(
-                self._ndim, self._npoints
-            )
+            indexes = np.mgrid[
+                0:self._shape[0], 0:self._shape[1], 0:self._shape[2]
+            ]
+            self._ndindex = indexes.reshape((indexes.shape[0], -1)).T
         return self._ndindex
 
     @property
     def ndcoords(self):
         """List the physical coordinates of this gridded space samples."""
         if self._coords is None:
-            self._coords = np.tensordot(
-                self._affine,
-                np.vstack((self.ndindex, np.ones((1, self._npoints)))),
-                axes=1,
-            )[:3, ...]
+            self._coords = self.ras(self.ndindex)
         return self._coords
 
     def ras(self, ijk):
         """Get RAS+ coordinates from input indexes."""
-        return _apply_affine(ijk, self._affine, self._ndim)
+        return _apply_affine(ijk, self._affine, self._ndim).T
 
     def index(self, x):
         """Get the image array's indexes corresponding to coordinates."""
-        return _apply_affine(x, self._inverse, self._ndim)
+        return _apply_affine(x, self._inverse, self._ndim).T
 
     def _to_hdf5(self, group):
         group.attrs["Type"] = "image"

nitransforms/io/itk.py

@@ -1,4 +1,5 @@
 """Read/write ITK transforms."""
+
 import warnings
 import numpy as np
 from scipy.io import loadmat as _read_mat, savemat as _save_mat
@@ -138,8 +139,7 @@ def from_matlab_dict(cls, mdict, index=0):
         sa = tf.structarr
 
         affine = mdict.get(
-            "AffineTransform_double_3_3",
-            mdict.get("AffineTransform_float_3_3")
+            "AffineTransform_double_3_3", mdict.get("AffineTransform_float_3_3")
         )
 
         if affine is None:
@@ -337,7 +337,7 @@ def from_image(cls, imgobj):
         hdr = imgobj.header.copy()
         shape = hdr.get_data_shape()
 
-        if len(shape) != 5 or shape[-2] != 1 or not shape[-1] in (2, 3):
+        if len(shape) != 5 or shape[-2] != 1 or shape[-1] not in (2, 3):
             raise TransformFileError(
                 'Displacements field "%s" does not come from ITK.'
                 % imgobj.file_map["image"].filename
@@ -348,9 +348,11 @@ def from_image(cls, imgobj):
             hdr.set_intent("vector")
 
         field = np.squeeze(np.asanyarray(imgobj.dataobj))
-        field[..., (0, 1)] *= -1.0
-
-        return imgobj.__class__(field, imgobj.affine, hdr)
+        affine = imgobj.affine
+        midindex = (np.array(field.shape[:3]) - 1) * 0.5
+        offset = (LPS @ affine - affine) @ (*midindex, 1)
+        affine[:3, 3] += offset[:3]
+        return imgobj.__class__(np.flip(field, axis=(0, 1)), imgobj.affine, hdr)
 
     @classmethod
     def to_image(cls, imgobj):
@@ -359,10 +361,10 @@ def to_image(cls, imgobj):
         hdr = imgobj.header.copy()
         hdr.set_intent("vector")
 
-        warp_data = imgobj.get_fdata().reshape(imgobj.shape[:3] + (1, imgobj.shape[-1]))
-        warp_data[..., (0, 1)] *= -1
-
-        return imgobj.__class__(warp_data, imgobj.affine, hdr)
+        field = imgobj.get_fdata()
+        field = field.transpose(2, 1, 0, 3)[..., None, :]
+        field[..., (0, 1)] *= 1.0
+        return imgobj.__class__(field, LPS @ imgobj.affine, hdr)
 
 
 class ITKCompositeH5:
@@ -410,21 +412,16 @@ def from_h5obj(cls, fileobj, check=True, only_linear=False):
                 directions = np.reshape(_fixed[9:], (3, 3))
                 affine = from_matvec(directions * zooms, offset)
                 # ITK uses Fortran ordering, like NIfTI, but with the vector dimension first
-                field = np.moveaxis(
-                    np.reshape(
-                        xfm[f"{typo_fallback}Parameters"], (3, *shape.astype(int)), order='F'
-                    ),
-                    0,
-                    -1,
-                )
-                field[..., (0, 1)] *= -1.0
+                # In practice, this seems to work (see issue #171)
+                field = np.reshape(
+                    xfm[f"{typo_fallback}Parameters"], (*shape.astype(int), 3)
+                ).transpose(2, 1, 0, 3)
+
                 hdr = Nifti1Header()
                 hdr.set_intent("vector")
                 hdr.set_data_dtype("float")
 
-                xfm_list.append(
-                    Nifti1Image(field.astype("float"), LPS @ affine, hdr)
-                )
+                xfm_list.append(Nifti1Image(field.astype("float"), affine, hdr))
                 continue
 
             raise TransformIOError(

nitransforms/nonlinear.py

@@ -65,50 +65,47 @@ def __init__(self, field=None, is_deltas=True, reference=None):
         <DenseFieldTransform[3D] (57, 67, 56)>
 
         """
+
         if field is None and reference is None:
-            raise TransformError("DenseFieldTransforms require a spatial reference")
+            raise TransformError("cannot initialize field")
 
         super().__init__()
 
-        self._is_deltas = is_deltas
+        if field is not None:
+            field = _ensure_image(field)
+            # Extract data if nibabel object otherwise assume numpy array
+            _data = np.squeeze(
+                np.asanyarray(field.dataobj)
+                if hasattr(field, "dataobj")
+                else field.copy()
+            )
 
         try:
             self.reference = ImageGrid(reference if reference is not None else field)
         except AttributeError:
             raise TransformError(
-                "Field must be a spatial image if reference is not provided"
+                "field must be a spatial image if reference is not provided"
                 if reference is None
-                else "Reference is not a spatial image"
+                else "reference is not a spatial image"
             )
 
         fieldshape = (*self.reference.shape, self.reference.ndim)
-        if field is not None:
-            field = _ensure_image(field)
-            self._field = np.squeeze(
-                np.asanyarray(field.dataobj) if hasattr(field, "dataobj") else field
-            )
-            if fieldshape != self._field.shape:
-                raise TransformError(
-                    f"Shape of the field ({'x'.join(str(i) for i in self._field.shape)}) "
-                    f"doesn't match that of the reference({'x'.join(str(i) for i in fieldshape)})"
-                )
-        else:
-            self._field = np.zeros(fieldshape, dtype="float32")
-            self._is_deltas = True
-
-        if self._field.shape[-1] != self.ndim:
+        if field is None:
+            _data = np.zeros(fieldshape)
+        elif fieldshape != _data.shape:
             raise TransformError(
-                "The number of components of the field (%d) does not match "
-                "the number of dimensions (%d)" % (self._field.shape[-1], self.ndim)
+                f"Shape of the field ({'x'.join(str(i) for i in _data.shape)}) "
+                f"doesn't match that of the reference({'x'.join(str(i) for i in fieldshape)})"
             )
 
+        self._is_deltas = is_deltas
+        self._field = self.reference.ndcoords.reshape(fieldshape)
+
         if self.is_deltas:
-            self._deltas = (
-                self._field.copy()
-            )  # IMPORTANT: you don't want to update deltas
-            # Convert from displacements (deltas) to deformations fields
-            # (just add its origin to each delta vector)
-            self._field += self.reference.ndcoords.T.reshape(fieldshape)
+            self._deltas = _data.copy()
+            self._field += self._deltas
+        else:
+            self._field = _data.copy()
 
     def __repr__(self):
         """Beautify the python representation."""
@@ -153,7 +150,7 @@ def map(self, x, inverse=False):
         ...     test_dir / "someones_displacement_field.nii.gz",
         ...     is_deltas=False,
         ... )
-        >>> xfm.map([-6.5, -36., -19.5]).tolist()
+        >>> xfm.map([[-6.5, -36., -19.5]]).tolist()
         [[0.0, -0.47516798973083496, 0.0]]
 
         >>> xfm.map([[-6.5, -36., -19.5], [-1., -41.5, -11.25]]).tolist()
@@ -170,8 +167,8 @@ def map(self, x, inverse=False):
         ...     test_dir / "someones_displacement_field.nii.gz",
         ...     is_deltas=True,
         ... )
-        >>> xfm.map([[-6.5, -36., -19.5], [-1., -41.5, -11.25]]).tolist()
-        [[-6.5, -36.47516632080078, -19.5], [-1.0, -42.03835678100586, -11.25]]
+        >>> xfm.map([[-6.5, -36., -19.5], [-1., -41.5, -11.25]]).tolist()  # doctest: +ELLIPSIS
+        [[-6.5, -36.475..., -19.5], [-1.0, -42.038..., -11.25]]
 
         >>> np.array_str(
         ...     xfm.map([[-6.7, -36.3, -19.2], [-1., -41.5, -11.25]]),
@@ -185,18 +182,19 @@ def map(self, x, inverse=False):
         if inverse is True:
             raise NotImplementedError
 
-        ijk = self.reference.index(x)
+        ijk = self.reference.index(np.array(x, dtype="float32"))
         indexes = np.round(ijk).astype("int")
+        ongrid = np.where(np.linalg.norm(ijk - indexes, axis=1) < 1e-3)[0]
 
-        if np.all(np.abs(ijk - indexes) < 1e-3):
-            indexes = tuple(tuple(i) for i in indexes)
-            return self._field[indexes]
+        if ongrid.size == np.shape(x)[0]:
+            # return self._field[*indexes.T, :]  # From Python 3.11
+            return self._field[tuple(indexes.T) + (np.s_[:],)]
 
-        new_map = np.vstack(
+        mapped_coords = np.vstack(
             tuple(
                 map_coordinates(
                     self._field[..., i],
-                    ijk,
+                    ijk.T,
                     order=3,
                     mode="constant",
                     cval=np.nan,
@@ -207,8 +205,8 @@ def map(self, x, inverse=False):
         ).T
 
         # Set NaN values back to the original coordinates value = no displacement
-        new_map[np.isnan(new_map)] = np.array(x)[np.isnan(new_map)]
-        return new_map
+        mapped_coords[np.isnan(mapped_coords)] = np.array(x)[np.isnan(mapped_coords)]
+        return mapped_coords
 
     def __matmul__(self, b):
         """

nitransforms/resampling.py

@@ -253,7 +253,7 @@ def apply(
     serialize_4d = n_resamplings >= serialize_nvols
 
     targets = None
-    ref_ndcoords = _ref.ndcoords.T
+    ref_ndcoords = _ref.ndcoords
     if hasattr(transform, "to_field") and callable(transform.to_field):
         targets = ImageGrid(spatialimage).index(
             _as_homogeneous(
@@ -271,11 +271,8 @@ def apply(
             else targets
         )
 
-    if targets.ndim == 3:
-        targets = np.rollaxis(targets, targets.ndim - 1, 0)
-    else:
-        assert targets.ndim == 2
-        targets = targets[np.newaxis, ...]
+    if targets.ndim == 2:
+        targets = targets.T[np.newaxis, ...]
 
     if serialize_4d:
         data = (
@@ -290,6 +287,9 @@ def apply(
             (len(ref_ndcoords), n_resamplings), dtype=input_dtype, order="F"
         )
 
+        if targets.ndim == 3:
+            targets = np.rollaxis(targets, targets.ndim - 1, 1)
+
         resampled = asyncio.run(
             _apply_serial(
                 data,
@@ -311,6 +311,9 @@ def apply(
     else:
         data = np.asanyarray(spatialimage.dataobj, dtype=input_dtype)
 
+        if targets.ndim == 3:
+            targets = np.rollaxis(targets, targets.ndim - 1, 0)
+
         if data_nvols == 1 and xfm_nvols == 1:
             targets = np.squeeze(targets)
             assert targets.ndim == 2
@@ -320,15 +323,19 @@ def apply(
 
         if xfm_nvols > 1:
             assert targets.ndim == 3
-            n_time, n_dim, n_vox = targets.shape
+
+            # Targets must have shape (n_dim x n_time x n_vox)
+            n_dim, n_time, n_vox = targets.shape
             # Reshape to (3, n_time x n_vox)
-            ijk_targets = np.rollaxis(targets, 0, 2).reshape((n_dim, -1))
+            ijk_targets = targets.reshape((n_dim, -1))
             time_row = np.repeat(np.arange(n_time), n_vox)[None, :]
 
             # Now targets is (4, n_vox x n_time), with indexes (t, i, j, k)
             # t is the slowest-changing axis, so we put it first
             targets = np.vstack((time_row, ijk_targets))
             data = np.rollaxis(data, data.ndim - 1, 0)
+        else:
+            targets = targets.T
 
         resampled = ndi.map_coordinates(
             data,

nitransforms/tests/test_base.py

@@ -55,20 +55,24 @@ def test_ImageGrid(get_testdata, image_orientation):
 
     assert np.allclose(np.squeeze(img.ras(ijk[0])), xyz[0])
     assert np.allclose(np.round(img.index(xyz[0])), ijk[0])
-    assert np.allclose(img.ras(ijk).T, xyz)
-    assert np.allclose(np.round(img.index(xyz)).T, ijk)
+    assert np.allclose(img.ras(ijk), xyz)
+    assert np.allclose(np.round(img.index(xyz)), ijk)
 
     # nd index / coords
     idxs = img.ndindex
     coords = img.ndcoords
     assert len(idxs.shape) == len(coords.shape) == 2
-    assert idxs.shape[0] == coords.shape[0] == img.ndim == 3
-    assert idxs.shape[1] == coords.shape[1] == img.npoints == np.prod(im.shape)
+    assert idxs.shape[1] == coords.shape[1] == img.ndim == 3
+    assert idxs.shape[0] == coords.shape[0] == img.npoints == np.prod(im.shape)
 
     img2 = ImageGrid(img)
     assert img2 == img
     assert (img2 != img) is False
 
+    # Test indexing round trip
+    np.testing.assert_allclose(img.ndcoords, img.ras(img.ndindex))
+    np.testing.assert_allclose(img.ndindex, np.round(img.index(img.ndcoords)))
+
 
 def test_ImageGrid_utils(tmpdir, testdata_path, get_testdata):
     """Check that images can be objects or paths and equality."""