FIX: Inconsistent flips of ITK's displacements fields - including loading from h5

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/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
@@ -347,10 +347,9 @@ def from_image(cls, imgobj):
             warnings.warn("Incorrect intent identified.")
             hdr.set_intent("vector")
 
-        field = np.squeeze(np.asanyarray(imgobj.dataobj))
-        field[..., (0, 1)] *= -1.0
+        field = np.squeeze(np.asanyarray(imgobj.dataobj)).transpose(2, 1, 0, 3)
 
-        return imgobj.__class__(field, imgobj.affine, hdr)
+        return imgobj.__class__(field, LPS @ imgobj.affine, hdr)
 
     @classmethod
     def to_image(cls, imgobj):
@@ -359,10 +358,8 @@ 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)
+        warp_data = imgobj.get_fdata().transpose(2, 1, 0, 3)[..., None, :]
+        return imgobj.__class__(warp_data, LPS @ imgobj.affine, hdr)
 
 
 class ITKCompositeH5:
@@ -410,21 +407,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

@@ -285,7 +285,7 @@ def to_x5(self, metadata=None):
         )
 
     @classmethod
-    def from_filename(cls, filename, fmt="X5", x5_position=0):
+    def from_filename(cls, filename, is_deltas=True, fmt="X5", x5_position=0):
         _factory = {
             "afni": io.afni.AFNIDisplacementsField,
             "itk": io.itk.ITKDisplacementsField,
@@ -299,7 +299,7 @@ def from_filename(cls, filename, fmt="X5", x5_position=0):
         if fmt == "X5":
             return from_x5(load_x5(filename), x5_position=x5_position)
 
-        return cls(_factory[fmt.lower()].from_filename(filename))
+        return cls(_factory[fmt.lower()].from_filename(filename), is_deltas=is_deltas)
 
 
 load = DenseFieldTransform.from_filename

nitransforms/resampling.py

@@ -254,22 +254,22 @@ def apply(
 
     targets = None
     ref_ndcoords = _ref.ndcoords.T
-    if hasattr(transform, "to_field") and callable(transform.to_field):
-        targets = ImageGrid(spatialimage).index(
-            _as_homogeneous(
-                transform.to_field(reference=reference).map(ref_ndcoords),
-                dim=_ref.ndim,
-            )
-        )
-    else:
-        # Targets' shape is (Nt, 3, Nv) with Nv = Num. voxels, Nt = Num. timepoints.
-        targets = (
-            ImageGrid(spatialimage).index(
-                _as_homogeneous(transform.map(ref_ndcoords), dim=_ref.ndim)
-            )
-            if targets is None
-            else targets
+    # if hasattr(transform, "to_field") and callable(transform.to_field):
+    #     targets = ImageGrid(spatialimage).index(
+    #         _as_homogeneous(
+    #             transform.to_field(reference=reference).map(ref_ndcoords),
+    #             dim=_ref.ndim,
+    #         )
+    #     )
+    # else:
+    # Targets' shape is (Nt, 3, Nv) with Nv = Num. voxels, Nt = Num. timepoints.
+    targets = (
+        ImageGrid(spatialimage).index(
+            _as_homogeneous(transform.map(ref_ndcoords), dim=_ref.ndim)
         )
+        if targets is None
+        else targets
+    )
 
     if targets.ndim == 3:
         targets = np.rollaxis(targets, targets.ndim - 1, 0)