From d624b2aed2d7de331a737cf25105c289d395dad9 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 24 Jul 2025 16:50:53 +0200 Subject: [PATCH 1/6] fix: proper ordering of displacements fields when reading and writing Transposes data as suggested by @feilong when reading and writing displacements fields. Tests have been added to ensure round-trip consistency. Related-to: #171. Co-Authored-By: Feilong Ma --- env.yml | 2 ++ nitransforms/io/itk.py | 28 +++++++++-------------- nitransforms/tests/test_resampling.py | 33 ++++++++++++++++----------- pyproject.toml | 2 ++ 4 files changed, 35 insertions(+), 30 deletions(-) diff --git a/env.yml b/env.yml index d550959b..f4382b52 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 b2a6df8a..1be51854 100644 --- a/nitransforms/io/itk.py +++ b/nitransforms/io/itk.py @@ -347,10 +347,10 @@ def from_image(cls, imgobj): warnings.warn("Incorrect intent identified.") hdr.set_intent("vector") - field = np.squeeze(np.asanyarray(imgobj.dataobj)) + field = np.squeeze(np.asanyarray(imgobj.dataobj)).transpose(2, 1, 0, 3) field[..., (0, 1)] *= -1.0 - return imgobj.__class__(field, imgobj.affine, hdr) + return imgobj.__class__(field, LPS @ imgobj.affine, hdr) @classmethod def to_image(cls, imgobj): @@ -359,10 +359,9 @@ 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) + field = imgobj.get_fdata().transpose(2, 1, 0, 3)[..., None, :] + field[..., (0, 1)] *= -1.0 + return imgobj.__class__(field, LPS @ imgobj.affine, hdr) class ITKCompositeH5: @@ -410,21 +409,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/tests/test_resampling.py b/nitransforms/tests/test_resampling.py index ca56a7e8..b2c75f52 100644 --- a/nitransforms/tests/test_resampling.py +++ b/nitransforms/tests/test_resampling.py @@ -175,19 +175,20 @@ def test_displacements_field1( msk.to_filename("mask.nii.gz") fieldmap = np.zeros( - (*nii.shape[:3], 1, 3) if sw_tool != "fsl" else (*nii.shape[:3], 3), + (*nii.shape[:3], 1, 3) if sw_tool == "afni" else (*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") - xfm_fname = "warp.nii.gz" field = nb.Nifti1Image(fieldmap, nii.affine, _hdr) - field.to_filename(xfm_fname) + + xfm_fname = "warp.nii.gz" + if sw_tool == "itk": + io.itk.ITKDisplacementsField.to_filename(field, xfm_fname) + else: + field.to_filename(xfm_fname) xfm = nitnl.load(xfm_fname, fmt=sw_tool) @@ -197,7 +198,8 @@ def test_displacements_field1( reference=tmp_path / "mask.nii.gz", moving=tmp_path / "mask.nii.gz", output=tmp_path / "resampled_brainmask.nii.gz", - extra="--output-data-type uchar" if sw_tool == "itk" else "", + extra="", + # extra="--output-data-type uchar" if sw_tool == "itk" else "", ) # skip test if command is not available on host @@ -207,14 +209,18 @@ 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") 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 + # Calculate xor between both: + sw_mask = np.asanyarray(sw_moved_mask.dataobj, dtype=bool) brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool) + percent_diff = (sw_mask != brainmask)[5:-5, 5:-5, 5:-5].sum() / brainmask.size + + assert exit_code == 0 + 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]( @@ -236,8 +242,9 @@ def test_displacements_field1( diff = np.asanyarray( sw_moved.dataobj, dtype=sw_moved.get_data_dtype() ) - np.asanyarray(nt_moved.dataobj, dtype=nt_moved.get_data_dtype()) - # A certain tolerance is necessary because of resampling at borders - assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL_LINEAR + + # Drop samples close to the border of the image + assert np.sqrt((diff[5:-5, 5:-5, 5:-5] ** 2).mean()) < 1e-6 @pytest.mark.xfail( diff --git a/pyproject.toml b/pyproject.toml index c4a1b8e4..0aa12f7a 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]"] From b2a794bf972c8f2656244fb8291266eb67e8b978 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 25 Jul 2025 06:39:26 +0200 Subject: [PATCH 2/6] fix: remove sign flip as it seems incorrect --- nitransforms/io/itk.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/nitransforms/io/itk.py b/nitransforms/io/itk.py index 1be51854..316a5673 100644 --- a/nitransforms/io/itk.py +++ b/nitransforms/io/itk.py @@ -348,8 +348,6 @@ def from_image(cls, imgobj): hdr.set_intent("vector") field = np.squeeze(np.asanyarray(imgobj.dataobj)).transpose(2, 1, 0, 3) - field[..., (0, 1)] *= -1.0 - return imgobj.__class__(field, LPS @ imgobj.affine, hdr) @classmethod @@ -360,7 +358,6 @@ def to_image(cls, imgobj): hdr.set_intent("vector") field = imgobj.get_fdata().transpose(2, 1, 0, 3)[..., None, :] - field[..., (0, 1)] *= -1.0 return imgobj.__class__(field, LPS @ imgobj.affine, hdr) From ae6c3309c3f78d6580e615aa8bf6f29d8c360eee Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 25 Jul 2025 08:27:24 +0200 Subject: [PATCH 3/6] enh: update tests accordingly --- nitransforms/io/itk.py | 8 +- nitransforms/tests/test_nonlinear.py | 128 +++++++++++++++++++++++++++ 2 files changed, 134 insertions(+), 2 deletions(-) diff --git a/nitransforms/io/itk.py b/nitransforms/io/itk.py index 316a5673..c67ecbbe 100644 --- a/nitransforms/io/itk.py +++ b/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) diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index a1142574..f50341a8 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/nitransforms/tests/test_nonlinear.py @@ -3,6 +3,9 @@ """Tests of nonlinear transforms.""" import os +from subprocess import check_call +import shutil + import pytest import numpy as np @@ -234,3 +237,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]}" From 6a365771c0ea6390183b638be57e16bc978eb29a Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 25 Jul 2025 18:11:03 +0200 Subject: [PATCH 4/6] fix: write tests to avoid regressions Resolves: #267. --- nitransforms/tests/test_nonlinear.py | 11 ++++++ nitransforms/tests/test_resampling.py | 51 ++++++++++++++++++++++----- 2 files changed, 53 insertions(+), 9 deletions(-) diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index f50341a8..8de59446 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/nitransforms/tests/test_nonlinear.py @@ -20,6 +20,17 @@ rng = np.random.default_rng() +SOME_TEST_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_displacements_init(): identity1 = DenseFieldTransform( np.zeros((10, 10, 10, 3)), diff --git a/nitransforms/tests/test_resampling.py b/nitransforms/tests/test_resampling.py index b2c75f52..9056b395 100644 --- a/nitransforms/tests/test_resampling.py +++ b/nitransforms/tests/test_resampling.py @@ -150,13 +150,13 @@ def test_apply_linear_transform( @pytest.mark.xfail( - reason="Disable while #266 is developed.", + reason="GH-267: disabled while debugging", 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)]) -def test_displacements_field1( +def test_apply_displacements_field1( tmp_path, get_testdata, get_testmask, @@ -190,16 +190,17 @@ def test_displacements_field1( else: field.to_filename(xfm_fname) - xfm = nitnl.load(xfm_fname, fmt=sw_tool) + # xfm = nitnl.load(xfm_fname, fmt=sw_tool) + xfm = nitnl.DenseFieldTransform(fieldmap, reference=nii) + ants_output = tmp_path / "ants_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", - extra="", - # extra="--output-data-type uchar" if sw_tool == "itk" else "", + output=ants_output, + extra="--output-data-type uchar" if sw_tool == "itk" else "", ) # skip test if command is not available on host @@ -209,11 +210,13 @@ def test_displacements_field1( # resample mask exit_code = check_call([cmd], shell=True) - sw_moved_mask = nb.load("resampled_brainmask.nii.gz") + assert exit_code == 0 + sw_moved_mask = nb.load(ants_output) nt_moved_mask = apply(xfm, msk, order=0) - # Calculate xor between both: - sw_mask = np.asanyarray(sw_moved_mask.dataobj, dtype=bool) + nt_moved_mask.to_filename(tmp_path / "nit_brainmask.nii.gz") + + assert np.sqrt((diff**2).mean()) < RMSE_TOL_LINEAR brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool) percent_diff = (sw_mask != brainmask)[5:-5, 5:-5, 5:-5].sum() / brainmask.size @@ -403,3 +406,33 @@ def test_apply_4d(serialize_4d): data = np.asanyarray(moved.dataobj) idxs = [tuple(np.argwhere(data[..., i])[0]) for i in range(nvols)] assert idxs == [(9 - i, 2, 2) for i in range(nvols)] + + +@pytest.mark.xfail( + reason="GH-267: disabled while debugging", + strict=False, +) +def test_apply_bspline(tmp_path, testdata_path): + """Cross-check B-Splines and deformation field.""" + os.chdir(str(tmp_path)) + + img_name = testdata_path / "someones_anatomy.nii.gz" + disp_name = testdata_path / "someones_displacement_field.nii.gz" + bs_name = testdata_path / "someones_bspline_coefficients.nii.gz" + + bsplxfm = nitnl.BSplineFieldTransform(bs_name, reference=img_name) + dispxfm = nitnl.DenseFieldTransform(disp_name) + + out_disp = apply(dispxfm, img_name) + out_bspl = apply(bsplxfm, img_name) + + out_disp.to_filename("resampled_field.nii.gz") + out_bspl.to_filename("resampled_bsplines.nii.gz") + + assert ( + np.sqrt( + (out_disp.get_fdata(dtype="float32") - out_bspl.get_fdata(dtype="float32")) + ** 2 + ).mean() + < 0.2 + ) From c033fce2f2de0eca566b386a6afb532fc1dea25a Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 25 Jul 2025 20:16:13 +0200 Subject: [PATCH 5/6] fix: make nitransforms/tests/test_nonlinear.py::test_densefield_map_vs_ants pass --- nitransforms/io/itk.py | 8 +- nitransforms/tests/test_nonlinear.py | 297 +++++++++++++++------------ 2 files changed, 165 insertions(+), 140 deletions(-) diff --git a/nitransforms/io/itk.py b/nitransforms/io/itk.py index c67ecbbe..d860009e 100644 --- a/nitransforms/io/itk.py +++ b/nitransforms/io/itk.py @@ -348,9 +348,11 @@ def from_image(cls, imgobj): hdr.set_intent("vector") 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) + 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) @classmethod def to_image(cls, imgobj): diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index 8de59446..2029dbbf 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/nitransforms/tests/test_nonlinear.py @@ -16,21 +16,11 @@ DenseFieldTransform, ) from nitransforms.tests.utils import get_points +from nitransforms.io.itk import ITKDisplacementsField rng = np.random.default_rng() -SOME_TEST_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_displacements_init(): identity1 = DenseFieldTransform( np.zeros((10, 10, 10, 3)), @@ -80,7 +70,6 @@ def test_bsplines_references(testdata_path): reference=testdata_path / "someones_anatomy.nii.gz", ) - @pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", "oblique"]) @pytest.mark.parametrize("ongrid", [True, False]) def test_densefield_map(get_testdata, image_orientation, ongrid): @@ -173,6 +162,165 @@ def test_densefield_map_vs_bspline(tmp_path, testdata_path): np.testing.assert_allclose(dispxfm._field, bsplxfm._field, atol=1e-1, rtol=1e-4) +@pytest.mark.parametrize("ongrid", [True, False]) +def test_densefield_map_vs_ants(testdata_path, tmp_path, ongrid): + """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") + + nii = ITKDisplacementsField.from_filename(warpfile) + + # Get sampling indices + coords_xyz, points_ijk, grid_xyz, shape, ref_affine, reference, subsample = ( + get_points(nii, ongrid, npoints=5, rng=rng) + ) + coords_map = grid_xyz.reshape(*shape, 3) + + csvin = tmp_path / "fixed_coords.csv" + csvout = tmp_path / "moving_coords.csv" + np.savetxt(csvin, coords_xyz, delimiter=",", header="x,y,z", comments="") + + 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(nii, reference=reference) + mapped = xfm.map(coords_xyz) + + if ongrid: + ants_mapped_xyz = ants_pts.reshape(*shape, 3) + nit_mapped_xyz = mapped.reshape(*shape, 3) + + nb.Nifti1Image(coords_map, ref_affine, None).to_filename( + tmp_path / "baseline_field.nii.gz" + ) + + nb.Nifti1Image(ants_mapped_xyz, ref_affine, None).to_filename( + tmp_path / "ants_deformation_xyz.nii.gz" + ) + nb.Nifti1Image(nit_mapped_xyz, ref_affine, None).to_filename( + tmp_path / "nit_deformation_xyz.nii.gz" + ) + nb.Nifti1Image(ants_mapped_xyz - coords_map, ref_affine, None).to_filename( + tmp_path / "ants_deltas_xyz.nii.gz" + ) + nb.Nifti1Image(nit_mapped_xyz - coords_map, ref_affine, None).to_filename( + tmp_path / "nit_deltas_xyz.nii.gz" + ) + + atol = 0 if ongrid else 1e-2 + rtol = 1e-4 if ongrid else 1e-6 + assert np.allclose(mapped, ants_pts, atol=atol, rtol=rtol) + + +@pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", "oblique"]) +@pytest.mark.parametrize("ongrid", [True, False]) +def test_constant_field_vs_ants(tmp_path, get_testdata, image_orientation, ongrid): + """Create a constant displacement field and compare mappings.""" + + nii = get_testdata[image_orientation] + + # Get sampling indices + coords_xyz, points_ijk, grid_xyz, shape, ref_affine, reference, subsample = ( + get_points(nii, ongrid, npoints=5, rng=rng) + ) + + coords_map = grid_xyz.reshape(*shape, 3) + gold_mapped_xyz = coords_map + deltas + + deltas = 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(deltas, 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) + + # Ensure transform (xfm_orig) and ITK roundtrip (itk_xfm) are equivalent + xfm_orig = DenseFieldTransform(deltas, reference=reference) + np.testing.assert_allclose(xfm_orig.reference.shape, itk_xfm.reference.shape) + np.testing.assert_allclose(ref_affine, xfm_orig.reference.affine) + np.testing.assert_allclose(xfm_orig.reference.affine, itk_xfm.reference.affine) + np.testing.assert_allclose(xfm_orig._field, itk_xfm._field) + + # Ensure deltas and mapped grid are equivalent + grid_mapped_xyz = itk_xfm.map(grid_xyz).reshape(*shape, -1) + orig_grid_mapped_xyz = xfm_orig.map(grid_xyz).reshape(*shape, -1) + + # Check apparent healthiness of mapping + np.testing.assert_array_equal(orig_grid_mapped_xyz, grid_mapped_xyz) + np.testing.assert_array_equal(gold_mapped_xyz, orig_grid_mapped_xyz) + np.testing.assert_array_equal(gold_mapped_xyz, grid_mapped_xyz) + + csvout = tmp_path / "mapped_xyz.csv" + csvin = tmp_path / "coords_xyz.csv" + np.savetxt(csvin, coords_xyz, delimiter=",", header="x,y,z", comments="") + + 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 + + nb.Nifti1Image(grid_mapped_xyz, ref_affine, None).to_filename( + tmp_path / "grid_mapped.nii.gz" + ) + nb.Nifti1Image(coords_map, ref_affine, None).to_filename( + tmp_path / "baseline_field.nii.gz" + ) + nb.Nifti1Image(gold_mapped_xyz, ref_affine, None).to_filename( + tmp_path / "gold_mapped_xyz.nii.gz" + ) + + if ongrid: + ants_pts = ants_pts.reshape(*shape, 3) + + nb.Nifti1Image(ants_pts, ref_affine, None).to_filename( + tmp_path / "ants_mapped_xyz.nii.gz" + ) + np.testing.assert_array_equal(gold_mapped_xyz, ants_pts) + np.testing.assert_array_equal(deltas, ants_pts - coords_map) + else: + ants_deltas = ants_pts - coords_xyz + deltas_xyz = deltas.reshape(-1, 3)[subsample] + gold_xyz = coords_xyz + deltas_xyz + np.testing.assert_array_equal(gold_xyz, ants_pts) + np.testing.assert_array_equal(deltas_xyz, 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]}" + + @pytest.mark.parametrize("is_deltas", [True, False]) def test_densefield_oob_resampling(is_deltas): """Ensure mapping outside the field returns input coordinates.""" @@ -248,128 +396,3 @@ 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]}" From 6d32a6fbe4993b2bc72e570be6c9cb68d765e137 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 8 Aug 2025 16:34:22 +0200 Subject: [PATCH 6/6] enh: ITK's coordinate system deciphered --- nitransforms/io/itk.py | 55 ++++++++++++++++++++++------ nitransforms/tests/test_nonlinear.py | 13 ++++--- 2 files changed, 50 insertions(+), 18 deletions(-) diff --git a/nitransforms/io/itk.py b/nitransforms/io/itk.py index d860009e..1b68996e 100644 --- a/nitransforms/io/itk.py +++ b/nitransforms/io/itk.py @@ -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,12 +355,31 @@ 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): @@ -360,11 +387,15 @@ def to_image(cls, imgobj): 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: diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index 2029dbbf..852c905a 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/nitransforms/tests/test_nonlinear.py @@ -172,7 +172,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 @@ -201,13 +201,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" ) @@ -236,7 +237,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( ( @@ -245,6 +245,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"