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_io.py

@@ -706,6 +706,41 @@ def test_itk_disp_load_intent():
 
 
 # Added tests for displacements fields orientations (ANTs/ITK)
+def test_itk_densefield_read(testdata_path, tmp_path):
+    """Map points with DenseFieldTransform and compare to ANTs."""
+    warpfile = (
+        testdata_path
+        / "regressions"
+        / ("01_ants_t1_to_mniComposite_DisplacementFieldTransform.nii.gz")
+    )
+    if not warpfile.exists():
+        pytest.skip("Composite transform test data not available")
+
+    nib_from_itk = nb.load(warpfile)
+    nit_from_itk = itk.ITKDisplacementsField.from_filename(warpfile)
+    itk_from_nit = itk.ITKDisplacementsField.to_image(nit_from_itk)
+
+    # Enable for debugging
+    # nib_from_itk.to_filename(tmp_path / "nib_from_itk.nii.gz")
+    # itk_from_nit.to_filename(tmp_path / "itk_from_nit.nii.gz")
+    # nit_from_itk.to_filename(tmp_path / "nit_from_itk.nii.gz")
+
+    # import pdb; pdb.set_trace()
+
+    np.testing.assert_allclose(
+        nib_from_itk.dataobj, itk_from_nit.dataobj, rtol=1e-5, atol=1e-5
+    )
+    np.testing.assert_allclose(
+        nit_from_itk.dataobj, np.squeeze(itk_from_nit.dataobj), rtol=1e-5, atol=1e-5
+    )
+    np.testing.assert_allclose(
+        nib_from_itk.affine, itk_from_nit.affine, rtol=1e-5, atol=1e-5
+    )
+    np.testing.assert_allclose(
+        nit_from_itk.affine, itk_from_nit.affine, rtol=1e-5, atol=1e-5
+    )
+
+
 @pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", "oblique"])
 @pytest.mark.parametrize("field_is_random", [False, True])
 def test_itk_displacements(tmp_path, get_testdata, image_orientation, field_is_random):
@@ -745,6 +780,15 @@ def test_itk_displacements(tmp_path, get_testdata, image_orientation, field_is_r
     itk_nit_nii = itk.ITKDisplacementsField.from_filename(itk_file)
     itk_nii = nb.load(itk_file)
 
+    # Test round trip
+    assert itk_nit_nii.shape == field.shape
+    np.testing.assert_allclose(itk_nit_nii.dataobj, field)
+    np.testing.assert_allclose(
+        itk_nit_nii.dataobj.transpose(2, 1, 0, 3)[..., None, :], nit_nii.dataobj
+    )
+    np.testing.assert_allclose(itk_nit_nii.affine, ref_affine)
+    np.testing.assert_allclose(LPS @ itk_nit_nii.affine, nit_nii.affine)
+
     assert nit_nii.shape == itk_nii.shape, (
         "ITK-generated and nitransforms-generated field shapes are different"
     )
@@ -754,14 +798,9 @@ def test_itk_displacements(tmp_path, get_testdata, image_orientation, field_is_r
     # Check ITK-generated field has LPS-rotated affine
     np.testing.assert_allclose(itk_nii.affine, LPS @ ref_affine)
     # Test ITK-generated dataobject vs. original field
-    np.testing.assert_allclose(itk_nii.dataobj, field.transpose(2, 1, 0, 3)[..., None, :])
-
-    # Test round trip
-    assert itk_nit_nii.shape == field.shape
-    np.testing.assert_allclose(itk_nit_nii.dataobj, field)
-    np.testing.assert_allclose(itk_nit_nii.dataobj.transpose(2, 1, 0, 3)[..., None, :], nit_nii.dataobj)
-    np.testing.assert_allclose(itk_nit_nii.affine, ref_affine)
-    np.testing.assert_allclose(LPS @ itk_nit_nii.affine, nit_nii.affine)
+    np.testing.assert_allclose(
+        itk_nii.dataobj, field.transpose(2, 1, 0, 3)[..., None, :]
+    )
 
 
 # Added tests for h5 orientation bug

nitransforms/tests/test_nonlinear.py

@@ -189,7 +189,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
@@ -218,13 +218,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"
         )
@@ -253,7 +254,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(
         (
@@ -262,6 +262,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"