Merge pull request #90 from oesteban/enh/88-collapse-linear-xfms

ENH: Add an ``.asaffine()`` member to ``TransformChain``

Author: oesteban

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):

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):

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)

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)

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,
+    )