diff --git a/env.yml b/env.yml index d550959..f4382b5 100644 --- a/env.yml +++ b/env.yml @@ -24,6 +24,8 @@ dependencies: - nitime=0.10 - scikit-image=0.22 - scikit-learn=1.4 + # SimpleITK, so build doesn't complain about building scikit from sources + - simpleitk=2.4 # Utilities - graphviz=9.0 - pandoc=3.1 diff --git a/nitransforms/io/itk.py b/nitransforms/io/itk.py index afabfd9..9960e94 100644 --- a/nitransforms/io/itk.py +++ b/nitransforms/io/itk.py @@ -1,4 +1,5 @@ """Read/write ITK transforms.""" + import warnings import numpy as np from scipy.io import loadmat as _read_mat, savemat as _save_mat @@ -138,8 +139,7 @@ def from_matlab_dict(cls, mdict, index=0): sa = tf.structarr affine = mdict.get( - "AffineTransform_double_3_3", - mdict.get("AffineTransform_float_3_3") + "AffineTransform_double_3_3", mdict.get("AffineTransform_float_3_3") ) if affine is None: @@ -337,7 +337,7 @@ def from_image(cls, imgobj): hdr = imgobj.header.copy() shape = hdr.get_data_shape() - if len(shape) != 5 or shape[-2] != 1 or not shape[-1] in (2, 3): + if len(shape) != 5 or shape[-2] != 1 or shape[-1] not in (2, 3): raise TransformFileError( 'Displacements field "%s" does not come from ITK.' % imgobj.file_map["image"].filename @@ -347,10 +347,9 @@ def from_image(cls, imgobj): warnings.warn("Incorrect intent identified.") hdr.set_intent("vector") - field = np.squeeze(np.asanyarray(imgobj.dataobj)) - field[..., (0, 1)] *= -1.0 + field = np.squeeze(np.asanyarray(imgobj.dataobj)).transpose(2, 1, 0, 3) - return imgobj.__class__(field, imgobj.affine, hdr) + return imgobj.__class__(field, LPS @ imgobj.affine, hdr) @classmethod def to_image(cls, imgobj): @@ -359,10 +358,8 @@ def to_image(cls, imgobj): hdr = imgobj.header.copy() hdr.set_intent("vector") - warp_data = imgobj.get_fdata().reshape(imgobj.shape[:3] + (1, imgobj.shape[-1])) - warp_data[..., (0, 1)] *= -1 - - return imgobj.__class__(warp_data, imgobj.affine, hdr) + warp_data = imgobj.get_fdata().transpose(2, 1, 0, 3)[..., None, :] + return imgobj.__class__(warp_data, LPS @ imgobj.affine, hdr) class ITKCompositeH5: @@ -410,21 +407,16 @@ def from_h5obj(cls, fileobj, check=True, only_linear=False): directions = np.reshape(_fixed[9:], (3, 3)) affine = from_matvec(directions * zooms, offset) # ITK uses Fortran ordering, like NIfTI, but with the vector dimension first - field = np.moveaxis( - np.reshape( - xfm[f"{typo_fallback}Parameters"], (3, *shape.astype(int)), order='F' - ), - 0, - -1, - ) - field[..., (0, 1)] *= -1.0 + # In practice, this seems to work (see issue #171) + field = np.reshape( + xfm[f"{typo_fallback}Parameters"], (*shape.astype(int), 3) + ).transpose(2, 1, 0, 3) + hdr = Nifti1Header() hdr.set_intent("vector") hdr.set_data_dtype("float") - xfm_list.append( - Nifti1Image(field.astype("float"), LPS @ affine, hdr) - ) + xfm_list.append(Nifti1Image(field.astype("float"), affine, hdr)) continue raise TransformIOError( diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py index 24e043c..6b1d942 100644 --- a/nitransforms/nonlinear.py +++ b/nitransforms/nonlinear.py @@ -285,7 +285,7 @@ def to_x5(self, metadata=None): ) @classmethod - def from_filename(cls, filename, fmt="X5", x5_position=0): + def from_filename(cls, filename, is_deltas=True, fmt="X5", x5_position=0): _factory = { "afni": io.afni.AFNIDisplacementsField, "itk": io.itk.ITKDisplacementsField, @@ -299,7 +299,7 @@ def from_filename(cls, filename, fmt="X5", x5_position=0): if fmt == "X5": return from_x5(load_x5(filename), x5_position=x5_position) - return cls(_factory[fmt.lower()].from_filename(filename)) + return cls(_factory[fmt.lower()].from_filename(filename), is_deltas=is_deltas) load = DenseFieldTransform.from_filename diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 98ef445..7b102e6 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -254,22 +254,22 @@ def apply( targets = None ref_ndcoords = _ref.ndcoords.T - if hasattr(transform, "to_field") and callable(transform.to_field): - targets = ImageGrid(spatialimage).index( - _as_homogeneous( - transform.to_field(reference=reference).map(ref_ndcoords), - dim=_ref.ndim, - ) - ) - else: - # Targets' shape is (Nt, 3, Nv) with Nv = Num. voxels, Nt = Num. timepoints. - targets = ( - ImageGrid(spatialimage).index( - _as_homogeneous(transform.map(ref_ndcoords), dim=_ref.ndim) - ) - if targets is None - else targets + # if hasattr(transform, "to_field") and callable(transform.to_field): + # targets = ImageGrid(spatialimage).index( + # _as_homogeneous( + # transform.to_field(reference=reference).map(ref_ndcoords), + # dim=_ref.ndim, + # ) + # ) + # else: + # Targets' shape is (Nt, 3, Nv) with Nv = Num. voxels, Nt = Num. timepoints. + targets = ( + ImageGrid(spatialimage).index( + _as_homogeneous(transform.map(ref_ndcoords), dim=_ref.ndim) ) + if targets is None + else targets + ) if targets.ndim == 3: targets = np.rollaxis(targets, targets.ndim - 1, 0) diff --git a/nitransforms/tests/test_io.py b/nitransforms/tests/test_io.py index 7058cdc..d6dfefe 100644 --- a/nitransforms/tests/test_io.py +++ b/nitransforms/tests/test_io.py @@ -1,6 +1,7 @@ # emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: """I/O test cases.""" + import os from subprocess import check_call from io import StringIO @@ -10,11 +11,13 @@ import pytest from h5py import File as H5File +import SimpleITK as sitk import nibabel as nb from nibabel.eulerangles import euler2mat from nibabel.affines import from_matvec from scipy.io import loadmat from nitransforms.linear import Affine +from nitransforms import nonlinear as nitnl from nitransforms.io import ( afni, fsl, @@ -68,11 +71,13 @@ def test_volume_group_voxel_ordering(): def test_VG_from_LTA(data_path): """Check the affine interpolation from volume geometries.""" # affine manually clipped after running mri_info on the image - oracle = np.loadtxt(StringIO("""\ + oracle = np.loadtxt( + StringIO("""\ -3.0000 0.0000 -0.0000 91.3027 -0.0000 2.0575 -2.9111 -25.5251 0.0000 2.1833 2.7433 -105.0820 - 0.0000 0.0000 0.0000 1.0000""")) + 0.0000 0.0000 0.0000 1.0000""") + ) lta_text = "\n".join( (data_path / "bold-to-t1w.lta").read_text().splitlines()[13:21] @@ -419,10 +424,17 @@ def test_afni_Displacements(): @pytest.mark.parametrize("only_linear", [True, False]) -@pytest.mark.parametrize("h5_path,nxforms", [ - (_datadir / "affine-antsComposite.h5", 1), - (_testdir / "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5", 2), -]) +@pytest.mark.parametrize( + "h5_path,nxforms", + [ + (_datadir / "affine-antsComposite.h5", 1), + ( + _testdir + / "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5", + 2, + ), + ], +) def test_itk_h5(tmpdir, only_linear, h5_path, nxforms): """Test displacements fields.""" assert ( @@ -434,7 +446,9 @@ def test_itk_h5(tmpdir, only_linear, h5_path, nxforms): ) ) ) - == nxforms if not only_linear else 1 + == nxforms + if not only_linear + else 1 ) with pytest.raises(TransformFileError): @@ -465,24 +479,33 @@ def test_regressions(file_type, test_file, data_path): file_type.from_filename(data_path / "regressions" / test_file) -@pytest.mark.parametrize("parameters", [ - {"x": 0.1, "y": 0.03, "z": 0.002}, - {"x": 0.001, "y": 0.3, "z": 0.002}, - {"x": 0.01, "y": 0.03, "z": 0.2}, -]) +@pytest.mark.parametrize( + "parameters", + [ + {"x": 0.1, "y": 0.03, "z": 0.002}, + {"x": 0.001, "y": 0.3, "z": 0.002}, + {"x": 0.01, "y": 0.03, "z": 0.2}, + ], +) @pytest.mark.parametrize("dir_x", (-1, 1)) @pytest.mark.parametrize("dir_y", (-1, 1)) @pytest.mark.parametrize("dir_z", (1, -1)) -@pytest.mark.parametrize("swapaxes", [ - None, (0, 1), (1, 2), (0, 2), -]) +@pytest.mark.parametrize( + "swapaxes", + [ + None, + (0, 1), + (1, 2), + (0, 2), + ], +) def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, dir_z): tmpdir.chdir() img, R = _generate_reoriented( testdata_path / "someones_anatomy.nii.gz", (dir_x, dir_y, dir_z), swapaxes, - parameters + parameters, ) img.to_filename("orig.nii.gz") @@ -507,9 +530,8 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, "orig.nii.gz", ) - diff = ( - np.asanyarray(img.dataobj, dtype="uint8") - - np.asanyarray(nt3drefit.dataobj, dtype="uint8") + diff = np.asanyarray(img.dataobj, dtype="uint8") - np.asanyarray( + nt3drefit.dataobj, dtype="uint8" ) assert np.sqrt((diff[10:-10, 10:-10, 10:-10] ** 2).mean()) < 0.1 @@ -522,14 +544,15 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, "deob_3drefit.nii.gz", ) - diff = ( - np.asanyarray(img.dataobj, dtype="uint8") - - np.asanyarray(nt_undo3drefit.dataobj, dtype="uint8") + diff = np.asanyarray(img.dataobj, dtype="uint8") - np.asanyarray( + nt_undo3drefit.dataobj, dtype="uint8" ) assert np.sqrt((diff[10:-10, 10:-10, 10:-10] ** 2).mean()) < 0.1 # Check the target grid by 3dWarp and the affine & size interpolated by NiTransforms - cmd = f"3dWarp -verb -deoblique -NN -prefix {tmpdir}/deob.nii.gz {tmpdir}/orig.nii.gz" + cmd = ( + f"3dWarp -verb -deoblique -NN -prefix {tmpdir}/deob.nii.gz {tmpdir}/orig.nii.gz" + ) assert check_call([cmd], shell=True) == 0 deobnii = nb.load("deob.nii.gz") @@ -540,11 +563,12 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, # Check resampling in deobliqued grid ntdeobnii = apply( - Affine(np.eye(4), reference=deobnii.__class__( - np.zeros(deobshape, dtype="uint8"), - deobaff, - deobnii.header - )), + Affine( + np.eye(4), + reference=deobnii.__class__( + np.zeros(deobshape, dtype="uint8"), deobaff, deobnii.header + ), + ), img, order=0, ) @@ -559,9 +583,8 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, ) mask = np.asanyarray(ntdeobmask.dataobj, dtype=bool) - diff = ( - np.asanyarray(deobnii.dataobj, dtype="uint8") - - np.asanyarray(ntdeobnii.dataobj, dtype="uint8") + diff = np.asanyarray(deobnii.dataobj, dtype="uint8") - np.asanyarray( + ntdeobnii.dataobj, dtype="uint8" ) assert np.sqrt((diff[mask] ** 2).mean()) < 0.1 @@ -591,7 +614,7 @@ def _generate_reoriented(path, directions, swapaxes, parameters): aff = np.diag((*directions, 1)) @ aff for ax in range(3): - if (directions[ax] == -1): + if directions[ax] == -1: aff[ax, 3] = last_xyz[ax] data = np.flip(data, ax) @@ -621,16 +644,15 @@ def test_itk_linear_h5(tmpdir, data_path, testdata_path): assert len(h5xfm.xforms) == 1 # File loadable with single affine object - itk.ITKLinearTransform.from_filename( - data_path / "affine-antsComposite.h5" - ) + itk.ITKLinearTransform.from_filename(data_path / "affine-antsComposite.h5") with open(data_path / "affine-antsComposite.h5", "rb") as f: itk.ITKLinearTransform.from_fileobj(f) # Exercise only_linear itk.ITKCompositeH5.from_filename( - testdata_path / "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5", + testdata_path + / "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5", only_linear=True, ) @@ -673,55 +695,66 @@ def test_itk_linear_h5(tmpdir, data_path, testdata_path): with pytest.raises(TransformIOError): itk.ITKLinearTransform.from_filename("test.h5") -# Added tests for h5 orientation bug +# 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): + """Exercise I/O of ITK displacements fields.""" + nii = get_testdata[image_orientation] -@pytest.mark.xfail( - reason="GH-137/GH-171: displacement field dimension order is wrong", - strict=False, -) -def test_itk_h5_field_order(tmp_path): - """Displacement fields stored in row-major order should fail to round-trip.""" - shape = (3, 4, 5) - vals = np.arange(np.prod(shape), dtype=float).reshape(shape) - field = np.stack([vals, vals + 100, vals + 200], axis=0) - - params = field.reshape(-1, order="C") - fixed = np.array(list(shape) + [0, 0, 0] + [1, 1, 1] + list(np.eye(3).ravel()), dtype=float) - fname = tmp_path / "field.h5" - with H5File(fname, "w") as f: - grp = f.create_group("TransformGroup") - grp.create_group("0")["TransformType"] = np.array([b"CompositeTransform_double_3_3"]) - g1 = grp.create_group("1") - g1["TransformType"] = np.array([b"DisplacementFieldTransform_float_3_3"]) - g1["TransformFixedParameters"] = fixed - g1["TransformParameters"] = params - - img = itk.ITKCompositeH5.from_filename(fname)[0] - expected = np.moveaxis(field, 0, -1) - expected[..., (0, 1)] *= -1 - assert np.allclose(img.get_fdata(), expected) + # Create a reference centered at the origin with various axis orders/flips + 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, )) + + nit_nii = itk.ITKDisplacementsField.to_image( + nb.Nifti1Image(field, ref_affine, None) + ) + + itk_file = tmp_path / "itk_displacements.nii.gz" + itk_img = sitk.GetImageFromArray(field, isVector=True) + itk_img.SetOrigin(tuple(ref_affine[:3, 3])) + zooms = np.sqrt((ref_affine[:3, :3] ** 2).sum(0)) + itk_img.SetSpacing(tuple(zooms)) + direction = (ref_affine[:3, :3] / zooms).ravel() + itk_img.SetDirection(tuple(direction)) + sitk.WriteImage(itk_img, str(itk_file)) + itk_nit_nii = itk.ITKDisplacementsField.from_filename(itk_file) + assert itk_nit_nii.shape == field.shape + np.testing.assert_allclose(itk_nit_nii.affine, ref_affine) + np.testing.assert_allclose(itk_nit_nii.dataobj, field) + + itk_nii = nb.load(itk_file) + assert nit_nii.shape == itk_nii.shape + np.testing.assert_allclose(itk_nii.dataobj, nit_nii.dataobj) + np.testing.assert_allclose(itk_nii.affine, nit_nii.affine) + + +# Added tests for h5 orientation bug (#167) def _load_composite_testdata(data_path): """Return the composite HDF5 and displacement field from regressions.""" h5file = data_path / "regressions" / "ants_t1_to_mniComposite.h5" # Generated using # CompositeTransformUtil --disassemble ants_t1_to_mniComposite.h5 \ # ants_t1_to_mniComposite - warpfile = data_path / "regressions" / ( - "01_ants_t1_to_mniComposite_DisplacementFieldTransform.nii.gz" + warpfile = ( + data_path + / "regressions" + / ("01_ants_t1_to_mniComposite_DisplacementFieldTransform.nii.gz") ) if not (h5file.exists() and warpfile.exists()): pytest.skip("Composite transform test data not available") return h5file, warpfile -@pytest.mark.xfail( - reason="GH-137/GH-171: displacement field dimension order is wrong", - strict=False, -) -def test_itk_h5_displacement_mismatch(testdata_path): +def test_itk_h5_and_displacement_equivalence(testdata_path): """Composite displacements should match the standalone field""" h5file, warpfile = _load_composite_testdata(testdata_path) xforms = itk.ITKCompositeH5.from_filename(h5file) @@ -732,49 +765,21 @@ def test_itk_h5_displacement_mismatch(testdata_path): np.asanyarray(field_h5.dataobj), np.asanyarray(field_img.dataobj) ) + 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, 0.0, -11.0], + ] + ) -def test_itk_h5_transpose_fix(testdata_path): - """Check the displacement field orientation explicitly. - - ITK stores displacement fields with the vector dimension leading in - Fortran (column-major) order [1]_. Transposing the parameters from the HDF5 - composite file accordingly should match the standalone displacement image. + # Load the displacements field + xfm_h5 = nitnl.DenseFieldTransform(xforms[1]) + mapped_h5 = xfm_h5.map(points) - References - ---------- - .. [1] ITK Software Guide. https://itk.org/ItkSoftwareGuide.pdf - """ - h5file, warpfile = _load_composite_testdata(testdata_path) + xfm_warp = nitnl.DenseFieldTransform(field_img) + mapped_warp = xfm_warp.map(points) - with H5File(h5file, "r") as f: - group = f["TransformGroup"]["2"] - size = group["TransformFixedParameters"][:3].astype(int) - params = group["TransformParameters"][:].reshape(*size, 3) - - img = nb.load(warpfile) - ref = np.squeeze(np.asanyarray(img.dataobj)) - - np.testing.assert_array_equal(params.transpose(2, 1, 0, 3), ref) - - -def test_itk_h5_field_order_fortran(tmp_path): - """Verify Fortran-order displacement fields load correctly""" - shape = (3, 4, 5) - vals = np.arange(np.prod(shape), dtype=float).reshape(shape) - field = np.stack([vals, vals + 100, vals + 200], axis=0) - - params = field.reshape(-1, order="F") - fixed = np.array(list(shape) + [0, 0, 0] + [1, 1, 1] + list(np.eye(3).ravel()), dtype=float) - fname = tmp_path / "field_f.h5" - with H5File(fname, "w") as f: - grp = f.create_group("TransformGroup") - grp.create_group("0")["TransformType"] = np.array([b"CompositeTransform_double_3_3"]) - g1 = grp.create_group("1") - g1["TransformType"] = np.array([b"DisplacementFieldTransform_float_3_3"]) - g1["TransformFixedParameters"] = fixed - g1["TransformParameters"] = params - - img = itk.ITKCompositeH5.from_filename(fname)[0] - expected = np.moveaxis(field, 0, -1) - expected[..., (0, 1)] *= -1 - assert np.allclose(img.get_fdata(), expected) + np.testing.assert_array_equal(mapped_h5, mapped_warp) diff --git a/nitransforms/tests/test_manip.py b/nitransforms/tests/test_manip.py index 4d7b3a5..73344df 100644 --- a/nitransforms/tests/test_manip.py +++ b/nitransforms/tests/test_manip.py @@ -2,6 +2,9 @@ # vi: set ft=python sts=4 ts=4 sw=4 et: """Tests of nonlinear transforms.""" +from subprocess import check_call +import shutil + import pytest import numpy as np @@ -11,7 +14,7 @@ from ..manip import TransformChain from ..linear import Affine from ..nonlinear import DenseFieldTransform -from ..io import x5 +from ..io import x5, itk FMT = {"lta": "fs", "tfm": "itk"} @@ -86,3 +89,42 @@ def test_transformchain_x5_roundtrip(tmp_path): assert len(loaded0) == len(chain) assert len(loaded1) == len(chain) assert np.allclose(chain.map([[0, 0, 0]]), loaded1.map([[0, 0, 0]])) + + +@pytest.mark.xfail( + reason="TransformChain.map() doesn't work properly (discovered with GH-167)", + strict=False, +) +def test_composite_h5_map_against_ants(testdata_path, tmp_path): + """Map points with NiTransforms and compare to ANTs.""" + h5file = testdata_path / "regressions" / "ants_t1_to_mniComposite.h5" + if not h5file.exists(): + pytest.skip(f"Necessary file <{h5file}> missing.") + + 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, 0.0, -11.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 {h5file}" + 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 + + xforms = itk.ITKCompositeH5.from_filename(h5file) + chain = Affine(xforms[0].to_ras()) + DenseFieldTransform(xforms[1]) + mapped = chain.map(points) + + assert np.allclose(mapped, ants_pts, atol=1e-6) diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index 936a62f..dc45cf9 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/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)]) @@ -243,3 +248,117 @@ 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.linspace(-50, 50, num=np.prod(shape)), + np.linspace(-80, 80, num=np.prod(shape)), + np.zeros(np.prod(shape)) + )).reshape(shape + (3, )) + fieldnii = nb.Nifti1Image(field, ref_affine, None) + + warpfile = tmp_path / "const_disp.nii.gz" + ITKDisplacementsField.to_filename(fieldnii, warpfile) + xfm = DenseFieldTransform(fieldnii) + xfm2 = DenseFieldTransform(ITKDisplacementsField.from_filename(warpfile)) + + np.testing.assert_allclose(xfm.reference.affine, xfm2.reference.affine) + + 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, 0.0, -11.0], + ] + ) + + if gridpoints: + coords = xfm.reference.ndcoords + points = (ref_affine @ np.vstack((coords, np.ones((1, coords.shape[1]))))).T[:, :3] + + 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 + + import pdb; pdb.set_trace() + + 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]}" + + mapped = xfm.map(points) + np.testing.assert_array_equal(np.round(mapped, 3), 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]}" diff --git a/nitransforms/tests/test_resampling.py b/nitransforms/tests/test_resampling.py index b65bf57..c0f2998 100644 --- a/nitransforms/tests/test_resampling.py +++ b/nitransforms/tests/test_resampling.py @@ -177,16 +177,21 @@ def test_displacements_field1( fieldmap[..., axis] = -10.0 _hdr = nii.header.copy() + affine = nii.affine.copy() if sw_tool in ("itk",): _hdr.set_intent("vector") + affine = io.itk.LPS @ affine _hdr.set_data_dtype("float32") + + field = nb.Nifti1Image(fieldmap, affine, _hdr) 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) + np.testing.assert_array_equal(xfm._deltas, np.squeeze(field.dataobj)) + # Then apply the transform and cross-check with software cmd = APPLY_NONLINEAR_CMD[sw_tool]( transform=os.path.abspath(xfm_fname), @@ -226,12 +231,9 @@ def test_displacements_field1( sw_moved = nb.load("resampled.nii.gz") nt_moved = apply(xfm, nii, order=0) - nt_moved.set_data_dtype(nii.get_data_dtype()) 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()) + diff = sw_moved.get_fdata() - nt_moved.get_fdata() + # A certain tolerance is necessary because of resampling at borders assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL_LINEAR diff --git a/pyproject.toml b/pyproject.toml index c4a1b8e..0aa12f7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,6 +41,8 @@ test = [ "pytest-xdist >= 2.5", "coverage[toml] >= 5.2.1", "nitransforms[niftiext]", + "SimpleITK ~= 2.4", + "scikit-build", ] # Aliases niftiexts = ["nitransforms[niftiext]"]