wip

Author: oesteban

nitransforms/io/itk.py

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

nitransforms/nonlinear.py

@@ -188,6 +188,7 @@ def map(self, x, inverse=False):
         ijk = self.reference.index(x)
         indexes = np.round(ijk).astype("int")
 
+        import pdb; pdb.set_trace()
         if np.all(np.abs(ijk - indexes) < 1e-3):
             indexes = tuple(tuple(i) for i in indexes)
             return self._field[indexes]

nitransforms/tests/test_io.py

@@ -710,7 +710,8 @@ def test_itk_disp_load_intent():
 
 # Added tests for displacements fields orientations (ANTs/ITK)
 @pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", "oblique"])
-def test_itk_displacements(tmp_path, get_testdata, image_orientation):
+@pytest.mark.parametrize("field_is_random", [False, True])
+def test_itk_displacements(tmp_path, get_testdata, image_orientation, field_is_random):
     """Exercise I/O of ITK displacements fields."""
 
     nii = get_testdata[image_orientation]
@@ -719,13 +720,17 @@ def test_itk_displacements(tmp_path, get_testdata, image_orientation):
     shape = nii.shape
     ref_affine = nii.affine.copy()
 
-    field = np.hstack(
-        (
-            np.linspace(-50, 50, num=np.prod(shape)),
-            np.linspace(-80, 80, num=np.prod(shape)),
-            np.zeros(np.prod(shape)),
-        )
-    ).reshape(shape + (3,))
+    field = (
+        np.hstack(
+            (
+                np.linspace(-50, 50, num=np.prod(shape)),
+                np.linspace(-80, 80, num=np.prod(shape)),
+                np.zeros(np.prod(shape)),
+            )
+        ).reshape(shape + (3,))
+        if not field_is_random
+        else np.random.normal(size=shape + (3,))
+    )
 
     nit_nii = itk.ITKDisplacementsField.to_image(
         nb.Nifti1Image(field, ref_affine, None)

nitransforms/tests/test_nonlinear.py

@@ -3,10 +3,15 @@
 """Tests of nonlinear transforms."""
 
 import os
+from subprocess import check_call
+import shutil
+
+import SimpleITK as sitk
 import pytest
 
 import numpy as np
 import nibabel as nb
+from nibabel.affines import from_matvec
 from nitransforms.resampling import apply
 from nitransforms.base import TransformError
 from nitransforms.io.base import TransformFileError
@@ -15,7 +20,7 @@
     DenseFieldTransform,
 )
 from nitransforms import io
-from ..io.itk import ITKDisplacementsField
+from nitransforms.io.itk import ITKDisplacementsField
 
 
 @pytest.mark.parametrize("size", [(20, 20, 20), (20, 20, 20, 3)])
@@ -34,16 +39,6 @@ def test_displacements_bad_sizes(size):
         DenseFieldTransform(nb.Nifti1Image(np.zeros(size), np.eye(4), None))
 
 
-def test_itk_disp_load_intent():
-    """Checks whether the NIfTI intent is fixed."""
-    with pytest.warns(UserWarning):
-        field = ITKDisplacementsField.from_image(
-            nb.Nifti1Image(np.zeros((20, 20, 20, 1, 3)), np.eye(4), None)
-        )
-
-    assert field.header.get_intent()[0] == "vector"
-
-
 def test_displacements_init():
     identity1 = DenseFieldTransform(
         np.zeros((10, 10, 10, 3)),
@@ -67,6 +62,30 @@ def test_displacements_init():
         )
 
 
+@pytest.mark.parametrize("is_deltas", [True, False])
+def test_densefield_oob_resampling(is_deltas):
+    """Ensure mapping outside the field returns input coordinates."""
+    ref = nb.Nifti1Image(np.zeros((2, 2, 2), dtype="uint8"), np.eye(4))
+
+    if is_deltas:
+        field = nb.Nifti1Image(np.ones((2, 2, 2, 3), dtype="float32"), np.eye(4))
+    else:
+        grid = np.stack(
+            np.meshgrid(*[np.arange(2) for _ in range(3)], indexing="ij"),
+            axis=-1,
+        ).astype("float32")
+        field = nb.Nifti1Image(grid + 1.0, np.eye(4))
+
+    xfm = DenseFieldTransform(field, is_deltas=is_deltas, reference=ref)
+
+    points = np.array([[-1.0, -1.0, -1.0], [0.5, 0.5, 0.5], [3.0, 3.0, 3.0]])
+    mapped = xfm.map(points)
+
+    assert np.allclose(mapped[0], points[0])
+    assert np.allclose(mapped[2], points[2])
+    assert np.allclose(mapped[1], points[1] + 1)
+
+
 def test_bsplines_init():
     with pytest.raises(TransformError):
         BSplineFieldTransform(
@@ -122,76 +141,6 @@ def test_bspline(tmp_path, testdata_path):
     )
 
 
-@pytest.mark.parametrize("is_deltas", [True, False])
-def test_densefield_x5_roundtrip(tmp_path, is_deltas):
-    """Ensure dense field transforms roundtrip via X5."""
-    ref = nb.Nifti1Image(np.zeros((2, 2, 2), dtype="uint8"), np.eye(4))
-    disp = nb.Nifti1Image(np.random.rand(2, 2, 2, 3).astype("float32"), np.eye(4))
-
-    xfm = DenseFieldTransform(disp, is_deltas=is_deltas, reference=ref)
-
-    node = xfm.to_x5(metadata={"GeneratedBy": "pytest"})
-    assert node.type == "nonlinear"
-    assert node.subtype == "densefield"
-    assert node.representation == "displacements" if is_deltas else "deformations"
-    assert node.domain.size == ref.shape
-    assert node.metadata["GeneratedBy"] == "pytest"
-
-    fname = tmp_path / "test.x5"
-    io.x5.to_filename(fname, [node])
-
-    xfm2 = DenseFieldTransform.from_filename(fname, fmt="X5")
-
-    assert xfm2.reference.shape == ref.shape
-    assert np.allclose(xfm2.reference.affine, ref.affine)
-    assert xfm == xfm2
-
-
-def test_bspline_to_x5(tmp_path):
-    """Check BSpline transforms export to X5."""
-    coeff = nb.Nifti1Image(np.zeros((2, 2, 2, 3), dtype="float32"), np.eye(4))
-    ref = nb.Nifti1Image(np.zeros((2, 2, 2), dtype="uint8"), np.eye(4))
-
-    xfm = BSplineFieldTransform(coeff, reference=ref)
-    node = xfm.to_x5(metadata={"tool": "pytest"})
-    assert node.type == "nonlinear"
-    assert node.subtype == "bspline"
-    assert node.representation == "coefficients"
-    assert node.metadata["tool"] == "pytest"
-
-    fname = tmp_path / "bspline.x5"
-    io.x5.to_filename(fname, [node])
-
-    xfm2 = BSplineFieldTransform.from_filename(fname, fmt="X5")
-    assert np.allclose(xfm._coeffs, xfm2._coeffs)
-    assert xfm2.reference.shape == ref.shape
-    assert np.allclose(xfm2.reference.affine, ref.affine)
-
-
-@pytest.mark.parametrize("is_deltas", [True, False])
-def test_densefield_oob_resampling(is_deltas):
-    """Ensure mapping outside the field returns input coordinates."""
-    ref = nb.Nifti1Image(np.zeros((2, 2, 2), dtype="uint8"), np.eye(4))
-
-    if is_deltas:
-        field = nb.Nifti1Image(np.ones((2, 2, 2, 3), dtype="float32"), np.eye(4))
-    else:
-        grid = np.stack(
-            np.meshgrid(*[np.arange(2) for _ in range(3)], indexing="ij"),
-            axis=-1,
-        ).astype("float32")
-        field = nb.Nifti1Image(grid + 1.0, np.eye(4))
-
-    xfm = DenseFieldTransform(field, is_deltas=is_deltas, reference=ref)
-
-    points = np.array([[-1.0, -1.0, -1.0], [0.5, 0.5, 0.5], [3.0, 3.0, 3.0]])
-    mapped = xfm.map(points)
-
-    assert np.allclose(mapped[0], points[0])
-    assert np.allclose(mapped[2], points[2])
-    assert np.allclose(mapped[1], points[1] + 1)
-
-
 def test_bspline_map_gridpoints():
     """BSpline mapping matches dense field on grid points."""
     ref = nb.Nifti1Image(np.zeros((5, 5, 5), dtype="uint8"), np.eye(4))
@@ -243,3 +192,128 @@ def manual_map(x):
     pts = np.array([[1.2, 1.5, 2.0], [3.3, 1.7, 2.4]])
     expected = np.vstack([manual_map(p) for p in pts])
     assert np.allclose(bspline.map(pts), expected, atol=1e-6)
+
+
+def test_densefield_map_against_ants(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")
+
+    points = np.array(
+        [
+            [0.0, 0.0, 0.0],
+            [1.0, 2.0, 3.0],
+            [10.0, -10.0, 5.0],
+            [-5.0, 7.0, -2.0],
+            [-12.0, 12.0, 0.0],
+        ]
+    )
+    csvin = tmp_path / "points.csv"
+    np.savetxt(csvin, points, delimiter=",", header="x,y,z", comments="")
+
+    csvout = tmp_path / "out.csv"
+    cmd = f"antsApplyTransformsToPoints -d 3 -i {csvin} -o {csvout} -t {warpfile}"
+    exe = cmd.split()[0]
+    if not shutil.which(exe):
+        pytest.skip(f"Command {exe} not found on host")
+    check_call(cmd, shell=True)
+
+    ants_res = np.genfromtxt(csvout, delimiter=",", names=True)
+    ants_pts = np.vstack([ants_res[n] for n in ("x", "y", "z")]).T
+
+    xfm = DenseFieldTransform(ITKDisplacementsField.from_filename(warpfile))
+    mapped = xfm.map(points)
+
+    assert np.allclose(mapped, ants_pts, atol=1e-6)
+
+
+@pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", "oblique"])
+@pytest.mark.parametrize("gridpoints", [True, False])
+def test_constant_field_vs_ants(tmp_path, get_testdata, image_orientation, gridpoints):
+    """Create a constant displacement field and compare mappings."""
+
+    nii = get_testdata[image_orientation]
+
+    # Create a reference centered at the origin with various axis orders/flips
+    shape = nii.shape
+    ref_affine = nii.affine.copy()
+
+    field = np.hstack((
+        np.zeros(np.prod(shape)),
+        np.linspace(-80, 80, num=np.prod(shape)),
+        np.linspace(-50, 50, num=np.prod(shape)),
+    )).reshape(shape + (3, ))
+    fieldnii = nb.Nifti1Image(field, ref_affine, None)
+
+    warpfile = tmp_path / "itk_transform.nii.gz"
+    ITKDisplacementsField.to_filename(fieldnii, warpfile)
+
+    # Ensure direct (xfm) and ITK roundtrip (itk_xfm) are equivalent
+    xfm = DenseFieldTransform(fieldnii)
+    itk_xfm = DenseFieldTransform(ITKDisplacementsField.from_filename(warpfile))
+
+    assert xfm == itk_xfm
+    np.testing.assert_allclose(xfm.reference.affine, itk_xfm.reference.affine)
+    np.testing.assert_allclose(ref_affine, itk_xfm.reference.affine)
+    np.testing.assert_allclose(xfm.reference.shape, itk_xfm.reference.shape)
+    np.testing.assert_allclose(xfm._field, itk_xfm._field)
+
+    points = (
+        xfm.reference.ndcoords.T if gridpoints
+        else np.array(
+            [
+                [0.0, 0.0, 0.0],
+                [1.0, 2.0, 3.0],
+                [10.0, -10.0, 5.0],
+                [-5.0, 7.0, -2.0],
+                [12.0, 0.0, -11.0],
+            ]
+        )
+    )
+
+    mapped = xfm.map(points)
+    nit_deltas = mapped - points
+
+    if gridpoints:
+        np.testing.assert_array_equal(field, nit_deltas.reshape(*shape, -1))
+
+    csvin = tmp_path / "points.csv"
+    np.savetxt(csvin, points, delimiter=",", header="x,y,z", comments="")
+
+    csvout = tmp_path / "out.csv"
+    cmd = f"antsApplyTransformsToPoints -d 3 -i {csvin} -o {csvout} -t {warpfile}"
+    exe = cmd.split()[0]
+    if not shutil.which(exe):
+        pytest.skip(f"Command {exe} not found on host")
+    check_call(cmd, shell=True)
+
+    ants_res = np.genfromtxt(csvout, delimiter=",", names=True)
+    ants_pts = np.vstack([ants_res[n] for n in ("x", "y", "z")]).T
+
+    # if gridpoints:
+    #     ants_field = ants_pts.reshape(shape + (3, ))
+    #     diff = xfm._field[..., 0] - ants_field[..., 0]
+    #     mask = np.argwhere(np.abs(diff) > 1e-2)[:, 0]
+    #     assert len(mask) == 0, f"A total of {len(mask)}/{ants_pts.shape[0]} contained errors:\n{diff[mask]}"
+
+    #     diff = xfm._field[..., 1] - ants_field[..., 1]
+    #     mask = np.argwhere(np.abs(diff) > 1e-2)[:, 0]
+    #     assert len(mask) == 0, f"A total of {len(mask)}/{ants_pts.shape[0]} contained errors:\n{diff[mask]}"
+
+    #     diff = xfm._field[..., 2] - ants_field[..., 2]
+    #     mask = np.argwhere(np.abs(diff) > 1e-2)[:, 0]
+    #     assert len(mask) == 0, f"A total of {len(mask)}/{ants_pts.shape[0]} contained errors:\n{diff[mask]}"
+
+    ants_deltas = ants_pts - points
+    np.testing.assert_array_equal(nit_deltas, ants_deltas)
+    np.testing.assert_array_equal(mapped, ants_pts)
+
+    diff = mapped - ants_pts
+    mask = np.argwhere(np.abs(diff) > 1e-2)[:, 0]
+
+    assert len(mask) == 0, f"A total of {len(mask)}/{ants_pts.shape[0]} contained errors:\n{diff[mask]}"

nitransforms/tests/test_resampling.py

@@ -188,6 +188,8 @@ def test_displacements_field1(
 
     xfm = nitnl.load(xfm_fname, fmt=sw_tool)
 
+    import pdb; pdb.set_trace()
+
     # Then apply the transform and cross-check with software
     cmd = APPLY_NONLINEAR_CMD[sw_tool](
         transform=os.path.abspath(xfm_fname),
@@ -243,7 +245,7 @@ def test_displacements_field1(
     assert np.sqrt((diff[5:-5, 5:-5, 5:-5] ** 2).mean()) < 1e-6
 
 
-@pytest.mark.parametrize("sw_tool", ["itk", "afni"])
+@pytest.mark.parametrize("sw_tool", ["afni"])
 def test_displacements_field2(tmp_path, testdata_path, sw_tool):
     """Check a translation-only field on one or more axes, different image orientations."""
     os.chdir(str(tmp_path))
@@ -275,6 +277,7 @@ def test_displacements_field2(tmp_path, testdata_path, sw_tool):
     nt_moved = apply(xfm, img_fname, order=0)
     nt_moved.to_filename("nt_resampled.nii.gz")
     sw_moved.set_data_dtype(nt_moved.get_data_dtype())
+
     diff = np.asanyarray(
         sw_moved.dataobj, dtype=sw_moved.get_data_dtype()
     ) - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype())