enh: ITK's coordinate system deciphered

Author: oesteban

nitransforms/io/itk.py

@@ -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"]
         )
@@ -329,7 +329,15 @@ def from_h5obj(cls, fileobj, check=True):
 
 
 class ITKDisplacementsField(DisplacementsField):
-    """A data structure representing displacements fields."""
+    """
+    A data structure representing displacements fields.
+
+    Note that ITK's world coordinates are LPS:
+    http://sourceforge.net/p/advants/discussion/840261/thread/2a1e9307/.
+    This translates into the fact that `antsApplyTransformsToPoints` expects LPS coordinates,
+    and therefore, points must correct for the x and y directions before feeding into the tool.
+
+    """
 
     @classmethod
     def from_image(cls, imgobj):
@@ -347,24 +355,47 @@ def from_image(cls, imgobj):
             warnings.warn("Incorrect intent identified.")
             hdr.set_intent("vector")
 
-        field = np.squeeze(np.asanyarray(imgobj.dataobj))
-        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)
+        affine, qcode = hdr.get_qform(coded=True)
+
+        if qcode != 1:
+            warnings.warn(
+                "Displacements field has qform code={qcode}, which is not ITK-like. "
+                "Setting it to 1 (aligned with the image)."
+            )
+            affine = imgobj.affine
+            hdr.set_qform(imgobj.affine, code=1)
+            hdr.set_sform(imgobj.affine, code=0)
+
+        # ITK uses LPS coordinates, so first we patch the affine's offset
+        # This patch was developed under gh-266, by adapting
+        # nitransforms/tests/test_nonlinear.py::test_densefield_map_vs_ants[True]
+        # until the same delta maps were obtained with `antsApplyTransformsToPoints`
+        # and our NiTransforms transform's `.map()`.
+        mid_ijk = 0.5 * (np.array(imgobj.shape[:3]) - 1)
+        offset = (affine - LPS @ affine) @ (*mid_ijk, 1.0)
+        affine[:3, 3] -= offset[:3]
+
+        # Create a NIfTI image with the patched affine
+        data = np.squeeze(imgobj.get_fdata(dtype="float32"))
+        refmap = imgobj.__class__(np.flip(data, axis=(0, 1)), affine, hdr)
+
+        return refmap
 
     @classmethod
     def to_image(cls, imgobj):
         """Export a displacements field from a nibabel object."""
 
         hdr = imgobj.header.copy()
         hdr.set_intent("vector")
-
-        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)
+        hdr.set_data_dtype("float32")
+        affine = LPS @ imgobj.affine @ LPS
+        hdr.set_qform(affine, code=1)
+        hdr.set_sform(affine, code=0)
+        hdr.set_xyzt_units("mm", None)
+
+        return imgobj.__class__(
+            imgobj.get_fdata(dtype="float32")[..., None, :], affine, hdr
+        )
 
 
 class ITKCompositeH5:

nitransforms/tests/test_nonlinear.py

@@ -196,7 +196,7 @@ def test_densefield_map_vs_ants(testdata_path, tmp_path, ongrid):
     )
     if not warpfile.exists():
         pytest.skip("Composite transform test data not available")
-    
+
     nii = ITKDisplacementsField.from_filename(warpfile)
 
     # Get sampling indices
@@ -225,13 +225,14 @@ def test_densefield_map_vs_ants(testdata_path, tmp_path, ongrid):
         ants_mapped_xyz = ants_pts.reshape(*shape, 3)
         nit_mapped_xyz = mapped.reshape(*shape, 3)
 
+        nb.load(warpfile).to_filename(tmp_path / "original_ants_deltas.nii.gz")
+
         nb.Nifti1Image(coords_map, ref_affine, None).to_filename(
-            tmp_path / "baseline_field.nii.gz"
+            tmp_path / "baseline_positions.nii.gz"
         )
 
-        nb.Nifti1Image(ants_mapped_xyz, ref_affine, None).to_filename(
-            tmp_path / "ants_deformation_xyz.nii.gz"
-        )
+        nii.to_filename(tmp_path / "original_interpreted_deltas.nii.gz")
+
         nb.Nifti1Image(nit_mapped_xyz, ref_affine, None).to_filename(
             tmp_path / "nit_deformation_xyz.nii.gz"
         )
@@ -260,7 +261,6 @@ def test_constant_field_vs_ants(tmp_path, get_testdata, image_orientation, ongri
     )
 
     coords_map = grid_xyz.reshape(*shape, 3)
-    gold_mapped_xyz = coords_map + deltas
 
     deltas = np.hstack(
         (
@@ -269,6 +269,7 @@ def test_constant_field_vs_ants(tmp_path, get_testdata, image_orientation, ongri
             np.linspace(-50, 50, num=np.prod(shape)),
         )
     ).reshape(shape + (3,))
+    gold_mapped_xyz = coords_map + deltas
 
     fieldnii = nb.Nifti1Image(deltas, ref_affine, None)
     warpfile = tmp_path / "itk_transform.nii.gz"