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

@@ -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(

nitransforms/tests/test_resampling.py

@@ -175,19 +175,20 @@ def test_displacements_field1(
     msk.to_filename("mask.nii.gz")
 
     fieldmap = np.zeros(
-        (*nii.shape[:3], 1, 3) if sw_tool != "fsl" else (*nii.shape[:3], 3),
+        (*nii.shape[:3], 1, 3) if sw_tool == "afni" else (*nii.shape[:3], 3),
         dtype="float32",
     )
     fieldmap[..., axis] = -10.0
-
     _hdr = nii.header.copy()
-    if sw_tool in ("itk",):
-        _hdr.set_intent("vector")
     _hdr.set_data_dtype("float32")
 
-    xfm_fname = "warp.nii.gz"
     field = nb.Nifti1Image(fieldmap, nii.affine, _hdr)
-    field.to_filename(xfm_fname)
+
+    xfm_fname = "warp.nii.gz"
+    if sw_tool == "itk":
+        io.itk.ITKDisplacementsField.to_filename(field, xfm_fname)
+    else:
+        field.to_filename(xfm_fname)
 
     xfm = nitnl.load(xfm_fname, fmt=sw_tool)
 
@@ -197,7 +198,8 @@ def test_displacements_field1(
         reference=tmp_path / "mask.nii.gz",
         moving=tmp_path / "mask.nii.gz",
         output=tmp_path / "resampled_brainmask.nii.gz",
-        extra="--output-data-type uchar" if sw_tool == "itk" else "",
+        extra="",
+        # extra="--output-data-type uchar" if sw_tool == "itk" else "",
     )
 
     # skip test if command is not available on host
@@ -207,14 +209,18 @@ def test_displacements_field1(
 
     # resample mask
     exit_code = check_call([cmd], shell=True)
-    assert exit_code == 0
     sw_moved_mask = nb.load("resampled_brainmask.nii.gz")
     nt_moved_mask = apply(xfm, msk, order=0)
-    nt_moved_mask.set_data_dtype(msk.get_data_dtype())
-    diff = np.asanyarray(sw_moved_mask.dataobj) - np.asanyarray(nt_moved_mask.dataobj)
 
-    assert np.sqrt((diff**2).mean()) < RMSE_TOL_LINEAR
+    # Calculate xor between both:
+    sw_mask = np.asanyarray(sw_moved_mask.dataobj, dtype=bool)
     brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool)
+    percent_diff = (sw_mask != brainmask)[5:-5, 5:-5, 5:-5].sum() / brainmask.size
+
+    assert exit_code == 0
+    assert percent_diff < 1e-8, (
+        f"Resampled masks differed by {percent_diff * 100:0.2f}%."
+    )
 
     # Then apply the transform and cross-check with software
     cmd = APPLY_NONLINEAR_CMD[sw_tool](
@@ -236,8 +242,9 @@ def test_displacements_field1(
     diff = np.asanyarray(
         sw_moved.dataobj, dtype=sw_moved.get_data_dtype()
     ) - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype())
-    # A certain tolerance is necessary because of resampling at borders
-    assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL_LINEAR
+
+    # Drop samples close to the border of the image
+    assert np.sqrt((diff[5:-5, 5:-5, 5:-5] ** 2).mean()) < 1e-6
 
 
 @pytest.mark.xfail(

pyproject.toml

@@ -41,6 +41,8 @@ test = [
     "pytest-xdist >= 2.5",
     "coverage[toml] >= 5.2.1",
     "nitransforms[niftiext]",
+    "SimpleITK ~= 2.4",
+    "scikit-build",
 ]
 # Aliases
 niftiexts = ["nitransforms[niftiext]"]