fix: proper ordering of displacements fields when reading and writing

Transposes data as suggested by @feilong when reading and writing
displacements fields.

Tests have been added to ensure round-trip consistency.

Related-to: #171.

Co-Authored-By: Feilong Ma <[email protected]>

Author: oesteban

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:
@@ -196,7 +196,7 @@ def from_string(cls, string):
             lines = lines[1:]  # Drop banner with version
 
         parameters = np.eye(4, dtype="f4")
-        sa["index"] = int(lines[0][lines[0].index("T"):].split()[1])
+        sa["index"] = int(lines[0][lines[0].index("T") :].split()[1])
         sa["offset"] = np.genfromtxt(
             [lines[3].split(":")[-1].encode()], dtype=cls.dtype["offset"]
         )
@@ -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,10 @@ def from_image(cls, imgobj):
             warnings.warn("Incorrect intent identified.")
             hdr.set_intent("vector")
 
-        field = np.squeeze(np.asanyarray(imgobj.dataobj))
+        field = np.squeeze(np.asanyarray(imgobj.dataobj)).transpose(2, 1, 0, 3)
         field[..., (0, 1)] *= -1.0
 
-        return imgobj.__class__(field, imgobj.affine, hdr)
+        return imgobj.__class__(field, LPS @ imgobj.affine, hdr)
 
     @classmethod
     def to_image(cls, imgobj):
@@ -359,10 +359,9 @@ 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().transpose(2, 1, 0, 3)[..., None, :]
+        field[..., (0, 1)] *= -1.0
+        return imgobj.__class__(field, LPS @ imgobj.affine, hdr)
 
 
 class ITKCompositeH5:
@@ -410,21 +409,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(