diff --git a/nitransforms/base.py b/nitransforms/base.py index eb6c2785..3f0e1151 100644 --- a/nitransforms/base.py +++ b/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.""" diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py index fe0b18d3..1f5b14f8 100644 --- a/nitransforms/nonlinear.py +++ b/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 {}) diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 6ade3eff..6166386d 100644 --- a/nitransforms/resampling.py +++ b/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( diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index fe9c8d20..c049a4ee 100644 --- a/nitransforms/tests/test_base.py +++ b/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.""" diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index a1142574..aeb03e73 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/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.""" diff --git a/nitransforms/tests/test_resampling.py b/nitransforms/tests/test_resampling.py index ca56a7e8..d7538ea4 100644 --- a/nitransforms/tests/test_resampling.py +++ b/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."""