From 0d687c22410f67c0279f41f8a101cf28af2c9643 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 22 Jul 2025 17:37:17 +0200 Subject: [PATCH 01/14] enh: add failing test cases demonstrating #137 and #171 In particular, leverages @feilong's dataset he posted on #171 and adds it to the test dataset. Co-authored-by: Feilong Ma --- nitransforms/tests/test_io.py | 106 ++++++++++++++++++++++++++++++++++ 1 file changed, 106 insertions(+) diff --git a/nitransforms/tests/test_io.py b/nitransforms/tests/test_io.py index 0cc79d15..7058cdc2 100644 --- a/nitransforms/tests/test_io.py +++ b/nitransforms/tests/test_io.py @@ -672,3 +672,109 @@ def test_itk_linear_h5(tmpdir, data_path, testdata_path): # Error raised if the we try to use the single affine loader with pytest.raises(TransformIOError): itk.ITKLinearTransform.from_filename("test.h5") + +# Added tests for h5 orientation bug + + +@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) + + +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" + ) + 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): + """Composite displacements should match the standalone field""" + h5file, warpfile = _load_composite_testdata(testdata_path) + xforms = itk.ITKCompositeH5.from_filename(h5file) + field_h5 = xforms[1] + field_img = itk.ITKDisplacementsField.from_filename(warpfile) + + np.testing.assert_array_equal( + np.asanyarray(field_h5.dataobj), np.asanyarray(field_img.dataobj) + ) + + +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. + + References + ---------- + .. [1] ITK Software Guide. https://itk.org/ItkSoftwareGuide.pdf + """ + h5file, warpfile = _load_composite_testdata(testdata_path) + + 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) From df85a8754ff30f4c82fe718e2721814fdbf5e872 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 23 Jul 2025 16:55:05 +0200 Subject: [PATCH 02/14] ci: separate coverage from test results --- .circleci/config.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index f1bab746..190bf59f 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -136,21 +136,21 @@ jobs: name: Run unit tests no_output_timeout: 2h command: | - mkdir -p /tmp/tests/{artifacts,summaries} + mkdir -p /tmp/tests/{artifacts,summaries,coverage} docker run -u $( id -u ) -it --rm \ -w /src/nitransforms -v $PWD:/src/nitransforms \ -v /tmp/data/nitransforms-tests:/data -e TEST_DATA_HOME=/data \ - -e COVERAGE_FILE=/tmp/summaries/.pytest.coverage \ + -e COVERAGE_FILE=/tmp/coverage/.pytest.coverage \ -v /tmp/fslicense/license.txt:/opt/freesurfer/license.txt:ro \ -v /tmp/tests:/tmp nitransforms:latest \ pytest --junit-xml=/tmp/summaries/pytest.xml \ - --cov nitransforms --cov-report xml:/tmp/summaries/unittests.xml \ + --cov nitransforms --cov-report xml:/tmp/coverage/unittests.xml \ nitransforms/ - run: name: Submit unit test coverage command: | cd /tmp/src/nitransforms - python3 -m codecov --file /tmp/tests/summaries/unittests.xml \ + python3 -m codecov --file /tmp/tests/coverage/unittests.xml \ --flags unittests -e CIRCLE_JOB - run: name: Clean up tests directory From 314b69fb63672cb3d05314c9ecb9ace03d69bfa5 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 12 Aug 2025 20:12:31 +0200 Subject: [PATCH 03/14] fix: revise ndcoords and ndindex of spatial references --- nitransforms/base.py | 18 ++++---- nitransforms/nonlinear.py | 74 ++++++++++++++++----------------- nitransforms/resampling.py | 23 ++++++---- nitransforms/tests/test_base.py | 12 ++++-- 4 files changed, 66 insertions(+), 61 deletions(-) diff --git a/nitransforms/base.py b/nitransforms/base.py index 6e1634c6..eb6c2785 100644 --- a/nitransforms/base.py +++ b/nitransforms/base.py @@ -202,30 +202,26 @@ def inverse(self): def ndindex(self): """List the indexes corresponding to the space grid.""" if self._ndindex is None: - indexes = tuple([np.arange(s) for s in self._shape]) - self._ndindex = np.array(np.meshgrid(*indexes, indexing="ij")).reshape( - self._ndim, self._npoints - ) + indexes = np.mgrid[ + 0:self._shape[0], 0:self._shape[1], 0:self._shape[2] + ] + self._ndindex = indexes.reshape((indexes.shape[0], -1)).T return self._ndindex @property def ndcoords(self): """List the physical coordinates of this gridded space samples.""" if self._coords is None: - self._coords = np.tensordot( - self._affine, - np.vstack((self.ndindex, np.ones((1, self._npoints)))), - axes=1, - )[:3, ...] + self._coords = self.ras(self.ndindex) return self._coords def ras(self, ijk): """Get RAS+ coordinates from input indexes.""" - return _apply_affine(ijk, self._affine, self._ndim) + return _apply_affine(ijk, self._affine, self._ndim).T def index(self, x): """Get the image array's indexes corresponding to coordinates.""" - return _apply_affine(x, self._inverse, self._ndim) + return _apply_affine(x, self._inverse, self._ndim).T def _to_hdf5(self, group): group.attrs["Type"] = "image" diff --git a/nitransforms/nonlinear.py b/nitransforms/nonlinear.py index 24e043c2..fe0b18d3 100644 --- a/nitransforms/nonlinear.py +++ b/nitransforms/nonlinear.py @@ -65,50 +65,47 @@ def __init__(self, field=None, is_deltas=True, reference=None): """ + if field is None and reference is None: - raise TransformError("DenseFieldTransforms require a spatial reference") + raise TransformError("cannot initialize field") super().__init__() - self._is_deltas = is_deltas + if field is not None: + field = _ensure_image(field) + # Extract data if nibabel object otherwise assume numpy array + _data = np.squeeze( + np.asanyarray(field.dataobj) + if hasattr(field, "dataobj") + else field.copy() + ) try: self.reference = ImageGrid(reference if reference is not None else field) except AttributeError: raise TransformError( - "Field must be a spatial image if reference is not provided" + "field must be a spatial image if reference is not provided" if reference is None - else "Reference is not a spatial image" + else "reference is not a spatial image" ) fieldshape = (*self.reference.shape, self.reference.ndim) - if field is not None: - field = _ensure_image(field) - self._field = np.squeeze( - np.asanyarray(field.dataobj) if hasattr(field, "dataobj") else field - ) - if fieldshape != self._field.shape: - raise TransformError( - f"Shape of the field ({'x'.join(str(i) for i in self._field.shape)}) " - f"doesn't match that of the reference({'x'.join(str(i) for i in fieldshape)})" - ) - else: - self._field = np.zeros(fieldshape, dtype="float32") - self._is_deltas = True - - if self._field.shape[-1] != self.ndim: + if field is None: + _data = np.zeros(fieldshape) + elif fieldshape != _data.shape: raise TransformError( - "The number of components of the field (%d) does not match " - "the number of dimensions (%d)" % (self._field.shape[-1], self.ndim) + f"Shape of the field ({'x'.join(str(i) for i in _data.shape)}) " + f"doesn't match that of the reference({'x'.join(str(i) for i in fieldshape)})" ) + self._is_deltas = is_deltas + self._field = self.reference.ndcoords.reshape(fieldshape) + if self.is_deltas: - self._deltas = ( - self._field.copy() - ) # IMPORTANT: you don't want to update deltas - # Convert from displacements (deltas) to deformations fields - # (just add its origin to each delta vector) - self._field += self.reference.ndcoords.T.reshape(fieldshape) + self._deltas = _data.copy() + self._field += self._deltas + else: + self._field = _data.copy() def __repr__(self): """Beautify the python representation.""" @@ -153,7 +150,7 @@ def map(self, x, inverse=False): ... test_dir / "someones_displacement_field.nii.gz", ... is_deltas=False, ... ) - >>> xfm.map([-6.5, -36., -19.5]).tolist() + >>> xfm.map([[-6.5, -36., -19.5]]).tolist() [[0.0, -0.47516798973083496, 0.0]] >>> xfm.map([[-6.5, -36., -19.5], [-1., -41.5, -11.25]]).tolist() @@ -170,8 +167,8 @@ def map(self, x, inverse=False): ... test_dir / "someones_displacement_field.nii.gz", ... is_deltas=True, ... ) - >>> xfm.map([[-6.5, -36., -19.5], [-1., -41.5, -11.25]]).tolist() - [[-6.5, -36.47516632080078, -19.5], [-1.0, -42.03835678100586, -11.25]] + >>> xfm.map([[-6.5, -36., -19.5], [-1., -41.5, -11.25]]).tolist() # doctest: +ELLIPSIS + [[-6.5, -36.475..., -19.5], [-1.0, -42.038..., -11.25]] >>> np.array_str( ... xfm.map([[-6.7, -36.3, -19.2], [-1., -41.5, -11.25]]), @@ -185,18 +182,19 @@ def map(self, x, inverse=False): if inverse is True: raise NotImplementedError - ijk = self.reference.index(x) + ijk = self.reference.index(np.array(x, dtype="float32")) indexes = np.round(ijk).astype("int") + ongrid = np.where(np.linalg.norm(ijk - indexes, axis=1) < 1e-3)[0] - if np.all(np.abs(ijk - indexes) < 1e-3): - indexes = tuple(tuple(i) for i in indexes) - return self._field[indexes] + if ongrid.size == np.shape(x)[0]: + # return self._field[*indexes.T, :] # From Python 3.11 + return self._field[tuple(indexes.T) + (np.s_[:],)] - new_map = np.vstack( + mapped_coords = np.vstack( tuple( map_coordinates( self._field[..., i], - ijk, + ijk.T, order=3, mode="constant", cval=np.nan, @@ -207,8 +205,8 @@ def map(self, x, inverse=False): ).T # Set NaN values back to the original coordinates value = no displacement - new_map[np.isnan(new_map)] = np.array(x)[np.isnan(new_map)] - return new_map + mapped_coords[np.isnan(mapped_coords)] = np.array(x)[np.isnan(mapped_coords)] + return mapped_coords def __matmul__(self, b): """ diff --git a/nitransforms/resampling.py b/nitransforms/resampling.py index 98ef4454..6ade3eff 100644 --- a/nitransforms/resampling.py +++ b/nitransforms/resampling.py @@ -253,7 +253,7 @@ def apply( serialize_4d = n_resamplings >= serialize_nvols targets = None - ref_ndcoords = _ref.ndcoords.T + ref_ndcoords = _ref.ndcoords if hasattr(transform, "to_field") and callable(transform.to_field): targets = ImageGrid(spatialimage).index( _as_homogeneous( @@ -271,11 +271,8 @@ def apply( else targets ) - if targets.ndim == 3: - targets = np.rollaxis(targets, targets.ndim - 1, 0) - else: - assert targets.ndim == 2 - targets = targets[np.newaxis, ...] + if targets.ndim == 2: + targets = targets.T[np.newaxis, ...] if serialize_4d: data = ( @@ -290,6 +287,9 @@ def apply( (len(ref_ndcoords), n_resamplings), dtype=input_dtype, order="F" ) + if targets.ndim == 3: + targets = np.rollaxis(targets, targets.ndim - 1, 1) + resampled = asyncio.run( _apply_serial( data, @@ -311,6 +311,9 @@ def apply( else: data = np.asanyarray(spatialimage.dataobj, dtype=input_dtype) + if targets.ndim == 3: + targets = np.rollaxis(targets, targets.ndim - 1, 0) + if data_nvols == 1 and xfm_nvols == 1: targets = np.squeeze(targets) assert targets.ndim == 2 @@ -320,15 +323,19 @@ def apply( if xfm_nvols > 1: assert targets.ndim == 3 - n_time, n_dim, n_vox = targets.shape + + # Targets must have shape (n_dim x n_time x n_vox) + n_dim, n_time, n_vox = targets.shape # Reshape to (3, n_time x n_vox) - ijk_targets = np.rollaxis(targets, 0, 2).reshape((n_dim, -1)) + ijk_targets = targets.reshape((n_dim, -1)) time_row = np.repeat(np.arange(n_time), n_vox)[None, :] # Now targets is (4, n_vox x n_time), with indexes (t, i, j, k) # t is the slowest-changing axis, so we put it first targets = np.vstack((time_row, ijk_targets)) data = np.rollaxis(data, data.ndim - 1, 0) + else: + targets = targets.T resampled = ndi.map_coordinates( data, diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index 45611745..fe9c8d20 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -55,20 +55,24 @@ def test_ImageGrid(get_testdata, image_orientation): assert np.allclose(np.squeeze(img.ras(ijk[0])), xyz[0]) assert np.allclose(np.round(img.index(xyz[0])), ijk[0]) - assert np.allclose(img.ras(ijk).T, xyz) - assert np.allclose(np.round(img.index(xyz)).T, ijk) + assert np.allclose(img.ras(ijk), xyz) + assert np.allclose(np.round(img.index(xyz)), ijk) # nd index / coords idxs = img.ndindex coords = img.ndcoords assert len(idxs.shape) == len(coords.shape) == 2 - assert idxs.shape[0] == coords.shape[0] == img.ndim == 3 - assert idxs.shape[1] == coords.shape[1] == img.npoints == np.prod(im.shape) + assert idxs.shape[1] == coords.shape[1] == img.ndim == 3 + assert idxs.shape[0] == coords.shape[0] == img.npoints == np.prod(im.shape) 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.""" From 12777390fc033ed201bcaf69092c4b6d8592234b Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 12 Aug 2025 20:33:38 +0200 Subject: [PATCH 04/14] enh: update tests --- nitransforms/tests/test_io_itk.py | 41 ++++++++++ nitransforms/tests/test_nonlinear.py | 112 ++++++-------------------- nitransforms/tests/test_resampling.py | 8 ++ nitransforms/tests/test_x5.py | 50 +++++++++++- 4 files changed, 121 insertions(+), 90 deletions(-) create mode 100644 nitransforms/tests/test_io_itk.py diff --git a/nitransforms/tests/test_io_itk.py b/nitransforms/tests/test_io_itk.py new file mode 100644 index 00000000..f952531d --- /dev/null +++ b/nitransforms/tests/test_io_itk.py @@ -0,0 +1,41 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Test io module for ITK transforms.""" + +import pytest + +import numpy as np +import nibabel as nb + +from nitransforms.base import TransformError +from nitransforms.io.base import TransformFileError +from nitransforms.io.itk import ITKDisplacementsField +from nitransforms.nonlinear import ( + DenseFieldTransform, +) + + +@pytest.mark.parametrize("size", [(20, 20, 20), (20, 20, 20, 3)]) +def test_itk_disp_load(size): + """Checks field sizes.""" + with pytest.raises(TransformFileError): + ITKDisplacementsField.from_image( + nb.Nifti1Image(np.zeros(size), np.eye(4), None) + ) + + +@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.""" + with pytest.raises(TransformError): + DenseFieldTransform(nb.Nifti1Image(np.zeros(size), np.eye(4), None)) + + +def test_itk_disp_load_intent(): + """Checks whether the NIfTI intent is fixed.""" + with pytest.warns(UserWarning): + field = ITKDisplacementsField.from_image( + nb.Nifti1Image(np.zeros((20, 20, 20, 1, 3)), np.eye(4), None) + ) + + assert field.header.get_intent()[0] == "vector" diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index 936a62f6..76c1acaa 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/nitransforms/tests/test_nonlinear.py @@ -9,39 +9,10 @@ import nibabel as nb from nitransforms.resampling import apply from nitransforms.base import TransformError -from nitransforms.io.base import TransformFileError from nitransforms.nonlinear import ( BSplineFieldTransform, DenseFieldTransform, ) -from nitransforms import io -from ..io.itk import ITKDisplacementsField - - -@pytest.mark.parametrize("size", [(20, 20, 20), (20, 20, 20, 3)]) -def test_itk_disp_load(size): - """Checks field sizes.""" - with pytest.raises(TransformFileError): - ITKDisplacementsField.from_image( - nb.Nifti1Image(np.zeros(size), np.eye(4), None) - ) - - -@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.""" - with pytest.raises(TransformError): - DenseFieldTransform(nb.Nifti1Image(np.zeros(size), np.eye(4), None)) - - -def test_itk_disp_load_intent(): - """Checks whether the NIfTI intent is fixed.""" - with pytest.warns(UserWarning): - field = ITKDisplacementsField.from_image( - nb.Nifti1Image(np.zeros((20, 20, 20, 1, 3)), np.eye(4), None) - ) - - assert field.header.get_intent()[0] == "vector" def test_displacements_init(): @@ -96,7 +67,27 @@ def test_bsplines_references(testdata_path): ) +@pytest.mark.xfail( + reason="Disable while #266 is developed.", + strict=False, +) def test_bspline(tmp_path, testdata_path): + """ + Cross-check B-Splines and deformation field. + + This test is disabled and will be split into two separate tests. + The current implementation will be moved into test_resampling.py, + since that's what it actually tests. + + In GH-266, this test will be re-implemented by testing the equivalence + of the B-Spline and deformation field transforms by calling the + transform's `map()` method on points. + + """ + assert True + + +def test_map_bspline_vs_displacement(tmp_path, testdata_path): """Cross-check B-Splines and deformation field.""" os.chdir(str(tmp_path)) @@ -104,68 +95,11 @@ def test_bspline(tmp_path, testdata_path): disp_name = testdata_path / "someones_displacement_field.nii.gz" bs_name = testdata_path / "someones_bspline_coefficients.nii.gz" - bsplxfm = BSplineFieldTransform(bs_name, reference=img_name) + bsplxfm = BSplineFieldTransform(bs_name, reference=img_name).to_field() dispxfm = 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 - ) - - -@pytest.mark.parametrize("is_deltas", [True, False]) -def test_densefield_x5_roundtrip(tmp_path, is_deltas): - """Ensure dense field transforms roundtrip via X5.""" - ref = nb.Nifti1Image(np.zeros((2, 2, 2), dtype="uint8"), np.eye(4)) - disp = nb.Nifti1Image(np.random.rand(2, 2, 2, 3).astype("float32"), np.eye(4)) - - xfm = DenseFieldTransform(disp, is_deltas=is_deltas, reference=ref) - - node = xfm.to_x5(metadata={"GeneratedBy": "pytest"}) - assert node.type == "nonlinear" - assert node.subtype == "densefield" - assert node.representation == "displacements" if is_deltas else "deformations" - assert node.domain.size == ref.shape - assert node.metadata["GeneratedBy"] == "pytest" - - fname = tmp_path / "test.x5" - io.x5.to_filename(fname, [node]) - - xfm2 = DenseFieldTransform.from_filename(fname, fmt="X5") - - assert xfm2.reference.shape == ref.shape - assert np.allclose(xfm2.reference.affine, ref.affine) - assert xfm == xfm2 - - -def test_bspline_to_x5(tmp_path): - """Check BSpline transforms export to X5.""" - coeff = nb.Nifti1Image(np.zeros((2, 2, 2, 3), dtype="float32"), np.eye(4)) - ref = nb.Nifti1Image(np.zeros((2, 2, 2), dtype="uint8"), np.eye(4)) - - xfm = BSplineFieldTransform(coeff, reference=ref) - node = xfm.to_x5(metadata={"tool": "pytest"}) - assert node.type == "nonlinear" - assert node.subtype == "bspline" - assert node.representation == "coefficients" - assert node.metadata["tool"] == "pytest" - - fname = tmp_path / "bspline.x5" - io.x5.to_filename(fname, [node]) - - xfm2 = BSplineFieldTransform.from_filename(fname, fmt="X5") - assert np.allclose(xfm._coeffs, xfm2._coeffs) - assert xfm2.reference.shape == ref.shape - assert np.allclose(xfm2.reference.affine, ref.affine) + # Interpolating the field should be reasonably similar + np.testing.assert_allclose(dispxfm._field, bsplxfm._field, atol=1e-1, rtol=1e-4) @pytest.mark.parametrize("is_deltas", [True, False]) diff --git a/nitransforms/tests/test_resampling.py b/nitransforms/tests/test_resampling.py index b65bf579..ca56a7e8 100644 --- a/nitransforms/tests/test_resampling.py +++ b/nitransforms/tests/test_resampling.py @@ -149,6 +149,10 @@ 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)]) @@ -236,6 +240,10 @@ 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.""" diff --git a/nitransforms/tests/test_x5.py b/nitransforms/tests/test_x5.py index 89b49e06..a213c7cd 100644 --- a/nitransforms/tests/test_x5.py +++ b/nitransforms/tests/test_x5.py @@ -1,8 +1,10 @@ import numpy as np +import nibabel as nb import pytest from h5py import File as H5File -from ..io.x5 import X5Transform, X5Domain, to_filename, from_filename +from nitransforms.nonlinear import DenseFieldTransform, BSplineFieldTransform +from nitransforms.io.x5 import X5Transform, X5Domain, to_filename, from_filename def test_x5_transform_defaults(): @@ -75,3 +77,49 @@ def test_from_filename_invalid(tmp_path): with pytest.raises(TypeError): from_filename(fname) + + +@pytest.mark.parametrize("is_deltas", [True, False]) +def test_densefield_x5_roundtrip(tmp_path, is_deltas): + """Ensure dense field transforms roundtrip via X5.""" + ref = nb.Nifti1Image(np.zeros((2, 2, 2), dtype="uint8"), np.eye(4)) + disp = nb.Nifti1Image(np.random.rand(2, 2, 2, 3).astype("float32"), np.eye(4)) + + xfm = DenseFieldTransform(disp, is_deltas=is_deltas, reference=ref) + + node = xfm.to_x5(metadata={"GeneratedBy": "pytest"}) + assert node.type == "nonlinear" + assert node.subtype == "densefield" + assert node.representation == "displacements" if is_deltas else "deformations" + assert node.domain.size == ref.shape + assert node.metadata["GeneratedBy"] == "pytest" + + fname = tmp_path / "test.x5" + to_filename(fname, [node]) + + xfm2 = DenseFieldTransform.from_filename(fname, fmt="X5") + + assert xfm2.reference.shape == ref.shape + assert np.allclose(xfm2.reference.affine, ref.affine) + assert xfm == xfm2 + + +def test_bspline_to_x5(tmp_path): + """Check BSpline transforms export to X5.""" + coeff = nb.Nifti1Image(np.zeros((2, 2, 2, 3), dtype="float32"), np.eye(4)) + ref = nb.Nifti1Image(np.zeros((2, 2, 2), dtype="uint8"), np.eye(4)) + + xfm = BSplineFieldTransform(coeff, reference=ref) + node = xfm.to_x5(metadata={"tool": "pytest"}) + assert node.type == "nonlinear" + assert node.subtype == "bspline" + assert node.representation == "coefficients" + assert node.metadata["tool"] == "pytest" + + fname = tmp_path / "bspline.x5" + to_filename(fname, [node]) + + xfm2 = BSplineFieldTransform.from_filename(fname, fmt="X5") + assert np.allclose(xfm._coeffs, xfm2._coeffs) + assert xfm2.reference.shape == ref.shape + assert np.allclose(xfm2.reference.affine, ref.affine) From 6f5a113b7dfa35e9df3572e0557fc489aabff4f2 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 13 Aug 2025 00:38:48 +0200 Subject: [PATCH 05/14] rf: move tests to better locations --- nitransforms/tests/test_io.py | 389 +++------------------------ nitransforms/tests/test_io_itk.py | 363 +++++++++++++++++++++++-- nitransforms/tests/test_nonlinear.py | 7 + 3 files changed, 394 insertions(+), 365 deletions(-) diff --git a/nitransforms/tests/test_io.py b/nitransforms/tests/test_io.py index 7058cdc2..36216d75 100644 --- a/nitransforms/tests/test_io.py +++ b/nitransforms/tests/test_io.py @@ -1,19 +1,17 @@ # 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 import filecmp import shutil import numpy as np import pytest -from h5py import File as H5File 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.io import ( afni, @@ -27,14 +25,9 @@ FSLinearTransformArray as LTA, ) from nitransforms.io.base import LinearParameters, TransformIOError, TransformFileError -from nitransforms.conftest import _datadir, _testdir from nitransforms.resampling import apply -LPS = np.diag([-1, -1, 1, 1]) -ITK_MAT = LPS.dot(np.ones((4, 4)).dot(LPS)) - - def test_VolumeGeometry(tmpdir, get_testdata): vg = VG() assert vg["valid"] == 0 @@ -68,11 +61,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] @@ -290,108 +285,6 @@ def test_LinearList_common(tmpdir, data_path, sw, image_orientation, get_testdat ) -def test_ITKLinearTransform(tmpdir, testdata_path): - tmpdir.chdir() - - matlabfile = testdata_path / "ds-005_sub-01_from-T1_to-OASIS_affine.mat" - mat = loadmat(str(matlabfile)) - with open(str(matlabfile), "rb") as f: - itkxfm = itk.ITKLinearTransform.from_fileobj(f) - assert np.allclose( - itkxfm["parameters"][:3, :3].flatten(), - mat["AffineTransform_float_3_3"][:-3].flatten(), - ) - assert np.allclose(itkxfm["offset"], mat["fixed"].reshape((3,))) - - itkxfm = itk.ITKLinearTransform.from_filename(matlabfile) - assert np.allclose( - itkxfm["parameters"][:3, :3].flatten(), - mat["AffineTransform_float_3_3"][:-3].flatten(), - ) - assert np.allclose(itkxfm["offset"], mat["fixed"].reshape((3,))) - - # Test to_filename(textfiles) - itkxfm.to_filename("textfile.tfm") - with open("textfile.tfm") as f: - itkxfm2 = itk.ITKLinearTransform.from_fileobj(f) - assert np.allclose(itkxfm["parameters"], itkxfm2["parameters"]) - - # Test to_filename(matlab) - itkxfm.to_filename("copy.mat") - with open("copy.mat", "rb") as f: - itkxfm3 = itk.ITKLinearTransform.from_fileobj(f) - assert np.all(itkxfm["parameters"] == itkxfm3["parameters"]) - - rasmat = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0]) - itkxfm = itk.ITKLinearTransform.from_ras(rasmat) - assert np.allclose(itkxfm["parameters"], ITK_MAT * rasmat) - assert np.allclose(itkxfm.to_ras(), rasmat) - - -def test_ITKLinearTransformArray(tmpdir, data_path): - tmpdir.chdir() - - with open(str(data_path / "itktflist.tfm")) as f: - text = f.read() - f.seek(0) - itklist = itk.ITKLinearTransformArray.from_fileobj(f) - - itklistb = itk.ITKLinearTransformArray.from_filename(data_path / "itktflist.tfm") - assert itklist["nxforms"] == itklistb["nxforms"] - assert all( - [ - np.allclose(x1["parameters"], x2["parameters"]) - for x1, x2 in zip(itklist.xforms, itklistb.xforms) - ] - ) - - tmpdir.join("empty.mat").write("") - with pytest.raises(TransformFileError): - itklistb.from_filename("empty.mat") - - assert itklist["nxforms"] == 9 - assert text == itklist.to_string() - with pytest.raises(TransformFileError): - itk.ITKLinearTransformArray.from_string("\n".join(text.splitlines()[1:])) - - itklist.to_filename("copy.tfm") - with open("copy.tfm") as f: - copytext = f.read() - assert text == copytext - - itklist = itk.ITKLinearTransformArray( - xforms=[np.random.normal(size=(4, 4)) for _ in range(4)] - ) - - assert itklist["nxforms"] == 4 - assert itklist["xforms"][1].structarr["index"] == 1 - - with pytest.raises(KeyError): - itklist["invalid_key"] - - xfm = itklist["xforms"][1] - xfm["index"] = 1 - with open("extracted.tfm", "w") as f: - f.write(xfm.to_string()) - - with open("extracted.tfm") as f: - xfm2 = itk.ITKLinearTransform.from_fileobj(f) - assert np.allclose( - xfm.structarr["parameters"][:3, ...], xfm2.structarr["parameters"][:3, ...] - ) - - # ITK does not handle transforms lists in Matlab format - with pytest.raises(TransformFileError): - itklist.to_filename("matlablist.mat") - - with pytest.raises(TransformFileError): - xfm2 = itk.ITKLinearTransformArray.from_binary({}) - - with open("filename.mat", "ab") as f: - with pytest.raises(TransformFileError): - xfm2 = itk.ITKLinearTransformArray.from_fileobj(f) - - def test_LinearParameters(tmpdir): """Just pushes coverage up.""" tmpdir.join("file.txt").write("") @@ -418,46 +311,6 @@ def test_afni_Displacements(): afni.AFNIDisplacementsField.from_image(field) -@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), -]) -def test_itk_h5(tmpdir, only_linear, h5_path, nxforms): - """Test displacements fields.""" - assert ( - len( - list( - itk.ITKCompositeH5.from_filename( - h5_path, - only_linear=only_linear, - ) - ) - ) - == nxforms if not only_linear else 1 - ) - - with pytest.raises(TransformFileError): - list( - itk.ITKCompositeH5.from_filename( - h5_path.absolute().name.replace(".h5", ".x5"), - only_linear=only_linear, - ) - ) - - tmpdir.chdir() - shutil.copy(h5_path, "test.h5") - os.chmod("test.h5", 0o666) - - with H5File("test.h5", "r+") as h5file: - h5group = h5file["TransformGroup"] - xfm = h5group[list(h5group.keys())[1]] - xfm["TransformType"][0] = b"InventTransform" - - with pytest.raises(TransformIOError): - itk.ITKCompositeH5.from_filename("test.h5") - - @pytest.mark.parametrize( "file_type, test_file", [(LTA, "from-fsnative_to-scanner_mode-image.lta")] ) @@ -465,24 +318,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 +369,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 +383,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 +402,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 +422,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 +453,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) @@ -605,176 +467,3 @@ def _generate_reoriented(path, directions, swapaxes, parameters): hdr.set_qform(newaff, code=1) hdr.set_sform(newaff, code=1) return img.__class__(data, newaff, hdr), R - - -def test_itk_linear_h5(tmpdir, data_path, testdata_path): - """Check different lower-level loading options.""" - - # File loadable with transform array - h5xfm = itk.ITKLinearTransformArray.from_filename( - data_path / "affine-antsComposite.h5" - ) - assert len(h5xfm.xforms) == 1 - - with open(data_path / "affine-antsComposite.h5", "rb") as f: - h5xfm = itk.ITKLinearTransformArray.from_fileobj(f) - assert len(h5xfm.xforms) == 1 - - # File loadable with single affine object - 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", - only_linear=True, - ) - - tmpdir.chdir() - shutil.copy(data_path / "affine-antsComposite.h5", "test.h5") - os.chmod("test.h5", 0o666) - - with H5File("test.h5", "r+") as h5file: - h5group = h5file["TransformGroup"] - xfm = h5group.create_group("2") - xfm["TransformType"] = (b"AffineTransform", b"") - xfm["TransformParameters"] = np.zeros(12, dtype=float) - xfm["TransformFixedParameters"] = np.zeros(3, dtype=float) - - # File loadable with transform array - h5xfm = itk.ITKLinearTransformArray.from_filename("test.h5") - assert len(h5xfm.xforms) == 2 - - # File loadable with generalistic object (NOTE we directly access the list) - h5xfm = itk.ITKCompositeH5.from_filename("test.h5") - assert len(h5xfm) == 2 - - # Error raised if the we try to use the single affine loader - with pytest.raises(TransformIOError): - itk.ITKLinearTransform.from_filename("test.h5") - - shutil.copy(data_path / "affine-antsComposite.h5", "test.h5") - os.chmod("test.h5", 0o666) - - # Generate an empty h5 file - with H5File("test.h5", "r+") as h5file: - h5group = h5file["TransformGroup"] - del h5group["1"] - - # File loadable with generalistic object - h5xfm = itk.ITKCompositeH5.from_filename("test.h5") - assert len(h5xfm) == 0 - - # Error raised if the we try to use the single affine loader - with pytest.raises(TransformIOError): - itk.ITKLinearTransform.from_filename("test.h5") - -# Added tests for h5 orientation bug - - -@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) - - -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" - ) - 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): - """Composite displacements should match the standalone field""" - h5file, warpfile = _load_composite_testdata(testdata_path) - xforms = itk.ITKCompositeH5.from_filename(h5file) - field_h5 = xforms[1] - field_img = itk.ITKDisplacementsField.from_filename(warpfile) - - np.testing.assert_array_equal( - np.asanyarray(field_h5.dataobj), np.asanyarray(field_img.dataobj) - ) - - -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. - - References - ---------- - .. [1] ITK Software Guide. https://itk.org/ItkSoftwareGuide.pdf - """ - h5file, warpfile = _load_composite_testdata(testdata_path) - - 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) diff --git a/nitransforms/tests/test_io_itk.py b/nitransforms/tests/test_io_itk.py index f952531d..32423733 100644 --- a/nitransforms/tests/test_io_itk.py +++ b/nitransforms/tests/test_io_itk.py @@ -2,40 +2,373 @@ # vi: set ft=python sts=4 ts=4 sw=4 et: """Test io module for ITK transforms.""" +import os +import shutil + import pytest import numpy as np import nibabel as nb +from nibabel.eulerangles import euler2mat +from nibabel.affines import from_matvec +from scipy.io import loadmat +from h5py import File as H5File -from nitransforms.base import TransformError -from nitransforms.io.base import TransformFileError -from nitransforms.io.itk import ITKDisplacementsField -from nitransforms.nonlinear import ( - DenseFieldTransform, -) +from nitransforms.io.base import TransformIOError, TransformFileError + +from nitransforms.io import itk +from nitransforms.conftest import _datadir, _testdir + +LPS = np.diag([-1, -1, 1, 1]) +ITK_MAT = LPS.dot(np.ones((4, 4)).dot(LPS)) + + +def test_ITKLinearTransform(tmpdir, testdata_path): + tmpdir.chdir() + + matlabfile = testdata_path / "ds-005_sub-01_from-T1_to-OASIS_affine.mat" + mat = loadmat(str(matlabfile)) + with open(str(matlabfile), "rb") as f: + itkxfm = itk.ITKLinearTransform.from_fileobj(f) + assert np.allclose( + itkxfm["parameters"][:3, :3].flatten(), + mat["AffineTransform_float_3_3"][:-3].flatten(), + ) + assert np.allclose(itkxfm["offset"], mat["fixed"].reshape((3,))) + + itkxfm = itk.ITKLinearTransform.from_filename(matlabfile) + assert np.allclose( + itkxfm["parameters"][:3, :3].flatten(), + mat["AffineTransform_float_3_3"][:-3].flatten(), + ) + assert np.allclose(itkxfm["offset"], mat["fixed"].reshape((3,))) + + # Test to_filename(textfiles) + itkxfm.to_filename("textfile.tfm") + with open("textfile.tfm") as f: + itkxfm2 = itk.ITKLinearTransform.from_fileobj(f) + assert np.allclose(itkxfm["parameters"], itkxfm2["parameters"]) + + # Test to_filename(matlab) + itkxfm.to_filename("copy.mat") + with open("copy.mat", "rb") as f: + itkxfm3 = itk.ITKLinearTransform.from_fileobj(f) + assert np.all(itkxfm["parameters"] == itkxfm3["parameters"]) + + rasmat = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0]) + itkxfm = itk.ITKLinearTransform.from_ras(rasmat) + assert np.allclose(itkxfm["parameters"], ITK_MAT * rasmat) + assert np.allclose(itkxfm.to_ras(), rasmat) + + +def test_ITKLinearTransformArray(tmpdir, data_path): + tmpdir.chdir() + + with open(str(data_path / "itktflist.tfm")) as f: + text = f.read() + f.seek(0) + itklist = itk.ITKLinearTransformArray.from_fileobj(f) + + itklistb = itk.ITKLinearTransformArray.from_filename(data_path / "itktflist.tfm") + assert itklist["nxforms"] == itklistb["nxforms"] + assert all( + [ + np.allclose(x1["parameters"], x2["parameters"]) + for x1, x2 in zip(itklist.xforms, itklistb.xforms) + ] + ) + + tmpdir.join("empty.mat").write("") + with pytest.raises(TransformFileError): + itklistb.from_filename("empty.mat") + + assert itklist["nxforms"] == 9 + assert text == itklist.to_string() + with pytest.raises(TransformFileError): + itk.ITKLinearTransformArray.from_string("\n".join(text.splitlines()[1:])) + + itklist.to_filename("copy.tfm") + with open("copy.tfm") as f: + copytext = f.read() + assert text == copytext + + itklist = itk.ITKLinearTransformArray( + xforms=[np.random.normal(size=(4, 4)) for _ in range(4)] + ) + + assert itklist["nxforms"] == 4 + assert itklist["xforms"][1].structarr["index"] == 1 + + with pytest.raises(KeyError): + itklist["invalid_key"] + + xfm = itklist["xforms"][1] + xfm["index"] = 1 + with open("extracted.tfm", "w") as f: + f.write(xfm.to_string()) + + with open("extracted.tfm") as f: + xfm2 = itk.ITKLinearTransform.from_fileobj(f) + assert np.allclose( + xfm.structarr["parameters"][:3, ...], xfm2.structarr["parameters"][:3, ...] + ) + + # ITK does not handle transforms lists in Matlab format + with pytest.raises(TransformFileError): + itklist.to_filename("matlablist.mat") + + with pytest.raises(TransformFileError): + xfm2 = itk.ITKLinearTransformArray.from_binary({}) + + with open("filename.mat", "ab") as f: + with pytest.raises(TransformFileError): + xfm2 = itk.ITKLinearTransformArray.from_fileobj(f) @pytest.mark.parametrize("size", [(20, 20, 20), (20, 20, 20, 3)]) def test_itk_disp_load(size): """Checks field sizes.""" with pytest.raises(TransformFileError): - ITKDisplacementsField.from_image( + itk.ITKDisplacementsField.from_image( nb.Nifti1Image(np.zeros(size), np.eye(4), None) ) -@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.""" - with pytest.raises(TransformError): - DenseFieldTransform(nb.Nifti1Image(np.zeros(size), np.eye(4), None)) - - def test_itk_disp_load_intent(): """Checks whether the NIfTI intent is fixed.""" with pytest.warns(UserWarning): - field = ITKDisplacementsField.from_image( + field = itk.ITKDisplacementsField.from_image( nb.Nifti1Image(np.zeros((20, 20, 20, 1, 3)), np.eye(4), None) ) assert field.header.get_intent()[0] == "vector" + + +@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, + ), + ], +) +def test_itk_h5(tmpdir, only_linear, h5_path, nxforms): + """Test displacements fields.""" + assert ( + len( + list( + itk.ITKCompositeH5.from_filename( + h5_path, + only_linear=only_linear, + ) + ) + ) + == nxforms + if not only_linear + else 1 + ) + + with pytest.raises(TransformFileError): + list( + itk.ITKCompositeH5.from_filename( + h5_path.absolute().name.replace(".h5", ".x5"), + only_linear=only_linear, + ) + ) + + tmpdir.chdir() + shutil.copy(h5_path, "test.h5") + os.chmod("test.h5", 0o666) + + with H5File("test.h5", "r+") as h5file: + h5group = h5file["TransformGroup"] + xfm = h5group[list(h5group.keys())[1]] + xfm["TransformType"][0] = b"InventTransform" + + with pytest.raises(TransformIOError): + itk.ITKCompositeH5.from_filename("test.h5") + + +def test_itk_linear_h5(tmpdir, data_path, testdata_path): + """Check different lower-level loading options.""" + + # File loadable with transform array + h5xfm = itk.ITKLinearTransformArray.from_filename( + data_path / "affine-antsComposite.h5" + ) + assert len(h5xfm.xforms) == 1 + + with open(data_path / "affine-antsComposite.h5", "rb") as f: + h5xfm = itk.ITKLinearTransformArray.from_fileobj(f) + assert len(h5xfm.xforms) == 1 + + # File loadable with single affine object + 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", + only_linear=True, + ) + + tmpdir.chdir() + shutil.copy(data_path / "affine-antsComposite.h5", "test.h5") + os.chmod("test.h5", 0o666) + + with H5File("test.h5", "r+") as h5file: + h5group = h5file["TransformGroup"] + xfm = h5group.create_group("2") + xfm["TransformType"] = (b"AffineTransform", b"") + xfm["TransformParameters"] = np.zeros(12, dtype=float) + xfm["TransformFixedParameters"] = np.zeros(3, dtype=float) + + # File loadable with transform array + h5xfm = itk.ITKLinearTransformArray.from_filename("test.h5") + assert len(h5xfm.xforms) == 2 + + # File loadable with generalistic object (NOTE we directly access the list) + h5xfm = itk.ITKCompositeH5.from_filename("test.h5") + assert len(h5xfm) == 2 + + # Error raised if the we try to use the single affine loader + with pytest.raises(TransformIOError): + itk.ITKLinearTransform.from_filename("test.h5") + + shutil.copy(data_path / "affine-antsComposite.h5", "test.h5") + os.chmod("test.h5", 0o666) + + # Generate an empty h5 file + with H5File("test.h5", "r+") as h5file: + h5group = h5file["TransformGroup"] + del h5group["1"] + + # File loadable with generalistic object + h5xfm = itk.ITKCompositeH5.from_filename("test.h5") + assert len(h5xfm) == 0 + + # Error raised if the we try to use the single affine loader + with pytest.raises(TransformIOError): + itk.ITKLinearTransform.from_filename("test.h5") + + +# Added tests for h5 orientation bug +@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) + + +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") + ) + 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): + """Composite displacements should match the standalone field""" + h5file, warpfile = _load_composite_testdata(testdata_path) + xforms = itk.ITKCompositeH5.from_filename(h5file) + field_h5 = xforms[1] + field_img = itk.ITKDisplacementsField.from_filename(warpfile) + + np.testing.assert_array_equal( + np.asanyarray(field_h5.dataobj), np.asanyarray(field_img.dataobj) + ) + + +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. + + References + ---------- + .. [1] ITK Software Guide. https://itk.org/ItkSoftwareGuide.pdf + """ + h5file, warpfile = _load_composite_testdata(testdata_path) + + 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) diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index 76c1acaa..8510d993 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/nitransforms/tests/test_nonlinear.py @@ -38,6 +38,13 @@ def test_displacements_init(): ) +@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.""" + with pytest.raises(TransformError): + DenseFieldTransform(nb.Nifti1Image(np.zeros(size), np.eye(4), None)) + + def test_bsplines_init(): with pytest.raises(TransformError): BSplineFieldTransform( From 3fe53f9820c3ba0f0780462a92f798600687f6cb Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 13 Aug 2025 09:27:31 +0200 Subject: [PATCH 06/14] sty: run ruff on all --- nitransforms/__init__.py | 1 + nitransforms/cli.py | 4 +- nitransforms/conftest.py | 5 +- nitransforms/io/afni.py | 18 ++- nitransforms/io/base.py | 1 + nitransforms/io/fsl.py | 16 ++- nitransforms/io/itk.py | 14 +- nitransforms/surface.py | 186 ++++++++++++++----------- nitransforms/tests/test_cli.py | 5 +- nitransforms/tests/test_conversions.py | 12 +- nitransforms/tests/test_surface.py | 58 +++++--- nitransforms/tests/test_version.py | 1 + nitransforms/tests/utils.py | 1 + 13 files changed, 194 insertions(+), 128 deletions(-) diff --git a/nitransforms/__init__.py b/nitransforms/__init__.py index 4ded59a0..a65fb01f 100644 --- a/nitransforms/__init__.py +++ b/nitransforms/__init__.py @@ -16,6 +16,7 @@ transform """ + from . import linear, manip, nonlinear, surface from .linear import Affine, LinearTransformsMapping from .nonlinear import DenseFieldTransform diff --git a/nitransforms/cli.py b/nitransforms/cli.py index 8f8f5ce0..aed0c9e2 100644 --- a/nitransforms/cli.py +++ b/nitransforms/cli.py @@ -32,8 +32,8 @@ def cli_apply(pargs): xfm = ( nlinload(pargs.transform, fmt=fmt) - if pargs.nonlinear else - linload(pargs.transform, fmt=fmt) + if pargs.nonlinear + else linload(pargs.transform, fmt=fmt) ) # ensure a reference is set diff --git a/nitransforms/conftest.py b/nitransforms/conftest.py index 70680882..f627fee6 100644 --- a/nitransforms/conftest.py +++ b/nitransforms/conftest.py @@ -1,4 +1,5 @@ """py.test configuration.""" + import os from pathlib import Path import numpy as np @@ -64,7 +65,9 @@ def _reorient(path): newaff = imgaff.copy() newaff[0, 0] *= -1.0 newaff[0, 3] = imgaff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[0] - _data["LAS"] = nb.Nifti1Image(np.flip(np.asanyarray(img.dataobj), 0), newaff, img.header) + _data["LAS"] = nb.Nifti1Image( + np.flip(np.asanyarray(img.dataobj), 0), newaff, img.header + ) _data["LAS"].set_data_dtype(img.get_data_dtype()) newaff = imgaff.copy() newaff[0, 0] *= -1.0 diff --git a/nitransforms/io/afni.py b/nitransforms/io/afni.py index 7c66d434..e09b68a7 100644 --- a/nitransforms/io/afni.py +++ b/nitransforms/io/afni.py @@ -1,4 +1,5 @@ """Read/write AFNI's transforms.""" + from math import pi import numpy as np from nibabel.affines import ( @@ -132,15 +133,16 @@ def to_ras(self, moving=None, reference=None): """Return a nitransforms' internal RAS matrix.""" pre_rotation = post_rotation = np.eye(4) - if reference is not None and _is_oblique(ref_aff := _ensure_image(reference).affine): + if reference is not None and _is_oblique( + ref_aff := _ensure_image(reference).affine + ): pre_rotation = _cardinal_rotation(ref_aff, True) if moving is not None and _is_oblique(mov_aff := _ensure_image(moving).affine): post_rotation = _cardinal_rotation(mov_aff, False) - return np.stack([ - post_rotation @ (xfm.to_ras() @ pre_rotation) - for xfm in self.xforms - ]) + return np.stack( + [post_rotation @ (xfm.to_ras() @ pre_rotation) for xfm in self.xforms] + ) def to_string(self): """Convert to a string directly writeable to file.""" @@ -161,7 +163,9 @@ def from_ras(cls, ras, moving=None, reference=None): pre_rotation = post_rotation = np.eye(4) - if reference is not None and _is_oblique(ref_aff := _ensure_image(reference).affine): + if reference is not None and _is_oblique( + ref_aff := _ensure_image(reference).affine + ): pre_rotation = _cardinal_rotation(ref_aff, False) if moving is not None and _is_oblique(mov_aff := _ensure_image(moving).affine): post_rotation = _cardinal_rotation(mov_aff, True) @@ -198,7 +202,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 AFNI.' % imgobj.file_map["image"].filename diff --git a/nitransforms/io/base.py b/nitransforms/io/base.py index 3c923426..c9f7d17b 100644 --- a/nitransforms/io/base.py +++ b/nitransforms/io/base.py @@ -1,4 +1,5 @@ """Read/write linear transforms.""" + from pathlib import Path import numpy as np from nibabel import load as loadimg diff --git a/nitransforms/io/fsl.py b/nitransforms/io/fsl.py index f454227e..5cf79e6c 100644 --- a/nitransforms/io/fsl.py +++ b/nitransforms/io/fsl.py @@ -1,4 +1,5 @@ """Read/write FSL's transforms.""" + import os import warnings import numpy as np @@ -41,7 +42,9 @@ def from_ras(cls, ras, moving=None, reference=None): moving = reference if reference is None: - raise TransformIOError("Cannot build FSL linear transform without a reference") + raise TransformIOError( + "Cannot build FSL linear transform without a reference" + ) reference = _ensure_image(reference) moving = _ensure_image(moving) @@ -78,7 +81,9 @@ def from_string(cls, string): def to_ras(self, moving=None, reference=None): """Return a nitransforms internal RAS+ matrix.""" if reference is None: - raise TransformIOError("Cannot build FSL linear transform without a reference") + raise TransformIOError( + "Cannot build FSL linear transform without a reference" + ) if moving is None: warnings.warn( @@ -180,10 +185,11 @@ def from_image(cls, imgobj): hdr = imgobj.header.copy() shape = hdr.get_data_shape() - if len(shape) != 4 or not shape[-1] in (2, 3): + if len(shape) != 4 or shape[-1] not in (2, 3): raise TransformFileError( - 'Displacements field "%s" does not come from FSL.' % - imgobj.file_map['image'].filename) + 'Displacements field "%s" does not come from FSL.' + % imgobj.file_map["image"].filename + ) field = np.squeeze(np.asanyarray(imgobj.dataobj)) field[..., 0] *= -1.0 diff --git a/nitransforms/io/itk.py b/nitransforms/io/itk.py index afabfd98..b2a6df8a 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 @@ -412,7 +412,9 @@ def from_h5obj(cls, fileobj, check=True, only_linear=False): # 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' + xfm[f"{typo_fallback}Parameters"], + (3, *shape.astype(int)), + order="F", ), 0, -1, @@ -422,9 +424,7 @@ def from_h5obj(cls, fileobj, check=True, only_linear=False): 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"), LPS @ affine, hdr)) continue raise TransformIOError( diff --git a/nitransforms/surface.py b/nitransforms/surface.py index 7e1e7116..0189325d 100644 --- a/nitransforms/surface.py +++ b/nitransforms/surface.py @@ -7,6 +7,7 @@ # ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## """Surface transforms.""" + import pathlib import warnings import h5py @@ -15,12 +16,10 @@ from scipy import sparse from scipy.spatial import KDTree from scipy.spatial.distance import cdist -from nitransforms.base import ( - SurfaceMesh -) +from nitransforms.base import SurfaceMesh -class SurfaceTransformBase(): +class SurfaceTransformBase: """Generic surface transformation class""" def __init__(self, reference, moving, spherical=False): @@ -64,7 +63,7 @@ def moving(self, surface): @classmethod def from_filename(cls, reference_path, moving_path): """Create an Surface Index Transformation from a pair of surfaces with corresponding - vertices.""" + vertices.""" reference = SurfaceMesh(nb.load(reference_path)) moving = SurfaceMesh(nb.load(moving_path)) return cls(reference, moving) @@ -72,9 +71,9 @@ def from_filename(cls, reference_path, moving_path): class SurfaceCoordinateTransform(SurfaceTransformBase): """Represents surface transformations in which the indices correspond and the coordinates - differ. This could be two surfaces representing difference structures from the same - hemisphere, like white matter and pial, or it could be a sphere and a deformed sphere that - moves those coordinates to a different location.""" + differ. This could be two surfaces representing difference structures from the same + hemisphere, like white matter and pial, or it could be a sphere and a deformed sphere that + moves those coordinates to a different location.""" __slots__ = ("_reference", "_moving") @@ -90,8 +89,9 @@ def __init__(self, reference, moving): super().__init__(reference=SurfaceMesh(reference), moving=SurfaceMesh(moving)) if np.all(self._reference._triangles != self._moving._triangles): - raise ValueError("Both surfaces for an index transform must have corresponding" - " vertices.") + raise ValueError( + "Both surfaces for an index transform must have corresponding vertices." + ) def map(self, x, inverse=False): if not inverse: @@ -104,8 +104,10 @@ def map(self, x, inverse=False): s_tree = KDTree(source._coords) dists, matches = s_tree.query(x) if not np.allclose(dists, 0): - raise NotImplementedError("Mapping on surfaces not implemented for coordinates that" - " aren't vertices") + raise NotImplementedError( + "Mapping on surfaces not implemented for coordinates that" + " aren't vertices" + ) return dest._coords[matches] def __add__(self, other): @@ -123,11 +125,11 @@ def _to_hdf5(self, x5_root): xform = x5_root.create_group("Transform") xform.attrs["Type"] = "SurfaceCoordinateTransform" reference = xform.create_group("Reference") - reference['Coordinates'] = h5py.SoftLink('/0/Coordinates/0') - reference['Triangles'] = h5py.SoftLink('/0/Triangles/0') + reference["Coordinates"] = h5py.SoftLink("/0/Coordinates/0") + reference["Triangles"] = h5py.SoftLink("/0/Triangles/0") moving = xform.create_group("Moving") - moving['Coordinates'] = h5py.SoftLink('/0/Coordinates/1') - moving['Triangles'] = h5py.SoftLink('/0/Triangles/0') + moving["Coordinates"] = h5py.SoftLink("/0/Coordinates/1") + moving["Triangles"] = h5py.SoftLink("/0/Triangles/0") def to_filename(self, filename, fmt=None): """Store the transform.""" @@ -148,15 +150,19 @@ def to_filename(self, filename, fmt=None): return filename @classmethod - def from_filename(cls, filename=None, reference_path=None, moving_path=None, - fmt=None): + def from_filename( + cls, filename=None, reference_path=None, moving_path=None, fmt=None + ): """Load transform from file.""" if filename is None: if reference_path is None or moving_path is None: - raise ValueError("You must pass either a X5 file or a pair of reference and moving" - " surfaces.") - return cls(SurfaceMesh(nb.load(reference_path)), - SurfaceMesh(nb.load(moving_path))) + raise ValueError( + "You must pass either a X5 file or a pair of reference and moving" + " surfaces." + ) + return cls( + SurfaceMesh(nb.load(reference_path)), SurfaceMesh(nb.load(moving_path)) + ) if fmt is None: try: @@ -175,13 +181,11 @@ def from_filename(cls, filename=None, reference_path=None, moving_path=None, assert f.attrs["Format"] == "X5" xform = f["/0/Transform"] reference = SurfaceMesh.from_arrays( - xform['Reference']['Coordinates'], - xform['Reference']['Triangles'] + xform["Reference"]["Coordinates"], xform["Reference"]["Triangles"] ) moving = SurfaceMesh.from_arrays( - xform['Moving']['Coordinates'], - xform['Moving']['Triangles'] + xform["Moving"]["Coordinates"], xform["Moving"]["Triangles"] ) return cls(reference, moving) @@ -196,9 +200,9 @@ class SurfaceResampler(SurfaceTransformBase): Then apply the transformation to sphere_unproject_from """ - __slots__ = ("_reference", "_moving", "mat", 'interpolation_method') + __slots__ = ("_reference", "_moving", "mat", "interpolation_method") - def __init__(self, reference, moving, interpolation_method='barycentric', mat=None): + def __init__(self, reference, moving, interpolation_method="barycentric", mat=None): """Initialize the resampling. Parameters @@ -218,7 +222,7 @@ def __init__(self, reference, moving, interpolation_method='barycentric', mat=No self.reference.set_radius() self.moving.set_radius() - if interpolation_method not in ['barycentric']: + if interpolation_method not in ["barycentric"]: raise NotImplementedError(f"{interpolation_method} is not implemented.") self.interpolation_method = interpolation_method @@ -267,15 +271,23 @@ def __init__(self, reference, moving, interpolation_method='barycentric', mat=No else: self.mat = sparse.csr_array(mat) # validate shape of the provided matrix - if (mat.shape[0] != moving._npoints) or (mat.shape[1] != reference._npoints): - msg = "Shape of provided mat does not match expectations based on " \ - "dimensions of moving and reference. \n" + if (mat.shape[0] != moving._npoints) or ( + mat.shape[1] != reference._npoints + ): + msg = ( + "Shape of provided mat does not match expectations based on " + "dimensions of moving and reference. \n" + ) if mat.shape[0] != moving._npoints: - msg += f" mat has {mat.shape[0]} rows but moving has {moving._npoints} " \ - f"vertices. \n" + msg += ( + f" mat has {mat.shape[0]} rows but moving has {moving._npoints} " + f"vertices. \n" + ) if mat.shape[1] != reference._npoints: - msg += f" mat has {mat.shape[1]} columns but reference has" \ - f" {reference._npoints} vertices." + msg += ( + f" mat has {mat.shape[1]} columns but reference has" + f" {reference._npoints} vertices." + ) raise ValueError(msg) def __calculate_mat(self): @@ -318,31 +330,34 @@ def map(self, x): return x def __add__(self, other): - if (isinstance(other, SurfaceResampler) - and (other.interpolation_method == self.interpolation_method)): + if isinstance(other, SurfaceResampler) and ( + other.interpolation_method == self.interpolation_method + ): return self.__class__( self.reference, other.moving, - interpolation_method=self.interpolation_method + interpolation_method=self.interpolation_method, ) raise NotImplementedError def __invert__(self): return self.__class__( - self.moving, - self.reference, - interpolation_method=self.interpolation_method + self.moving, self.reference, interpolation_method=self.interpolation_method ) @SurfaceTransformBase.reference.setter def reference(self, surface): - raise ValueError("Don't modify the reference of an existing resampling." - "Create a new one instead.") + raise ValueError( + "Don't modify the reference of an existing resampling." + "Create a new one instead." + ) @SurfaceTransformBase.moving.setter def moving(self, surface): - raise ValueError("Don't modify the moving of an existing resampling." - "Create a new one instead.") + raise ValueError( + "Don't modify the moving of an existing resampling." + "Create a new one instead." + ) def apply(self, x, inverse=False, normalize="element"): """Apply the transform to surface data. @@ -406,18 +421,18 @@ def _to_hdf5(self, x5_root): triangles.create_dataset("1", data=self.moving._triangles) xform = x5_root.create_group("Transform") xform.attrs["Type"] = "SurfaceResampling" - xform.attrs['InterpolationMethod'] = self.interpolation_method + xform.attrs["InterpolationMethod"] = self.interpolation_method mat = xform.create_group("IndexWeights") mat.create_dataset("Data", data=self.mat.data) mat.create_dataset("Indices", data=self.mat.indices) mat.create_dataset("Indptr", data=self.mat.indptr) mat.create_dataset("Shape", data=self.mat.shape) reference = xform.create_group("Reference") - reference['Coordinates'] = h5py.SoftLink('/0/Coordinates/0') - reference['Triangles'] = h5py.SoftLink('/0/Triangles/0') + reference["Coordinates"] = h5py.SoftLink("/0/Coordinates/0") + reference["Triangles"] = h5py.SoftLink("/0/Triangles/0") moving = xform.create_group("Moving") - moving['Coordinates'] = h5py.SoftLink('/0/Coordinates/1') - moving['Triangles'] = h5py.SoftLink('/0/Triangles/1') + moving["Coordinates"] = h5py.SoftLink("/0/Coordinates/1") + moving["Triangles"] = h5py.SoftLink("/0/Triangles/1") def to_filename(self, filename, fmt=None): """Store the transform.""" @@ -438,18 +453,28 @@ def to_filename(self, filename, fmt=None): return filename @classmethod - def from_filename(cls, filename=None, reference_path=None, moving_path=None, - fmt=None, interpolation_method=None): + def from_filename( + cls, + filename=None, + reference_path=None, + moving_path=None, + fmt=None, + interpolation_method=None, + ): """Load transform from file.""" if filename is None: if reference_path is None or moving_path is None: - raise ValueError("You must pass either a X5 file or a pair of reference and moving" - " surfaces.") + raise ValueError( + "You must pass either a X5 file or a pair of reference and moving" + " surfaces." + ) if interpolation_method is None: - interpolation_method = 'barycentric' - return cls(SurfaceMesh(nb.load(reference_path)), - SurfaceMesh(nb.load(moving_path)), - interpolation_method=interpolation_method) + interpolation_method = "barycentric" + return cls( + SurfaceMesh(nb.load(reference_path)), + SurfaceMesh(nb.load(moving_path)), + interpolation_method=interpolation_method, + ) if fmt is None: try: @@ -468,7 +493,7 @@ def from_filename(cls, filename=None, reference_path=None, moving_path=None, assert f.attrs["Format"] == "X5" xform = f["/0/Transform"] try: - iws = xform['IndexWeights'] + iws = xform["IndexWeights"] mat = sparse.csr_matrix( (iws["Data"][()], iws["Indices"][()], iws["Indptr"][()]), shape=iws["Shape"][()], @@ -476,43 +501,42 @@ def from_filename(cls, filename=None, reference_path=None, moving_path=None, except KeyError: mat = None reference = SurfaceMesh.from_arrays( - xform['Reference']['Coordinates'], - xform['Reference']['Triangles'] + xform["Reference"]["Coordinates"], xform["Reference"]["Triangles"] ) moving = SurfaceMesh.from_arrays( - xform['Moving']['Coordinates'], - xform['Moving']['Triangles'] + xform["Moving"]["Coordinates"], xform["Moving"]["Triangles"] ) - interpolation_method = xform.attrs['InterpolationMethod'] - return cls(reference, moving, interpolation_method=interpolation_method, mat=mat) + interpolation_method = xform.attrs["InterpolationMethod"] + return cls( + reference, moving, interpolation_method=interpolation_method, mat=mat + ) def _points_to_triangles(points, triangles): - """Implementation that vectorizes project of a point to a set of triangles. from: https://stackoverflow.com/a/32529589 """ - with np.errstate(all='ignore'): + with np.errstate(all="ignore"): # Unpack triangle points p0, p1, p2 = np.asarray(triangles).swapaxes(0, 1) # Calculate triangle edges e0 = p1 - p0 e1 = p2 - p0 - a = np.einsum('...i,...i', e0, e0) - b = np.einsum('...i,...i', e0, e1) - c = np.einsum('...i,...i', e1, e1) + a = np.einsum("...i,...i", e0, e0) + b = np.einsum("...i,...i", e0, e1) + c = np.einsum("...i,...i", e1, e1) # Calculate determinant and denominator det = a * c - b * b - inv_det = 1. / det + inv_det = 1.0 / det denom = a - 2 * b + c # Project to the edges p = p0 - points[:, np.newaxis] - d = np.einsum('...i,...i', e0, p) - e = np.einsum('...i,...i', e1, p) + d = np.einsum("...i,...i", e0, p) + e = np.einsum("...i,...i", e1, p) u = b * e - c * d v = b * d - a * e @@ -595,13 +619,15 @@ def _barycentric_weights(vecs, coords): det = coords[0] * vecs[3, 0] + coords[1] * vecs[3, 1] + coords[2] * vecs[3, 2] if det == 0: if vecs[3, 0] == 0 and vecs[3, 1] == 0 and vecs[3, 2] == 0: - warnings.warn("Zero cross product of two edges: " - "The three vertices are in the same line.") + warnings.warn( + "Zero cross product of two edges: " + "The three vertices are in the same line." + ) else: print(vecs[3]) y = coords - vecs[0] u, v = np.linalg.lstsq(vecs[1:3].T, y, rcond=None)[0] - t = 1. + t = 1.0 else: uu = coords[0] * vecs[4, 0] + coords[1] * vecs[4, 1] + coords[2] * vecs[4, 2] vv = coords[0] * vecs[5, 0] + coords[1] * vecs[5, 1] + coords[2] * vecs[5, 2] @@ -609,7 +635,7 @@ def _barycentric_weights(vecs, coords): v = vv / det tt = vecs[0, 0] * vecs[3, 0] + vecs[0, 1] * vecs[3, 1] + vecs[0, 2] * vecs[3, 2] t = tt / det - w = 1. - (u + v) + w = 1.0 - (u + v) return w, u, v, t @@ -637,7 +663,6 @@ def _find_weights(point, close_tris, d_tree): # Make sure point is actually inside triangle enclosing = True if np.all((point > closest_tri).sum(0) != 3): - enclosing = False _, ct_idxs = d_tree.query(closest_tri) a = closest_tri[0] @@ -646,7 +671,6 @@ def _find_weights(point, close_tris, d_tree): vecs = np.vstack([a, e1, e2, np.cross(e1, e2), np.cross(e2, a), np.cross(a, e1)]) res = {} res[ct_idxs[0]], res[ct_idxs[1]], res[ct_idxs[2]], _ = _barycentric_weights( - vecs, - point.squeeze() + vecs, point.squeeze() ) return res, enclosing diff --git a/nitransforms/tests/test_cli.py b/nitransforms/tests/test_cli.py index 58867131..40b911c5 100644 --- a/nitransforms/tests/test_cli.py +++ b/nitransforms/tests/test_cli.py @@ -6,8 +6,11 @@ from ..cli import cli_apply, main as ntcli if os.getenv("PYTEST_XDIST_WORKER"): - breaks_on_xdist = pytest.mark.skip(reason="xdist is active; rerun without to run this test.") + breaks_on_xdist = pytest.mark.skip( + reason="xdist is active; rerun without to run this test." + ) else: + def breaks_on_xdist(test): return test diff --git a/nitransforms/tests/test_conversions.py b/nitransforms/tests/test_conversions.py index 1229d49a..51f69aa0 100644 --- a/nitransforms/tests/test_conversions.py +++ b/nitransforms/tests/test_conversions.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: """Conversions between formats.""" + import numpy as np import pytest from .. import linear as _l @@ -109,7 +110,10 @@ def test_lta2fsl_conversions(data_path, fromto, testdata_path): ], ) def test_fsl2lta_conversions( - data_path, testdata_path, tmp_path, fromto, + data_path, + testdata_path, + tmp_path, + fromto, ): """Check conversions between formats.""" filename = f"from-{fromto[0]}_to-{fromto[1]}_mode-image" @@ -119,7 +123,7 @@ def test_fsl2lta_conversions( data_path / "regressions" / f"{filename}.fsl", reference=testdata_path / f"T1w_{fromto[0]}.nii.gz", moving=testdata_path / refname, - fmt="fsl" + fmt="fsl", ) fsl.to_filename( tmp_path / "test.lta", @@ -134,4 +138,6 @@ def test_fsl2lta_conversions( expected_fname = data_path / "regressions" / "".join((filename, ".lta")) exp_lta = LTA.from_filename(expected_fname) - assert np.allclose(converted_lta["xforms"][0]["m_L"], exp_lta["xforms"][0]["m_L"], atol=1e-4) + assert np.allclose( + converted_lta["xforms"][0]["m_L"], exp_lta["xforms"][0]["m_L"], atol=1e-4 + ) diff --git a/nitransforms/tests/test_surface.py b/nitransforms/tests/test_surface.py index a210583e..815e4838 100644 --- a/nitransforms/tests/test_surface.py +++ b/nitransforms/tests/test_surface.py @@ -8,7 +8,7 @@ from nitransforms.surface import ( SurfaceTransformBase, SurfaceCoordinateTransform, - SurfaceResampler + SurfaceResampler, ) # def test_surface_transform_npz(): @@ -42,6 +42,7 @@ # assert y_none.sum() != y_element.sum() # assert y_none.sum() != y_sum.sum() + def test_SurfaceTransformBase(testdata_path): # note these transformations are a bit of a weird use of surface transformation, but I'm # just testing the base class and the io @@ -49,7 +50,9 @@ def test_SurfaceTransformBase(testdata_path): testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_space-fsLR_desc-reg_sphere.surf.gii" ) - pial_path = testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_pial.surf.gii" + pial_path = ( + testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_pial.surf.gii" + ) sphere_reg = SurfaceMesh(nb.load(sphere_reg_path)) pial = SurfaceMesh(nb.load(pial_path)) @@ -78,7 +81,9 @@ def test_SurfaceCoordinateTransform(testdata_path): testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_space-fsLR_desc-reg_sphere.surf.gii" ) - pial_path = testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_pial.surf.gii" + pial_path = ( + testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_pial.surf.gii" + ) fslr_sphere_path = testdata_path / "tpl-fsLR_hemi-R_den-32k_sphere.surf.gii" sphere_reg = SurfaceMesh(nb.load(sphere_reg_path)) @@ -91,12 +96,15 @@ def test_SurfaceCoordinateTransform(testdata_path): # test loading from filenames sct = SurfaceCoordinateTransform(pial, sphere_reg) - sctf = SurfaceCoordinateTransform.from_filename(reference_path=pial_path, - moving_path=sphere_reg_path) + sctf = SurfaceCoordinateTransform.from_filename( + reference_path=pial_path, moving_path=sphere_reg_path + ) assert sct == sctf # test mapping - assert np.all(sct.map(sct.moving._coords[:100], inverse=True) == sct.reference._coords[:100]) + assert np.all( + sct.map(sct.moving._coords[:100], inverse=True) == sct.reference._coords[:100] + ) assert np.all(sct.map(sct.reference._coords[:100]) == sct.moving._coords[:100]) with pytest.raises(NotImplementedError): sct.map(sct.moving._coords[0]) @@ -119,7 +127,9 @@ def test_SurfaceCoordinateTransformIO(testdata_path, tmpdir): testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_space-fsLR_desc-reg_sphere.surf.gii" ) - pial_path = testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_pial.surf.gii" + pial_path = ( + testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_pial.surf.gii" + ) sct = SurfaceCoordinateTransform(pial_path, sphere_reg_path) fn = tempfile.mktemp(suffix=".h5") @@ -129,7 +139,6 @@ def test_SurfaceCoordinateTransformIO(testdata_path, tmpdir): def test_ProjectUnproject(testdata_path): - sphere_reg_path = ( testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_space-fsLR_desc-reg_sphere.surf.gii" @@ -140,10 +149,11 @@ def test_ProjectUnproject(testdata_path): / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_space-fsaverage_desc-reg_sphere.surf.gii" ) fslr_fsaverage_sphere_path = ( - testdata_path - / "tpl-fsLR_space-fsaverage_hemi-R_den-32k_sphere.surf.gii" + testdata_path / "tpl-fsLR_space-fsaverage_hemi-R_den-32k_sphere.surf.gii" + ) + pial_path = ( + testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_pial.surf.gii" ) - pial_path = testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_pial.surf.gii" # test project-unproject funcitonality projunproj = SurfaceResampler(sphere_reg_path, fslr_sphere_path) @@ -157,10 +167,7 @@ def test_ProjectUnproject(testdata_path): def test_SurfaceResampler(testdata_path, tmpdir): dif_tol = 0.001 - fslr_sphere_path = ( - testdata_path - / "tpl-fsLR_hemi-R_den-32k_sphere.surf.gii" - ) + fslr_sphere_path = testdata_path / "tpl-fsLR_hemi-R_den-32k_sphere.surf.gii" shape_path = ( testdata_path / "sub-sid000005_ses-budapest_acq-MPRAGE_hemi-R_thickness.shape.gii" @@ -190,7 +197,9 @@ def test_SurfaceResampler(testdata_path, tmpdir): moving = sphere_reg # compare results to what connectome workbench produces resampling = SurfaceResampler(reference, moving) - resampled_thickness = resampling.apply(subj_thickness.agg_data(), normalize='element') + resampled_thickness = resampling.apply( + subj_thickness.agg_data(), normalize="element" + ) ref_resampled = nb.load(ref_resampled_thickness_path).agg_data() max_dif = np.abs(resampled_thickness.astype(np.float32) - ref_resampled).max() @@ -218,13 +227,17 @@ def test_SurfaceResampler(testdata_path, tmpdir): assert np.allclose(resampling2.reference._coords, resampling.reference._coords) assert np.all(resampling2.moving._triangles == resampling.moving._triangles) - resampled_thickness2 = resampling2.apply(subj_thickness.agg_data(), normalize='element') + resampled_thickness2 = resampling2.apply( + subj_thickness.agg_data(), normalize="element" + ) assert np.all(resampled_thickness2 == resampled_thickness) # test loading with a csr assert isinstance(resampling.mat, sparse.csr_array) resampling2a = SurfaceResampler(reference, moving, mat=resampling.mat) - resampled_thickness2a = resampling2a.apply(subj_thickness.agg_data(), normalize='element') + resampled_thickness2a = resampling2a.apply( + subj_thickness.agg_data(), normalize="element" + ) assert np.all(resampled_thickness2a == resampled_thickness) with pytest.raises(ValueError): @@ -234,8 +247,11 @@ def test_SurfaceResampler(testdata_path, tmpdir): assert np.all(resampling.map(np.array([[0, 0, 0]])) == np.array([[0, 0, 0]])) # test loading from surfaces - resampling3 = SurfaceResampler.from_filename(reference_path=fslr_sphere_path, - moving_path=sphere_reg_path) + resampling3 = SurfaceResampler.from_filename( + reference_path=fslr_sphere_path, moving_path=sphere_reg_path + ) assert resampling3 == resampling - resampled_thickness3 = resampling3.apply(subj_thickness.agg_data(), normalize='element') + resampled_thickness3 = resampling3.apply( + subj_thickness.agg_data(), normalize="element" + ) assert np.all(resampled_thickness3 == resampled_thickness) diff --git a/nitransforms/tests/test_version.py b/nitransforms/tests/test_version.py index 48a70ecf..bc2a8a9f 100644 --- a/nitransforms/tests/test_version.py +++ b/nitransforms/tests/test_version.py @@ -1,4 +1,5 @@ """Test _version.py.""" + import sys from importlib import reload import nitransforms diff --git a/nitransforms/tests/utils.py b/nitransforms/tests/utils.py index 39dd8e11..e653113e 100644 --- a/nitransforms/tests/utils.py +++ b/nitransforms/tests/utils.py @@ -1,4 +1,5 @@ """Utilities for testing.""" + from pathlib import Path import numpy as np From 8bb8024aff20eec8676c0e040b9bcce1ee0037ba Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 13 Aug 2025 09:47:02 +0200 Subject: [PATCH 07/14] enh: add unit test on dense fields --- nitransforms/tests/test_nonlinear.py | 90 +++++++++++++++++++++++----- nitransforms/tests/utils.py | 27 ++++++++- 2 files changed, 101 insertions(+), 16 deletions(-) diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index 8510d993..0a208a96 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/nitransforms/tests/test_nonlinear.py @@ -13,6 +13,9 @@ BSplineFieldTransform, DenseFieldTransform, ) +from nitransforms.tests.utils import get_points + +rng = np.random.default_rng() def test_displacements_init(): @@ -74,24 +77,81 @@ def test_bsplines_references(testdata_path): ) -@pytest.mark.xfail( - reason="Disable while #266 is developed.", - strict=False, -) -def test_bspline(tmp_path, testdata_path): - """ - Cross-check B-Splines and deformation field. +@pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", "oblique"]) +@pytest.mark.parametrize("ongrid", [True, False]) +def test_densefield_map(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, rng=rng) + ) + + coords_map = grid_xyz.reshape(*shape, 3) + deltas = np.stack( + ( + np.zeros(np.prod(shape), dtype="float32").reshape(shape), + np.linspace(-80, 80, num=np.prod(shape), dtype="float32").reshape(shape), + np.linspace(-50, 50, num=np.prod(shape), dtype="float32").reshape(shape), + ), + axis=-1, + ) + + if ongrid: + atol = 1e-3 if image_orientation == "oblique" or not ongrid else 1e-7 + # Build an identity transform (deltas) + id_xfm_deltas = DenseFieldTransform(reference=reference) + np.testing.assert_array_equal(coords_map, id_xfm_deltas._field) + np.testing.assert_allclose(coords_xyz, id_xfm_deltas.map(coords_xyz)) + + # Build an identity transform (deformation) + id_xfm_field = DenseFieldTransform( + coords_map, is_deltas=False, reference=reference + ) + np.testing.assert_array_equal(coords_map, id_xfm_field._field) + np.testing.assert_allclose(coords_xyz, id_xfm_field.map(coords_xyz), atol=atol) - This test is disabled and will be split into two separate tests. - The current implementation will be moved into test_resampling.py, - since that's what it actually tests. + # Collapse to zero transform (deltas) + zero_xfm_deltas = DenseFieldTransform(-coords_map, reference=reference) + np.testing.assert_array_equal( + np.zeros_like(zero_xfm_deltas._field), zero_xfm_deltas._field + ) + np.testing.assert_allclose( + np.zeros_like(coords_xyz), zero_xfm_deltas.map(coords_xyz), atol=atol + ) - In GH-266, this test will be re-implemented by testing the equivalence - of the B-Spline and deformation field transforms by calling the - transform's `map()` method on points. + # Collapse to zero transform (deformation) + zero_xfm_field = DenseFieldTransform( + np.zeros_like(deltas), is_deltas=False, reference=reference + ) + np.testing.assert_array_equal( + np.zeros_like(zero_xfm_field._field), zero_xfm_field._field + ) + np.testing.assert_allclose( + np.zeros_like(coords_xyz), zero_xfm_field.map(coords_xyz), atol=atol + ) + + # Now let's apply a transform + xfm = DenseFieldTransform(deltas, reference=reference) + np.testing.assert_array_equal(deltas, xfm._deltas) + np.testing.assert_array_equal(coords_map + deltas, xfm._field) - """ - assert True + mapped = xfm.map(coords_xyz) + nit_deltas = mapped - coords_xyz + + if ongrid: + mapped_image = mapped.reshape(*shape, 3) + np.testing.assert_allclose(deltas + coords_map, mapped_image) + np.testing.assert_allclose(deltas, nit_deltas.reshape(*shape, 3), atol=1e-4) + np.testing.assert_allclose(xfm._field, mapped_image) + else: + ongrid_xyz = xfm.map(grid_xyz[subsample]) + assert ( + (np.linalg.norm(ongrid_xyz - mapped, axis=1) > 2).sum() + / ongrid_xyz.shape[0] + ) < 0.5 def test_map_bspline_vs_displacement(tmp_path, testdata_path): diff --git a/nitransforms/tests/utils.py b/nitransforms/tests/utils.py index e653113e..e3e8e4d9 100644 --- a/nitransforms/tests/utils.py +++ b/nitransforms/tests/utils.py @@ -2,8 +2,10 @@ from pathlib import Path import numpy as np +import nibabel as nb -from .. import linear as nbl +from nitransforms import linear as nbl +from nitransforms.base import ImageGrid def assert_affines_by_filename(affine1, affine2): @@ -26,3 +28,26 @@ def assert_affines_by_filename(affine1, affine2): xfm1 = np.loadtxt(str(affine1)) xfm2 = np.loadtxt(str(affine2)) assert np.allclose(xfm1, xfm2, atol=1e-04) + + +def get_points(reference_nii, ongrid, npoints=5000, rng=None): + """Get points in RAS space.""" + if rng is None: + rng = np.random.default_rng() + + # Get sampling indices + shape = reference_nii.shape[:3] + ref_affine = reference_nii.affine.copy() + reference = ImageGrid(nb.Nifti1Image(np.zeros(shape), ref_affine, None)) + grid_ijk = reference.ndindex + grid_xyz = reference.ras(grid_ijk) + + subsample = rng.choice(grid_ijk.shape[0], npoints) + points_ijk = grid_ijk.copy() if ongrid else grid_ijk[subsample] + coords_xyz = ( + grid_xyz + if ongrid + else reference.ras(points_ijk) + rng.normal(size=points_ijk.shape) + ) + + return coords_xyz, points_ijk, grid_xyz, shape, ref_affine, reference, subsample From f4e704ff89c02df361306b244ba5444a6ec7e6e1 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 13 Aug 2025 10:59:38 +0200 Subject: [PATCH 08/14] fix: remove apply from test_nonlinear --- nitransforms/tests/test_nonlinear.py | 14 ++------------ 1 file changed, 2 insertions(+), 12 deletions(-) diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index 0a208a96..813200d3 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/nitransforms/tests/test_nonlinear.py @@ -7,7 +7,6 @@ import numpy as np import nibabel as nb -from nitransforms.resampling import apply from nitransforms.base import TransformError from nitransforms.nonlinear import ( BSplineFieldTransform, @@ -62,17 +61,8 @@ def test_bsplines_references(testdata_path): testdata_path / "someones_bspline_coefficients.nii.gz" ).to_field() - with pytest.raises(TransformError): - apply( - BSplineFieldTransform( - testdata_path / "someones_bspline_coefficients.nii.gz" - ), - testdata_path / "someones_anatomy.nii.gz", - ) - - apply( - BSplineFieldTransform(testdata_path / "someones_bspline_coefficients.nii.gz"), - testdata_path / "someones_anatomy.nii.gz", + BSplineFieldTransform( + testdata_path / "someones_bspline_coefficients.nii.gz", reference=testdata_path / "someones_anatomy.nii.gz", ) From eeae750583920e1858eb4f19a2e1a4203f0e82cd Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 13 Aug 2025 11:00:04 +0200 Subject: [PATCH 09/14] fix: rename test --- nitransforms/tests/test_nonlinear.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index 813200d3..a1142574 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/nitransforms/tests/test_nonlinear.py @@ -144,8 +144,8 @@ def test_densefield_map(get_testdata, image_orientation, ongrid): ) < 0.5 -def test_map_bspline_vs_displacement(tmp_path, testdata_path): - """Cross-check B-Splines and deformation field.""" +def test_densefield_map_vs_bspline(tmp_path, testdata_path): + """Cross-check B-Splines and displacements field.""" os.chdir(str(tmp_path)) img_name = testdata_path / "someones_anatomy.nii.gz" From 5fda2d7425149efb332d3ea7920708ea2a15b059 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 13 Aug 2025 18:26:32 +0200 Subject: [PATCH 10/14] enh: enable displacements field tests --- nitransforms/base.py | 4 +-- nitransforms/nonlinear.py | 28 ++++++++++++++++ nitransforms/resampling.py | 1 + nitransforms/tests/test_base.py | 9 ++--- nitransforms/tests/test_resampling.py | 47 ++++++++++----------------- 5 files changed, 54 insertions(+), 35 deletions(-) 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_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.""" From 434f526bcc0259d37cc799e850b17ace9a2121ca Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 13 Aug 2025 20:45:03 +0200 Subject: [PATCH 11/14] enh: exercise ``.to_filename()`` --- nitransforms/tests/test_nonlinear.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) 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.""" From 9f5b7eb9a1567b2534adb57e54c4ddb732a793ca Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 13 Aug 2025 22:10:27 +0200 Subject: [PATCH 12/14] enh: implement ITK densefields checks vs ``antsApplyTransformsToPoints`` --- nitransforms/tests/test_io_itk.py | 149 +++++++++++++++++++++++------- 1 file changed, 114 insertions(+), 35 deletions(-) diff --git a/nitransforms/tests/test_io_itk.py b/nitransforms/tests/test_io_itk.py index 32423733..eb213e56 100644 --- a/nitransforms/tests/test_io_itk.py +++ b/nitransforms/tests/test_io_itk.py @@ -4,6 +4,7 @@ import os import shutil +from subprocess import check_call import pytest @@ -14,10 +15,19 @@ from scipy.io import loadmat from h5py import File as H5File -from nitransforms.io.base import TransformIOError, TransformFileError -from nitransforms.io import itk +from nitransforms.io.base import TransformIOError, TransformFileError +from nitransforms.io.itk import ( + ITKLinearTransform, + ITKLinearTransformArray, + ITKDisplacementsField, + ITKCompositeH5, +) +from nitransforms import nonlinear as nitnl from nitransforms.conftest import _datadir, _testdir +from nitransforms.tests.utils import get_points + +rng = np.random.default_rng() LPS = np.diag([-1, -1, 1, 1]) ITK_MAT = LPS.dot(np.ones((4, 4)).dot(LPS)) @@ -29,14 +39,14 @@ def test_ITKLinearTransform(tmpdir, testdata_path): matlabfile = testdata_path / "ds-005_sub-01_from-T1_to-OASIS_affine.mat" mat = loadmat(str(matlabfile)) with open(str(matlabfile), "rb") as f: - itkxfm = itk.ITKLinearTransform.from_fileobj(f) + itkxfm = ITKLinearTransform.from_fileobj(f) assert np.allclose( itkxfm["parameters"][:3, :3].flatten(), mat["AffineTransform_float_3_3"][:-3].flatten(), ) assert np.allclose(itkxfm["offset"], mat["fixed"].reshape((3,))) - itkxfm = itk.ITKLinearTransform.from_filename(matlabfile) + itkxfm = ITKLinearTransform.from_filename(matlabfile) assert np.allclose( itkxfm["parameters"][:3, :3].flatten(), mat["AffineTransform_float_3_3"][:-3].flatten(), @@ -46,17 +56,17 @@ def test_ITKLinearTransform(tmpdir, testdata_path): # Test to_filename(textfiles) itkxfm.to_filename("textfile.tfm") with open("textfile.tfm") as f: - itkxfm2 = itk.ITKLinearTransform.from_fileobj(f) + itkxfm2 = ITKLinearTransform.from_fileobj(f) assert np.allclose(itkxfm["parameters"], itkxfm2["parameters"]) # Test to_filename(matlab) itkxfm.to_filename("copy.mat") with open("copy.mat", "rb") as f: - itkxfm3 = itk.ITKLinearTransform.from_fileobj(f) + itkxfm3 = ITKLinearTransform.from_fileobj(f) assert np.all(itkxfm["parameters"] == itkxfm3["parameters"]) rasmat = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0]) - itkxfm = itk.ITKLinearTransform.from_ras(rasmat) + itkxfm = ITKLinearTransform.from_ras(rasmat) assert np.allclose(itkxfm["parameters"], ITK_MAT * rasmat) assert np.allclose(itkxfm.to_ras(), rasmat) @@ -67,9 +77,9 @@ def test_ITKLinearTransformArray(tmpdir, data_path): with open(str(data_path / "itktflist.tfm")) as f: text = f.read() f.seek(0) - itklist = itk.ITKLinearTransformArray.from_fileobj(f) + itklist = ITKLinearTransformArray.from_fileobj(f) - itklistb = itk.ITKLinearTransformArray.from_filename(data_path / "itktflist.tfm") + itklistb = ITKLinearTransformArray.from_filename(data_path / "itktflist.tfm") assert itklist["nxforms"] == itklistb["nxforms"] assert all( [ @@ -85,14 +95,14 @@ def test_ITKLinearTransformArray(tmpdir, data_path): assert itklist["nxforms"] == 9 assert text == itklist.to_string() with pytest.raises(TransformFileError): - itk.ITKLinearTransformArray.from_string("\n".join(text.splitlines()[1:])) + ITKLinearTransformArray.from_string("\n".join(text.splitlines()[1:])) itklist.to_filename("copy.tfm") with open("copy.tfm") as f: copytext = f.read() assert text == copytext - itklist = itk.ITKLinearTransformArray( + itklist = ITKLinearTransformArray( xforms=[np.random.normal(size=(4, 4)) for _ in range(4)] ) @@ -108,7 +118,7 @@ def test_ITKLinearTransformArray(tmpdir, data_path): f.write(xfm.to_string()) with open("extracted.tfm") as f: - xfm2 = itk.ITKLinearTransform.from_fileobj(f) + xfm2 = ITKLinearTransform.from_fileobj(f) assert np.allclose( xfm.structarr["parameters"][:3, ...], xfm2.structarr["parameters"][:3, ...] ) @@ -118,18 +128,18 @@ def test_ITKLinearTransformArray(tmpdir, data_path): itklist.to_filename("matlablist.mat") with pytest.raises(TransformFileError): - xfm2 = itk.ITKLinearTransformArray.from_binary({}) + xfm2 = ITKLinearTransformArray.from_binary({}) with open("filename.mat", "ab") as f: with pytest.raises(TransformFileError): - xfm2 = itk.ITKLinearTransformArray.from_fileobj(f) + xfm2 = ITKLinearTransformArray.from_fileobj(f) @pytest.mark.parametrize("size", [(20, 20, 20), (20, 20, 20, 3)]) def test_itk_disp_load(size): """Checks field sizes.""" with pytest.raises(TransformFileError): - itk.ITKDisplacementsField.from_image( + ITKDisplacementsField.from_image( nb.Nifti1Image(np.zeros(size), np.eye(4), None) ) @@ -137,7 +147,7 @@ def test_itk_disp_load(size): def test_itk_disp_load_intent(): """Checks whether the NIfTI intent is fixed.""" with pytest.warns(UserWarning): - field = itk.ITKDisplacementsField.from_image( + field = ITKDisplacementsField.from_image( nb.Nifti1Image(np.zeros((20, 20, 20, 1, 3)), np.eye(4), None) ) @@ -161,7 +171,7 @@ def test_itk_h5(tmpdir, only_linear, h5_path, nxforms): assert ( len( list( - itk.ITKCompositeH5.from_filename( + ITKCompositeH5.from_filename( h5_path, only_linear=only_linear, ) @@ -174,7 +184,7 @@ def test_itk_h5(tmpdir, only_linear, h5_path, nxforms): with pytest.raises(TransformFileError): list( - itk.ITKCompositeH5.from_filename( + ITKCompositeH5.from_filename( h5_path.absolute().name.replace(".h5", ".x5"), only_linear=only_linear, ) @@ -190,30 +200,28 @@ def test_itk_h5(tmpdir, only_linear, h5_path, nxforms): xfm["TransformType"][0] = b"InventTransform" with pytest.raises(TransformIOError): - itk.ITKCompositeH5.from_filename("test.h5") + ITKCompositeH5.from_filename("test.h5") def test_itk_linear_h5(tmpdir, data_path, testdata_path): """Check different lower-level loading options.""" # File loadable with transform array - h5xfm = itk.ITKLinearTransformArray.from_filename( - data_path / "affine-antsComposite.h5" - ) + h5xfm = ITKLinearTransformArray.from_filename(data_path / "affine-antsComposite.h5") assert len(h5xfm.xforms) == 1 with open(data_path / "affine-antsComposite.h5", "rb") as f: - h5xfm = itk.ITKLinearTransformArray.from_fileobj(f) + h5xfm = ITKLinearTransformArray.from_fileobj(f) assert len(h5xfm.xforms) == 1 # File loadable with single affine object - itk.ITKLinearTransform.from_filename(data_path / "affine-antsComposite.h5") + ITKLinearTransform.from_filename(data_path / "affine-antsComposite.h5") with open(data_path / "affine-antsComposite.h5", "rb") as f: - itk.ITKLinearTransform.from_fileobj(f) + ITKLinearTransform.from_fileobj(f) # Exercise only_linear - itk.ITKCompositeH5.from_filename( + ITKCompositeH5.from_filename( testdata_path / "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5", only_linear=True, @@ -231,16 +239,16 @@ def test_itk_linear_h5(tmpdir, data_path, testdata_path): xfm["TransformFixedParameters"] = np.zeros(3, dtype=float) # File loadable with transform array - h5xfm = itk.ITKLinearTransformArray.from_filename("test.h5") + h5xfm = ITKLinearTransformArray.from_filename("test.h5") assert len(h5xfm.xforms) == 2 # File loadable with generalistic object (NOTE we directly access the list) - h5xfm = itk.ITKCompositeH5.from_filename("test.h5") + h5xfm = ITKCompositeH5.from_filename("test.h5") assert len(h5xfm) == 2 # Error raised if the we try to use the single affine loader with pytest.raises(TransformIOError): - itk.ITKLinearTransform.from_filename("test.h5") + ITKLinearTransform.from_filename("test.h5") shutil.copy(data_path / "affine-antsComposite.h5", "test.h5") os.chmod("test.h5", 0o666) @@ -251,12 +259,12 @@ def test_itk_linear_h5(tmpdir, data_path, testdata_path): del h5group["1"] # File loadable with generalistic object - h5xfm = itk.ITKCompositeH5.from_filename("test.h5") + h5xfm = ITKCompositeH5.from_filename("test.h5") assert len(h5xfm) == 0 # Error raised if the we try to use the single affine loader with pytest.raises(TransformIOError): - itk.ITKLinearTransform.from_filename("test.h5") + ITKLinearTransform.from_filename("test.h5") # Added tests for h5 orientation bug @@ -285,7 +293,7 @@ def test_itk_h5_field_order(tmp_path): g1["TransformFixedParameters"] = fixed g1["TransformParameters"] = params - img = itk.ITKCompositeH5.from_filename(fname)[0] + img = ITKCompositeH5.from_filename(fname)[0] expected = np.moveaxis(field, 0, -1) expected[..., (0, 1)] *= -1 assert np.allclose(img.get_fdata(), expected) @@ -314,9 +322,9 @@ def _load_composite_testdata(data_path): def test_itk_h5_displacement_mismatch(testdata_path): """Composite displacements should match the standalone field""" h5file, warpfile = _load_composite_testdata(testdata_path) - xforms = itk.ITKCompositeH5.from_filename(h5file) + xforms = ITKCompositeH5.from_filename(h5file) field_h5 = xforms[1] - field_img = itk.ITKDisplacementsField.from_filename(warpfile) + field_img = ITKDisplacementsField.from_filename(warpfile) np.testing.assert_array_equal( np.asanyarray(field_h5.dataobj), np.asanyarray(field_img.dataobj) @@ -368,7 +376,78 @@ def test_itk_h5_field_order_fortran(tmp_path): g1["TransformFixedParameters"] = fixed g1["TransformParameters"] = params - img = itk.ITKCompositeH5.from_filename(fname)[0] + img = ITKCompositeH5.from_filename(fname)[0] expected = np.moveaxis(field, 0, -1) expected[..., (0, 1)] *= -1 assert np.allclose(img.get_fdata(), expected) + + +# Tests against ANTs' ``antsApplyTransformsToPoints`` +@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" + + # antsApplyTransformsToPoints wants LPS coordinates, see last post at + # http://sourceforge.net/p/advants/discussion/840261/thread/2a1e9307/ + lps_xyz = coords_xyz.copy() * (-1, -1, 1) + np.savetxt(csvin, lps_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) + + # antsApplyTransformsToPoints writes LPS coordinates, see last post at + # http://sourceforge.net/p/advants/discussion/840261/thread/2a1e9307/ + ants_pts = np.vstack([ants_res[n] for n in ("x", "y", "z")]).T * (-1, -1, 1) + + xfm = nitnl.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.load(warpfile).to_filename(tmp_path / "original_ants_deltas.nii.gz") + + nb.Nifti1Image(coords_map, ref_affine, None).to_filename( + tmp_path / "baseline_positions.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" + ) + 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" + ) + + # When transforming off-grid points, rounding errors are large + atol = 0 if ongrid else 1e-1 + rtol = 1e-4 if ongrid else 1e-3 + np.testing.assert_allclose(mapped, ants_pts, atol=atol, rtol=rtol) From d48a739318678da97f96843eba1488ee91bf4d44 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 14 Aug 2025 06:59:26 +0200 Subject: [PATCH 13/14] enh: add test with constant field --- nitransforms/tests/test_io_itk.py | 106 +++++++++++++++++++++++++++++- 1 file changed, 104 insertions(+), 2 deletions(-) diff --git a/nitransforms/tests/test_io_itk.py b/nitransforms/tests/test_io_itk.py index eb213e56..54f2e29f 100644 --- a/nitransforms/tests/test_io_itk.py +++ b/nitransforms/tests/test_io_itk.py @@ -27,8 +27,7 @@ from nitransforms.conftest import _datadir, _testdir from nitransforms.tests.utils import get_points -rng = np.random.default_rng() - +RNG_SEED = 202508140819 LPS = np.diag([-1, -1, 1, 1]) ITK_MAT = LPS.dot(np.ones((4, 4)).dot(LPS)) @@ -386,6 +385,8 @@ def test_itk_h5_field_order_fortran(tmp_path): @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.""" + + rng = np.random.default_rng(RNG_SEED) warpfile = ( testdata_path / "regressions" @@ -451,3 +452,104 @@ def test_densefield_map_vs_ants(testdata_path, tmp_path, ongrid): atol = 0 if ongrid else 1e-1 rtol = 1e-4 if ongrid else 1e-3 np.testing.assert_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.""" + + rng = np.random.default_rng(RNG_SEED) + + 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) + ) + + tol = ( + {"atol": 0, "rtol": 1e-4} + if image_orientation != "oblique" + else {"atol": 1e-4, "rtol": 1e-2} + ) + coords_map = grid_xyz.reshape(*shape, 3) + + 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,)) + gold_mapped_xyz = coords_map + deltas + + fieldnii = nb.Nifti1Image(deltas, ref_affine, None) + warpfile = tmp_path / "itk_transform.nii.gz" + + # Ensure direct (xfm) and ITK roundtrip (itk_xfm) are equivalent + xfm = nitnl.DenseFieldTransform(fieldnii) + xfm.to_filename(warpfile, fmt="itk") + itk_xfm = nitnl.DenseFieldTransform(ITKDisplacementsField.from_filename(warpfile)) + + 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, **tol) + if image_orientation != "oblique": + assert xfm == itk_xfm + + # Ensure deltas and mapped grid are equivalent + orig_grid_mapped_xyz = xfm.map(grid_xyz).reshape(*shape, -1) + np.testing.assert_allclose(gold_mapped_xyz, orig_grid_mapped_xyz) + + # Test ANTs mapping + grid_mapped_xyz = itk_xfm.map(grid_xyz).reshape(*shape, -1) + + # Check apparent healthiness of mapping + np.testing.assert_allclose(gold_mapped_xyz, grid_mapped_xyz, **tol) + np.testing.assert_allclose(orig_grid_mapped_xyz, grid_mapped_xyz, **tol) + + csvout = tmp_path / "mapped_xyz.csv" + csvin = tmp_path / "coords_xyz.csv" + # antsApplyTransformsToPoints wants LPS coordinates, see last post at + # http://sourceforge.net/p/advants/discussion/840261/thread/2a1e9307/ + lps_xyz = coords_xyz.copy() * (-1, -1, 1) + np.savetxt(csvin, lps_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) + # antsApplyTransformsToPoints writes LPS coordinates, see last post at + # http://sourceforge.net/p/advants/discussion/840261/thread/2a1e9307/ + ants_pts = np.vstack([ants_res[n] for n in ("x", "y", "z")]).T * (-1, -1, 1) + + 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_allclose(gold_mapped_xyz, ants_pts, rtol=1e-2, atol=1e-3) + np.testing.assert_allclose(deltas, ants_pts - coords_map, rtol=1e-2, atol=1e-3) + else: + # TODO Change test to norms and investigate extreme cases + # We're likely hitting OBB points (see gh-188) + # https://github.com/nipy/nitransforms/pull/188 + np.testing.assert_allclose( + xfm.map(coords_xyz) - coords_xyz, ants_pts - coords_xyz, rtol=1, atol=1 + ) From 7772372ead9da388493cef3e8d9f5a28da8baefe Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 14 Aug 2025 10:10:51 +0200 Subject: [PATCH 14/14] rel(25.0.1): Update CHANGES [skip ci] --- CHANGES.rst | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index 126b804e..ab5b1a92 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,23 @@ +25.0.1 (August 14, 2025) +======================== +A patch release addressing a critical issue in the ``ImageGrid`` class relating coordinates, +and bolstered the testing of ANTs' generated displacements fields by cross-comparing against +``antsApplyTransformsToPoints`` in several new tests. + +CHANGES +------- +* FIX: ``ImageGrid._coords`` was somehow overwritten + re-enable tests by @oesteban in https://github.com/nipy/nitransforms/pull/276 +* FIX: Revision of index and RAS coordinate grids generation by @oesteban in https://github.com/nipy/nitransforms/pull/271 +* ENH: Implement ITK densefields checks vs ``antsApplyTransformsToPoints`` by @oesteban in https://github.com/nipy/nitransforms/pull/277 +* ENH: Add unit test on dense fields (extracted from #266) by @oesteban in https://github.com/nipy/nitransforms/pull/274 +* RF: Move tests to better locations by @oesteban in https://github.com/nipy/nitransforms/pull/272 +* MNT: Minimal housekeeping of tests by @oesteban in https://github.com/nipy/nitransforms/pull/275 +* MNT: Fix coverage XML path in CircleCI by @oesteban in https://github.com/nipy/nitransforms/pull/265 +* MNT: Add test cases demonstrating ordering bug reading composite ITK's HDF5 files by @oesteban in https://github.com/nipy/nitransforms/pull/263 +* STY: Run ruff at the source root by @oesteban in https://github.com/nipy/nitransforms/pull/273 + +**Full Changelog**: https://github.com/nipy/nitransforms/compare/25.0.0...25.0.1 + 25.0.0 (July 22, 2025) ====================== A new major release introducing critical fixes and important new functionality. @@ -10,7 +30,7 @@ CHANGES ------- * FIX: BSpline mapping of individual points by @oesteban in https://github.com/nipy/nitransforms/pull/256 * FIX: Remove implementation of an abstract class by @oesteban in https://github.com/nipy/nitransforms/pull/255 -* FIX: Add test for `DenseFieldTransform` handling of OOB points by @oesteban in https://github.com/nipy/nitransforms/pull/254 +* FIX: Add test for ``DenseFieldTransform`` handling of OOB points by @oesteban in https://github.com/nipy/nitransforms/pull/254 * ENH: X5 read/write support of ``TransformChain`` by @oesteban in https://github.com/nipy/nitransforms/pull/253 * ENH: Loading of X5 (linear) transforms by @oesteban in https://github.com/nipy/nitransforms/pull/243 * ENH: Implement X5 representation and output to filesystem by @oesteban in https://github.com/nipy/nitransforms/pull/241