diff --git a/nitransforms/linear.py b/nitransforms/linear.py index f09aae72..1b8d4040 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -24,7 +24,7 @@ class Affine(TransformBase): """Represents linear transforms on image data.""" - __slots__ = ("_matrix", ) + __slots__ = ("_matrix", "_inverse") def __init__(self, matrix=None, reference=None): """ @@ -57,6 +57,7 @@ def __init__(self, matrix=None, reference=None): """ super().__init__(reference=reference) self._matrix = np.eye(4) + self._inverse = np.eye(4) if matrix is not None: matrix = np.array(matrix) @@ -72,6 +73,7 @@ def __init__(self, matrix=None, reference=None): # Normalize last row self._matrix[3, :] = (0, 0, 0, 1) + self._inverse = np.linalg.inv(self._matrix) def __eq__(self, other): """ @@ -90,6 +92,44 @@ def __eq__(self, other): warnings.warn("Affines are equal, but references do not match.") return _eq + def __invert__(self): + """ + Get the inverse of this transform. + + Example + ------- + >>> matrix = [[1, 0, 0, 4], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]] + >>> Affine(np.linalg.inv(matrix)) == ~Affine(matrix) + True + + """ + return self.__class__(self._inverse) + + def __matmul__(self, b): + """ + Compose two Affines. + + Example + ------- + >>> xfm1 = Affine([[1, 0, 0, 4], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) + >>> xfm1 @ ~xfm1 == Affine() + True + + >>> xfm1 = Affine([[1, 0, 0, 4], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]) + >>> xfm1 @ np.eye(4) == xfm1 + True + + """ + if not isinstance(b, self.__class__): + _b = self.__class__(b) + else: + _b = b + + retval = self.__class__(self.matrix.dot(_b.matrix)) + if _b.reference: + retval.reference = _b.reference + return retval + @property def matrix(self): """Access the internal representation of this affine.""" @@ -124,14 +164,14 @@ def map(self, x, inverse=False): affine = self._matrix coords = _as_homogeneous(x, dim=affine.shape[0] - 1).T if inverse is True: - affine = np.linalg.inv(self._matrix) + affine = self._inverse return affine.dot(coords).T[..., :-1] def _to_hdf5(self, x5_root): """Serialize this object into the x5 file format.""" xform = x5_root.create_dataset("Transform", data=[self._matrix]) xform.attrs["Type"] = "affine" - x5_root.create_dataset("Inverse", data=[np.linalg.inv(self._matrix)]) + x5_root.create_dataset("Inverse", data=[(~self).matrix]) if self._reference: self.reference._to_hdf5(x5_root.create_group("Reference")) @@ -175,7 +215,7 @@ def to_filename(self, filename, fmt="X5", moving=None): lt["dst"] = io.VolumeGeometry.from_image(moving) # However, the affine needs to be inverted # (i.e., it is not a pure "points" convention). - lt["m_L"] = np.linalg.inv(self.matrix) + lt["m_L"] = (~self).matrix # to make LTA file format lta = io.LinearTransformArray() lta["type"] = 1 # RAS2RAS @@ -234,6 +274,11 @@ def __init__(self, transforms, reference=None): [0., 1., 0., 2.], [0., 0., 1., 3.], [0., 0., 0., 1.]]) + >>> (~xfm)[0].matrix # doctest: +NORMALIZE_WHITESPACE + array([[ 1., 0., 0., -1.], + [ 0., 1., 0., -2.], + [ 0., 0., 1., -3.], + [ 0., 0., 0., 1.]]) """ super().__init__(reference=reference) @@ -245,6 +290,7 @@ def __init__(self, transforms, reference=None): ).matrix for xfm in transforms ], axis=0) + self._inverse = np.linalg.inv(self._matrix) def __getitem__(self, i): """Enable indexed access to the series of matrices.""" @@ -304,7 +350,7 @@ def map(self, x, inverse=False): affine = self.matrix coords = _as_homogeneous(x, dim=affine.shape[-1] - 1).T if inverse is True: - affine = np.linalg.inv(affine) + affine = self._inverse return np.swapaxes(affine.dot(coords), 1, 2) def to_filename(self, filename, fmt="X5", moving=None): diff --git a/nitransforms/manip.py b/nitransforms/manip.py index b3f3cb2d..76a9f90f 100644 --- a/nitransforms/manip.py +++ b/nitransforms/manip.py @@ -139,6 +139,13 @@ def map(self, x, inverse=False): return x + def asaffine(self): + """Combine a succession of linear transforms into one.""" + retval = self.transforms[-1] + for xfm in self.transforms[:-1][::-1]: + retval @= xfm + return retval + @classmethod def from_filename(cls, filename, fmt="X5", reference=None, moving=None): diff --git a/nitransforms/tests/test_base.py b/nitransforms/tests/test_base.py index ff92fcba..3d096f22 100644 --- a/nitransforms/tests/test_base.py +++ b/nitransforms/tests/test_base.py @@ -5,6 +5,7 @@ import h5py from ..base import SpatialReference, SampledSpatialData, ImageGrid, TransformBase +from .. import linear as nitl def test_SpatialReference(testdata_path): @@ -134,3 +135,11 @@ def test_SampledSpatialData(testdata_path): with pytest.raises(TypeError): gii = nb.gifti.GiftiImage() SampledSpatialData(gii) + + +def test_concatenation(testdata_path): + """Check concatenation of affines.""" + aff = nitl.Affine(reference=testdata_path / 'someones_anatomy.nii.gz') + x = [(0., 0., 0.), (1., 1., 1.), (-1., -1., -1.)] + assert np.all((aff + nitl.Affine())(x) == x) + assert np.all((aff + nitl.Affine())(x, inverse=True) == x) diff --git a/nitransforms/tests/test_linear.py b/nitransforms/tests/test_linear.py index 2a638b51..3045050d 100644 --- a/nitransforms/tests/test_linear.py +++ b/nitransforms/tests/test_linear.py @@ -11,7 +11,7 @@ import nibabel as nb from nibabel.eulerangles import euler2mat from nibabel.affines import from_matvec -from .. import linear as ntl +from .. import linear as nitl from .utils import assert_affines_by_filename TESTS_BORDER_TOLERANCE = 0.05 @@ -42,35 +42,35 @@ def test_linear_typeerrors1(matrix): """Exercise errors in Affine creation.""" with pytest.raises(TypeError): - ntl.Affine(matrix) + nitl.Affine(matrix) def test_linear_typeerrors2(data_path): """Exercise errors in Affine creation.""" with pytest.raises(TypeError): - ntl.Affine.from_filename(data_path / 'itktflist.tfm', fmt='itk') + nitl.Affine.from_filename(data_path / 'itktflist.tfm', fmt='itk') def test_linear_valueerror(): """Exercise errors in Affine creation.""" with pytest.raises(ValueError): - ntl.Affine(np.ones((4, 4))) + nitl.Affine(np.ones((4, 4))) def test_loadsave_itk(tmp_path, data_path, testdata_path): """Test idempotency.""" ref_file = testdata_path / 'someones_anatomy.nii.gz' - xfm = ntl.load(data_path / 'itktflist2.tfm', fmt='itk') - assert isinstance(xfm, ntl.LinearTransformsMapping) + xfm = nitl.load(data_path / 'itktflist2.tfm', fmt='itk') + assert isinstance(xfm, nitl.LinearTransformsMapping) xfm.reference = ref_file xfm.to_filename(tmp_path / 'transform-mapping.tfm', fmt='itk') assert (data_path / 'itktflist2.tfm').read_text() \ == (tmp_path / 'transform-mapping.tfm').read_text() - single_xfm = ntl.load(data_path / 'affine-LAS.itk.tfm', fmt='itk') - assert isinstance(single_xfm, ntl.Affine) - assert single_xfm == ntl.Affine.from_filename( + single_xfm = nitl.load(data_path / 'affine-LAS.itk.tfm', fmt='itk') + assert isinstance(single_xfm, nitl.Affine) + assert single_xfm == nitl.Affine.from_filename( data_path / 'affine-LAS.itk.tfm', fmt='itk') @@ -79,23 +79,23 @@ def test_loadsave_itk(tmp_path, data_path, testdata_path): def test_loadsave(tmp_path, data_path, testdata_path, fmt): """Test idempotency.""" ref_file = testdata_path / 'someones_anatomy.nii.gz' - xfm = ntl.load(data_path / 'itktflist2.tfm', fmt='itk') + xfm = nitl.load(data_path / 'itktflist2.tfm', fmt='itk') xfm.reference = ref_file fname = tmp_path / '.'.join(('transform-mapping', fmt)) xfm.to_filename(fname, fmt=fmt) - xfm == ntl.load(fname, fmt=fmt, reference=ref_file) + xfm == nitl.load(fname, fmt=fmt, reference=ref_file) xfm.to_filename(fname, fmt=fmt, moving=ref_file) - xfm == ntl.load(fname, fmt=fmt, reference=ref_file) + xfm == nitl.load(fname, fmt=fmt, reference=ref_file) ref_file = testdata_path / 'someones_anatomy.nii.gz' - xfm = ntl.load(data_path / 'affine-LAS.itk.tfm', fmt='itk') + xfm = nitl.load(data_path / 'affine-LAS.itk.tfm', fmt='itk') xfm.reference = ref_file fname = tmp_path / '.'.join(('single-transform', fmt)) xfm.to_filename(fname, fmt=fmt) - xfm == ntl.load(fname, fmt=fmt, reference=ref_file) + xfm == nitl.load(fname, fmt=fmt, reference=ref_file) xfm.to_filename(fname, fmt=fmt, moving=ref_file) - xfm == ntl.load(fname, fmt=fmt, reference=ref_file) + xfm == nitl.load(fname, fmt=fmt, reference=ref_file) @pytest.mark.xfail(reason="Not fully implemented") @@ -107,7 +107,7 @@ def test_linear_save(tmpdir, data_path, get_testdata, image_orientation, sw_tool img = get_testdata[image_orientation] # Generate test transform T = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0]) - xfm = ntl.Affine(T) + xfm = nitl.Affine(T) xfm.reference = img ext = '' @@ -140,7 +140,7 @@ def test_apply_linear_transform( img = get_testdata[image_orientation] # Generate test transform T = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0]) - xfm = ntl.Affine(T) + xfm = nitl.Affine(T) xfm.reference = img ext = '' @@ -172,11 +172,16 @@ def test_apply_linear_transform( # A certain tolerance is necessary because of resampling at borders assert (np.abs(diff) > 1e-3).sum() / diff.size < TESTS_BORDER_TOLERANCE + nt_moved = xfm.apply('img.nii.gz', order=0) + diff = sw_moved.get_fdata() - nt_moved.get_fdata() + # A certain tolerance is necessary because of resampling at borders + assert (np.abs(diff) > 1e-3).sum() / diff.size < TESTS_BORDER_TOLERANCE + def test_Affine_to_x5(tmpdir, testdata_path): """Test affine's operations.""" tmpdir.chdir() - aff = ntl.Affine() + aff = nitl.Affine() with h5py.File('xfm.x5', 'w') as f: aff._to_hdf5(f.create_group('Affine')) @@ -185,34 +190,42 @@ def test_Affine_to_x5(tmpdir, testdata_path): aff._to_hdf5(f.create_group('Affine')) -def test_concatenation(testdata_path): - """Check concatenation of affines.""" - aff = ntl.Affine(reference=testdata_path / 'someones_anatomy.nii.gz') - x = [(0., 0., 0.), (1., 1., 1.), (-1., -1., -1.)] - assert np.all((aff + ntl.Affine())(x) == x) - assert np.all((aff + ntl.Affine())(x, inverse=True) == x) - - def test_LinearTransformsMapping_apply(tmp_path, data_path, testdata_path): """Apply transform mappings.""" - hmc = ntl.load(data_path / 'hmc-itk.tfm', fmt='itk', - reference=testdata_path / 'sbref.nii.gz') - assert isinstance(hmc, ntl.LinearTransformsMapping) + hmc = nitl.load(data_path / 'hmc-itk.tfm', fmt='itk', + reference=testdata_path / 'sbref.nii.gz') + assert isinstance(hmc, nitl.LinearTransformsMapping) - # Test-case: realing functional data on to sbref + # Test-case: realign functional data on to sbref nii = hmc.apply(testdata_path / 'func.nii.gz', order=1, reference=testdata_path / 'sbref.nii.gz') assert nii.dataobj.shape[-1] == len(hmc) # Test-case: write out a fieldmap moved with head - hmcinv = ntl.LinearTransformsMapping( + hmcinv = nitl.LinearTransformsMapping( np.linalg.inv(hmc.matrix), reference=testdata_path / 'func.nii.gz') nii = hmcinv.apply(testdata_path / 'fmap.nii.gz', order=1) assert nii.dataobj.shape[-1] == len(hmc) # Ensure a ValueError is issued when trying to do weird stuff - hmc = ntl.LinearTransformsMapping(hmc.matrix[:1, ...]) + hmc = nitl.LinearTransformsMapping(hmc.matrix[:1, ...]) with pytest.raises(ValueError): hmc.apply(testdata_path / 'func.nii.gz', order=1, reference=testdata_path / 'sbref.nii.gz') + + +def test_mulmat_operator(testdata_path): + """Check the @ operator.""" + ref = testdata_path / 'someones_anatomy.nii.gz' + mat1 = np.diag([2., 2., 2., 1.]) + mat2 = from_matvec(np.eye(3), (4, 2, -1)) + aff = nitl.Affine(mat1, reference=ref) + + composed = aff @ mat2 + assert composed.reference is None + assert composed == nitl.Affine(mat1.dot(mat2)) + + composed = nitl.Affine(mat2) @ aff + assert composed.reference == aff.reference + assert composed == nitl.Affine(mat2.dot(mat1), reference=ref) diff --git a/nitransforms/tests/test_manip.py b/nitransforms/tests/test_manip.py index 5bfc6298..8b387bba 100644 --- a/nitransforms/tests/test_manip.py +++ b/nitransforms/tests/test_manip.py @@ -8,12 +8,15 @@ import numpy as np import nibabel as nb -from ..manip import load as _load +from ..manip import load as _load, TransformChain +from ..linear import Affine from .test_nonlinear import ( TESTS_BORDER_TOLERANCE, APPLY_NONLINEAR_CMD, ) +FMT = {"lta": "fs", "tfm": "itk"} + def test_itk_h5(tmp_path, testdata_path): """Check a translation-only field on one or more axes, different image orientations.""" @@ -51,3 +54,22 @@ def test_itk_h5(tmp_path, testdata_path): diff = sw_moved.get_fdata() - nt_moved.get_fdata() # A certain tolerance is necessary because of resampling at borders assert (np.abs(diff) > 1e-3).sum() / diff.size < TESTS_BORDER_TOLERANCE + + +@pytest.mark.parametrize("ext0", ["lta", "tfm"]) +@pytest.mark.parametrize("ext1", ["lta", "tfm"]) +@pytest.mark.parametrize("ext2", ["lta", "tfm"]) +def test_collapse_affines(tmp_path, data_path, ext0, ext1, ext2): + """Check whether affines are correctly collapsed.""" + chain = TransformChain([ + Affine.from_filename(data_path / "regressions" + / f"from-fsnative_to-scanner_mode-image.{ext0}", fmt=f"{FMT[ext0]}"), + Affine.from_filename(data_path / "regressions" + / f"from-scanner_to-bold_mode-image.{ext1}", fmt=f"{FMT[ext1]}"), + ]) + assert np.allclose( + chain.asaffine().matrix, + Affine.from_filename( + data_path / "regressions" / f"from-fsnative_to-bold_mode-image.{ext2}", + fmt=f"{FMT[ext2]}").matrix, + )