Skip to content

Commit 9f5b7eb

Browse files
oestebanoesteban
committed
enh: implement ITK densefields checks vs antsApplyTransformsToPoints 
1 parent 5031627 commit 9f5b7eb

File tree

1 file changed

+114
-35
lines changed

1 file changed

+114
-35
lines changed

‎nitransforms/tests/test_io_itk.py

Lines changed: 114 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44

55
import os
66
import shutil
7+
from subprocess import check_call
78

89
import pytest
910

@@ -14,10 +15,19 @@
1415
from scipy.io import loadmat
1516
from h5py import File as H5File
1617

17-
from nitransforms.io.base import TransformIOError, TransformFileError
1818

19-
from nitransforms.io import itk
19+
from nitransforms.io.base import TransformIOError, TransformFileError
20+
from nitransforms.io.itk import (
21+
ITKLinearTransform,
22+
ITKLinearTransformArray,
23+
ITKDisplacementsField,
24+
ITKCompositeH5,
25+
)
26+
from nitransforms import nonlinear as nitnl
2027
from nitransforms.conftest import _datadir, _testdir
28+
from nitransforms.tests.utils import get_points
29+
30+
rng = np.random.default_rng()
2131

2232
LPS = np.diag([-1, -1, 1, 1])
2333
ITK_MAT = LPS.dot(np.ones((4, 4)).dot(LPS))
@@ -29,14 +39,14 @@ def test_ITKLinearTransform(tmpdir, testdata_path):
2939
matlabfile = testdata_path / "ds-005_sub-01_from-T1_to-OASIS_affine.mat"
3040
mat = loadmat(str(matlabfile))
3141
with open(str(matlabfile), "rb") as f:
32-
itkxfm = itk.ITKLinearTransform.from_fileobj(f)
42+
itkxfm = ITKLinearTransform.from_fileobj(f)
3343
assert np.allclose(
3444
itkxfm["parameters"][:3, :3].flatten(),
3545
mat["AffineTransform_float_3_3"][:-3].flatten(),
3646
)
3747
assert np.allclose(itkxfm["offset"], mat["fixed"].reshape((3,)))
3848

39-
itkxfm = itk.ITKLinearTransform.from_filename(matlabfile)
49+
itkxfm = ITKLinearTransform.from_filename(matlabfile)
4050
assert np.allclose(
4151
itkxfm["parameters"][:3, :3].flatten(),
4252
mat["AffineTransform_float_3_3"][:-3].flatten(),
@@ -46,17 +56,17 @@ def test_ITKLinearTransform(tmpdir, testdata_path):
4656
# Test to_filename(textfiles)
4757
itkxfm.to_filename("textfile.tfm")
4858
with open("textfile.tfm") as f:
49-
itkxfm2 = itk.ITKLinearTransform.from_fileobj(f)
59+
itkxfm2 = ITKLinearTransform.from_fileobj(f)
5060
assert np.allclose(itkxfm["parameters"], itkxfm2["parameters"])
5161

5262
# Test to_filename(matlab)
5363
itkxfm.to_filename("copy.mat")
5464
with open("copy.mat", "rb") as f:
55-
itkxfm3 = itk.ITKLinearTransform.from_fileobj(f)
65+
itkxfm3 = ITKLinearTransform.from_fileobj(f)
5666
assert np.all(itkxfm["parameters"] == itkxfm3["parameters"])
5767

5868
rasmat = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0])
59-
itkxfm = itk.ITKLinearTransform.from_ras(rasmat)
69+
itkxfm = ITKLinearTransform.from_ras(rasmat)
6070
assert np.allclose(itkxfm["parameters"], ITK_MAT * rasmat)
6171
assert np.allclose(itkxfm.to_ras(), rasmat)
6272

@@ -67,9 +77,9 @@ def test_ITKLinearTransformArray(tmpdir, data_path):
6777
with open(str(data_path / "itktflist.tfm")) as f:
6878
text = f.read()
6979
f.seek(0)
70-
itklist = itk.ITKLinearTransformArray.from_fileobj(f)
80+
itklist = ITKLinearTransformArray.from_fileobj(f)
7181

72-
itklistb = itk.ITKLinearTransformArray.from_filename(data_path / "itktflist.tfm")
82+
itklistb = ITKLinearTransformArray.from_filename(data_path / "itktflist.tfm")
7383
assert itklist["nxforms"] == itklistb["nxforms"]
7484
assert all(
7585
[
@@ -85,14 +95,14 @@ def test_ITKLinearTransformArray(tmpdir, data_path):
8595
assert itklist["nxforms"] == 9
8696
assert text == itklist.to_string()
8797
with pytest.raises(TransformFileError):
88-
itk.ITKLinearTransformArray.from_string("\n".join(text.splitlines()[1:]))
98+
ITKLinearTransformArray.from_string("\n".join(text.splitlines()[1:]))
8999

90100
itklist.to_filename("copy.tfm")
91101
with open("copy.tfm") as f:
92102
copytext = f.read()
93103
assert text == copytext
94104

95-
itklist = itk.ITKLinearTransformArray(
105+
itklist = ITKLinearTransformArray(
96106
xforms=[np.random.normal(size=(4, 4)) for _ in range(4)]
97107
)
98108

@@ -108,7 +118,7 @@ def test_ITKLinearTransformArray(tmpdir, data_path):
108118
f.write(xfm.to_string())
109119

110120
with open("extracted.tfm") as f:
111-
xfm2 = itk.ITKLinearTransform.from_fileobj(f)
121+
xfm2 = ITKLinearTransform.from_fileobj(f)
112122
assert np.allclose(
113123
xfm.structarr["parameters"][:3, ...], xfm2.structarr["parameters"][:3, ...]
114124
)
@@ -118,26 +128,26 @@ def test_ITKLinearTransformArray(tmpdir, data_path):
118128
itklist.to_filename("matlablist.mat")
119129

120130
with pytest.raises(TransformFileError):
121-
xfm2 = itk.ITKLinearTransformArray.from_binary({})
131+
xfm2 = ITKLinearTransformArray.from_binary({})
122132

123133
with open("filename.mat", "ab") as f:
124134
with pytest.raises(TransformFileError):
125-
xfm2 = itk.ITKLinearTransformArray.from_fileobj(f)
135+
xfm2 = ITKLinearTransformArray.from_fileobj(f)
126136

127137

128138
@pytest.mark.parametrize("size", [(20, 20, 20), (20, 20, 20, 3)])
129139
def test_itk_disp_load(size):
130140
"""Checks field sizes."""
131141
with pytest.raises(TransformFileError):
132-
itk.ITKDisplacementsField.from_image(
142+
ITKDisplacementsField.from_image(
133143
nb.Nifti1Image(np.zeros(size), np.eye(4), None)
134144
)
135145

136146

137147
def test_itk_disp_load_intent():
138148
"""Checks whether the NIfTI intent is fixed."""
139149
with pytest.warns(UserWarning):
140-
field = itk.ITKDisplacementsField.from_image(
150+
field = ITKDisplacementsField.from_image(
141151
nb.Nifti1Image(np.zeros((20, 20, 20, 1, 3)), np.eye(4), None)
142152
)
143153

@@ -161,7 +171,7 @@ def test_itk_h5(tmpdir, only_linear, h5_path, nxforms):
161171
assert (
162172
len(
163173
list(
164-
itk.ITKCompositeH5.from_filename(
174+
ITKCompositeH5.from_filename(
165175
h5_path,
166176
only_linear=only_linear,
167177
)
@@ -174,7 +184,7 @@ def test_itk_h5(tmpdir, only_linear, h5_path, nxforms):
174184

175185
with pytest.raises(TransformFileError):
176186
list(
177-
itk.ITKCompositeH5.from_filename(
187+
ITKCompositeH5.from_filename(
178188
h5_path.absolute().name.replace(".h5", ".x5"),
179189
only_linear=only_linear,
180190
)
@@ -190,30 +200,28 @@ def test_itk_h5(tmpdir, only_linear, h5_path, nxforms):
190200
xfm["TransformType"][0] = b"InventTransform"
191201

192202
with pytest.raises(TransformIOError):
193-
itk.ITKCompositeH5.from_filename("test.h5")
203+
ITKCompositeH5.from_filename("test.h5")
194204

195205

196206
def test_itk_linear_h5(tmpdir, data_path, testdata_path):
197207
"""Check different lower-level loading options."""
198208

199209
# File loadable with transform array
200-
h5xfm = itk.ITKLinearTransformArray.from_filename(
201-
data_path / "affine-antsComposite.h5"
202-
)
210+
h5xfm = ITKLinearTransformArray.from_filename(data_path / "affine-antsComposite.h5")
203211
assert len(h5xfm.xforms) == 1
204212

205213
with open(data_path / "affine-antsComposite.h5", "rb") as f:
206-
h5xfm = itk.ITKLinearTransformArray.from_fileobj(f)
214+
h5xfm = ITKLinearTransformArray.from_fileobj(f)
207215
assert len(h5xfm.xforms) == 1
208216

209217
# File loadable with single affine object
210-
itk.ITKLinearTransform.from_filename(data_path / "affine-antsComposite.h5")
218+
ITKLinearTransform.from_filename(data_path / "affine-antsComposite.h5")
211219

212220
with open(data_path / "affine-antsComposite.h5", "rb") as f:
213-
itk.ITKLinearTransform.from_fileobj(f)
221+
ITKLinearTransform.from_fileobj(f)
214222

215223
# Exercise only_linear
216-
itk.ITKCompositeH5.from_filename(
224+
ITKCompositeH5.from_filename(
217225
testdata_path
218226
/ "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5",
219227
only_linear=True,
@@ -231,16 +239,16 @@ def test_itk_linear_h5(tmpdir, data_path, testdata_path):
231239
xfm["TransformFixedParameters"] = np.zeros(3, dtype=float)
232240

233241
# File loadable with transform array
234-
h5xfm = itk.ITKLinearTransformArray.from_filename("test.h5")
242+
h5xfm = ITKLinearTransformArray.from_filename("test.h5")
235243
assert len(h5xfm.xforms) == 2
236244

237245
# File loadable with generalistic object (NOTE we directly access the list)
238-
h5xfm = itk.ITKCompositeH5.from_filename("test.h5")
246+
h5xfm = ITKCompositeH5.from_filename("test.h5")
239247
assert len(h5xfm) == 2
240248

241249
# Error raised if the we try to use the single affine loader
242250
with pytest.raises(TransformIOError):
243-
itk.ITKLinearTransform.from_filename("test.h5")
251+
ITKLinearTransform.from_filename("test.h5")
244252

245253
shutil.copy(data_path / "affine-antsComposite.h5", "test.h5")
246254
os.chmod("test.h5", 0o666)
@@ -251,12 +259,12 @@ def test_itk_linear_h5(tmpdir, data_path, testdata_path):
251259
del h5group["1"]
252260

253261
# File loadable with generalistic object
254-
h5xfm = itk.ITKCompositeH5.from_filename("test.h5")
262+
h5xfm = ITKCompositeH5.from_filename("test.h5")
255263
assert len(h5xfm) == 0
256264

257265
# Error raised if the we try to use the single affine loader
258266
with pytest.raises(TransformIOError):
259-
itk.ITKLinearTransform.from_filename("test.h5")
267+
ITKLinearTransform.from_filename("test.h5")
260268

261269

262270
# Added tests for h5 orientation bug
@@ -285,7 +293,7 @@ def test_itk_h5_field_order(tmp_path):
285293
g1["TransformFixedParameters"] = fixed
286294
g1["TransformParameters"] = params
287295

288-
img = itk.ITKCompositeH5.from_filename(fname)[0]
296+
img = ITKCompositeH5.from_filename(fname)[0]
289297
expected = np.moveaxis(field, 0, -1)
290298
expected[..., (0, 1)] *= -1
291299
assert np.allclose(img.get_fdata(), expected)
@@ -314,9 +322,9 @@ def _load_composite_testdata(data_path):
314322
def test_itk_h5_displacement_mismatch(testdata_path):
315323
"""Composite displacements should match the standalone field"""
316324
h5file, warpfile = _load_composite_testdata(testdata_path)
317-
xforms = itk.ITKCompositeH5.from_filename(h5file)
325+
xforms = ITKCompositeH5.from_filename(h5file)
318326
field_h5 = xforms[1]
319-
field_img = itk.ITKDisplacementsField.from_filename(warpfile)
327+
field_img = ITKDisplacementsField.from_filename(warpfile)
320328

321329
np.testing.assert_array_equal(
322330
np.asanyarray(field_h5.dataobj), np.asanyarray(field_img.dataobj)
@@ -368,7 +376,78 @@ def test_itk_h5_field_order_fortran(tmp_path):
368376
g1["TransformFixedParameters"] = fixed
369377
g1["TransformParameters"] = params
370378

371-
img = itk.ITKCompositeH5.from_filename(fname)[0]
379+
img = ITKCompositeH5.from_filename(fname)[0]
372380
expected = np.moveaxis(field, 0, -1)
373381
expected[..., (0, 1)] *= -1
374382
assert np.allclose(img.get_fdata(), expected)
383+
384+
385+
# Tests against ANTs' ``antsApplyTransformsToPoints``
386+
@pytest.mark.parametrize("ongrid", [True, False])
387+
def test_densefield_map_vs_ants(testdata_path, tmp_path, ongrid):
388+
"""Map points with DenseFieldTransform and compare to ANTs."""
389+
warpfile = (
390+
testdata_path
391+
/ "regressions"
392+
/ ("01_ants_t1_to_mniComposite_DisplacementFieldTransform.nii.gz")
393+
)
394+
if not warpfile.exists():
395+
pytest.skip("Composite transform test data not available")
396+
397+
nii = ITKDisplacementsField.from_filename(warpfile)
398+
399+
# Get sampling indices
400+
coords_xyz, points_ijk, grid_xyz, shape, ref_affine, reference, subsample = (
401+
get_points(nii, ongrid, npoints=5, rng=rng)
402+
)
403+
coords_map = grid_xyz.reshape(*shape, 3)
404+
405+
csvin = tmp_path / "fixed_coords.csv"
406+
csvout = tmp_path / "moving_coords.csv"
407+
408+
# antsApplyTransformsToPoints wants LPS coordinates, see last post at
409+
# http://sourceforge.net/p/advants/discussion/840261/thread/2a1e9307/
410+
lps_xyz = coords_xyz.copy() * (-1, -1, 1)
411+
np.savetxt(csvin, lps_xyz, delimiter=",", header="x,y,z", comments="")
412+
413+
cmd = f"antsApplyTransformsToPoints -d 3 -i {csvin} -o {csvout} -t {warpfile}"
414+
exe = cmd.split()[0]
415+
if not shutil.which(exe):
416+
pytest.skip(f"Command {exe} not found on host")
417+
check_call(cmd, shell=True)
418+
419+
ants_res = np.genfromtxt(csvout, delimiter=",", names=True)
420+
421+
# antsApplyTransformsToPoints writes LPS coordinates, see last post at
422+
# http://sourceforge.net/p/advants/discussion/840261/thread/2a1e9307/
423+
ants_pts = np.vstack([ants_res[n] for n in ("x", "y", "z")]).T * (-1, -1, 1)
424+
425+
xfm = nitnl.DenseFieldTransform(nii, reference=reference)
426+
mapped = xfm.map(coords_xyz)
427+
428+
if ongrid:
429+
ants_mapped_xyz = ants_pts.reshape(*shape, 3)
430+
nit_mapped_xyz = mapped.reshape(*shape, 3)
431+
432+
nb.load(warpfile).to_filename(tmp_path / "original_ants_deltas.nii.gz")
433+
434+
nb.Nifti1Image(coords_map, ref_affine, None).to_filename(
435+
tmp_path / "baseline_positions.nii.gz"
436+
)
437+
438+
nii.to_filename(tmp_path / "original_interpreted_deltas.nii.gz")
439+
440+
nb.Nifti1Image(nit_mapped_xyz, ref_affine, None).to_filename(
441+
tmp_path / "nit_deformation_xyz.nii.gz"
442+
)
443+
nb.Nifti1Image(ants_mapped_xyz - coords_map, ref_affine, None).to_filename(
444+
tmp_path / "ants_deltas_xyz.nii.gz"
445+
)
446+
nb.Nifti1Image(nit_mapped_xyz - coords_map, ref_affine, None).to_filename(
447+
tmp_path / "nit_deltas_xyz.nii.gz"
448+
)
449+
450+
# When transforming off-grid points, rounding errors are large
451+
atol = 0 if ongrid else 1e-1
452+
rtol = 1e-4 if ongrid else 1e-3
453+
np.testing.assert_allclose(mapped, ants_pts, atol=atol, rtol=rtol)

0 commit comments

Comments
 (0)