Merge pull request #276 from nipy/fix/re-enable-tests

FIX: ``ImageGrid._coords`` was somehow overwritten + re-enable tests

Author: oesteban

nitransforms/base.py

@@ -206,14 +206,14 @@ def ndindex(self):
                 0:self._shape[0], 0:self._shape[1], 0:self._shape[2]
             ]
             self._ndindex = indexes.reshape((indexes.shape[0], -1)).T
-        return self._ndindex
+        return self._ndindex.copy()  # Return copies to disallow alteration
 
     @property
     def ndcoords(self):
         """List the physical coordinates of this gridded space samples."""
         if self._coords is None:
             self._coords = self.ras(self.ndindex)
-        return self._coords
+        return self._coords.copy()  # Return copies to disallow alteration
 
     def ras(self, ijk):
         """Get RAS+ coordinates from input indexes."""

nitransforms/nonlinear.py

@@ -257,6 +257,34 @@ def __eq__(self, other):
             warnings.warn("Fields are equal, but references do not match.")
         return _eq
 
+    def to_filename(self, filename, fmt="X5", moving=None, x5_inverse=False):
+        """Store the transform in the designated format."""
+
+        if fmt.upper() == "X5":
+            raise TypeError("Please use .to_x5()")
+
+        field = nb.Nifti1Image(
+            self._deltas if self.is_deltas else self._field,
+            self.reference.affine,
+            None,
+        )
+
+        if fmt.lower() == "afni":
+            from nitransforms.io.afni import AFNIDisplacementsField as FieldIOType
+
+        elif fmt.lower() in ("itk", "ants", "elastix"):
+            from nitransforms.io.itk import ITKDisplacementsField as FieldIOType
+
+        elif fmt.lower() == "fsl":
+            from nitransforms.io.fsl import FSLDisplacementsField as FieldIOType
+
+        else:
+            raise NotImplementedError(
+                f"Dense field of type '{fmt}' cannot be converted."
+            )
+
+        FieldIOType.to_image(field).to_filename(filename)
+
     def to_x5(self, metadata=None):
         """Return an :class:`~nitransforms.io.x5.X5Transform` representation."""
         metadata = {"WrittenBy": f"NiTransforms {__version__}"} | (metadata or {})

nitransforms/resampling.py

@@ -254,6 +254,7 @@ def apply(
 
     targets = None
     ref_ndcoords = _ref.ndcoords
+
     if hasattr(transform, "to_field") and callable(transform.to_field):
         targets = ImageGrid(spatialimage).index(
             _as_homogeneous(

nitransforms/tests/test_base.py

@@ -65,14 +65,15 @@ def test_ImageGrid(get_testdata, image_orientation):
     assert idxs.shape[1] == coords.shape[1] == img.ndim == 3
     assert idxs.shape[0] == coords.shape[0] == img.npoints == np.prod(im.shape)
 
+    # Test indexing round trip
+    np.testing.assert_allclose(coords, img.ras(idxs))
+    np.testing.assert_allclose(idxs, img.index(coords), rtol=1e-3, atol=1e-3)
+
+    # Test equality
     img2 = ImageGrid(img)
     assert img2 == img
     assert (img2 != img) is False
 
-    # Test indexing round trip
-    np.testing.assert_allclose(img.ndcoords, img.ras(img.ndindex))
-    np.testing.assert_allclose(img.ndindex, np.round(img.index(img.ndcoords)))
-
 
 def test_ImageGrid_utils(tmpdir, testdata_path, get_testdata):
     """Check that images can be objects or paths and equality."""

nitransforms/tests/test_nonlinear.py

@@ -40,6 +40,30 @@ def test_displacements_init():
         )
 
 
+@pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", "oblique"])
+@pytest.mark.parametrize("axis", [0, 1, 2, (0, 1), (1, 2), (0, 1, 2)])
+def test_displacements_to_filename(tmp_path, get_testdata, image_orientation, axis):
+    """Exercise to_filename."""
+
+    nii = get_testdata[image_orientation]
+    fieldmap = np.zeros((*nii.shape[:3], 3), dtype="float32")
+    fieldmap[..., axis] = -10.0
+
+    xfm = DenseFieldTransform(
+        fieldmap,
+        reference=nii,
+    )
+    xfm.to_filename(tmp_path / "warp_itk.nii.gz", fmt="itk")
+    xfm.to_filename(tmp_path / "warp_afni.nii.gz", fmt="afni")
+    xfm.to_filename(tmp_path / "warp_fsl.nii.gz", fmt="fsl")
+
+    with pytest.raises(NotImplementedError):
+        xfm.to_filename(tmp_path / "warp_freesurfer.nii.gz", fmt="fs")
+
+    with pytest.raises(TypeError):
+        xfm.to_filename(tmp_path / "warp.x5", fmt="X5")
+
+
 @pytest.mark.parametrize("size", [(20, 20, 20), (20, 20, 20, 2, 3), (20, 20, 20, 1, 4)])
 def test_displacements_bad_sizes(size):
     """Checks field sizes."""

nitransforms/tests/test_resampling.py

@@ -149,10 +149,6 @@ def test_apply_linear_transform(
     assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL_LINEAR
 
 
-@pytest.mark.xfail(
-    reason="Disable while #266 is developed.",
-    strict=False,
-)
 @pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", "oblique"])
 @pytest.mark.parametrize("sw_tool", ["itk", "afni"])
 @pytest.mark.parametrize("axis", [0, 1, 2, (0, 1), (1, 2), (0, 1, 2)])
@@ -174,29 +170,24 @@ def test_displacements_field1(
     nii.to_filename("reference.nii.gz")
     msk.to_filename("mask.nii.gz")
 
-    fieldmap = np.zeros(
-        (*nii.shape[:3], 1, 3) if sw_tool != "fsl" else (*nii.shape[:3], 3),
-        dtype="float32",
-    )
+    fieldmap = np.zeros((*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")
-
+    # Generate a transform file for the particular software
     xfm_fname = "warp.nii.gz"
-    field = nb.Nifti1Image(fieldmap, nii.affine, _hdr)
-    field.to_filename(xfm_fname)
-
-    xfm = nitnl.load(xfm_fname, fmt=sw_tool)
+    xfm = nitnl.DenseFieldTransform(
+        fieldmap,
+        reference=nii,
+    )
+    xfm.to_filename(xfm_fname, fmt=sw_tool)
 
+    tool_output = tmp_path / f"{sw_tool}_brainmask.nii.gz"
     # Then apply the transform and cross-check with software
     cmd = APPLY_NONLINEAR_CMD[sw_tool](
         transform=os.path.abspath(xfm_fname),
         reference=tmp_path / "mask.nii.gz",
         moving=tmp_path / "mask.nii.gz",
-        output=tmp_path / "resampled_brainmask.nii.gz",
+        output=tool_output,
         extra="--output-data-type uchar" if sw_tool == "itk" else "",
     )
 
@@ -208,26 +199,28 @@ 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")
+    sw_moved_mask = np.asanyarray(nb.load(tool_output).dataobj, dtype=bool)
     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
+    nt_moved_mask.to_filename(tmp_path / "nit_brainmask.nii.gz")
     brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool)
+    percent_diff = (sw_moved_mask != brainmask)[5:-5, 5:-5, 5:-5].sum() / brainmask.size
+
+    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](
         transform=os.path.abspath(xfm_fname),
         reference=tmp_path / "reference.nii.gz",
         moving=tmp_path / "reference.nii.gz",
-        output=tmp_path / "resampled.nii.gz",
+        output=tmp_path / f"{sw_tool}_resampled.nii.gz",
         extra="--output-data-type uchar" if sw_tool == "itk" else "",
     )
 
     exit_code = check_call([cmd], shell=True)
     assert exit_code == 0
-    sw_moved = nb.load("resampled.nii.gz")
+    sw_moved = nb.load(f"{sw_tool}_resampled.nii.gz")
 
     nt_moved = apply(xfm, nii, order=0)
     nt_moved.set_data_dtype(nii.get_data_dtype())
@@ -240,10 +233,6 @@ def test_displacements_field1(
     assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL_LINEAR
 
 
-@pytest.mark.xfail(
-    reason="Disable while #266 is developed.",
-    strict=False,
-)
 @pytest.mark.parametrize("sw_tool", ["itk", "afni"])
 def test_displacements_field2(tmp_path, testdata_path, sw_tool):
     """Check a translation-only field on one or more axes, different image orientations."""