Merge pull request #66 from oesteban/test/lta-itk-conversion

oesteban · web-flow · commit 7a9984173abd · 2020-03-19T15:33:55.000-07:00
ENH: Add ITK-LTA conversion test
diff --git a/nitransforms/io/lta.py b/nitransforms/io/lta.py
@@ -160,9 +160,38 @@ def set_type(self, new_type):
         )
 
     def to_ras(self, moving=None, reference=None):
-        """Return a nitransforms internal RAS+ matrix."""
+        """
+        Return a nitransforms' internal RAS+ array.
+
+        Seemingly, the matrix of an LTA is defined such that it
+        maps coordinates from the ``dest volume`` to the ``src volume``.
+        Therefore, without inversion, the LTA matrix is appropiate
+        to move the information from ``src volume`` into the
+        ``dest volume``'s grid.
+
+        .. important ::
+
+            The ``moving`` and ``reference`` parameters are dismissed
+            because ``VOX2VOX`` LTAs are converted to ``RAS2RAS`` type
+            before returning the RAS+ matrix, using the ``dest`` and
+            ``src`` contained in the LTA. Both arguments are kept for
+            API compatibility.
+
+        Parameters
+        ----------
+        moving : dismissed
+            The spatial reference of moving images.
+        reference : dismissed
+            The spatial reference of moving images.
+
+        Returns
+        -------
+        matrix : :obj:`numpy.ndarray`
+            The RAS+ affine matrix corresponding to the LTA.
+
+        """
         self.set_type(1)
-        return self.structarr['m_L']
+        return np.linalg.inv(self.structarr['m_L'])
 
     def to_string(self):
         """Convert this transform to text."""
diff --git 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", )
 
     def __init__(self, matrix=None, reference=None):
         """
@@ -40,7 +40,7 @@ def __init__(self, matrix=None, reference=None):
 
         Examples
         --------
-        >>> xfm = Affine(reference=datadir / 'someones_anatomy.nii.gz')
+        >>> xfm = Affine(reference=datadir / "someones_anatomy.nii.gz")
         >>> xfm.matrix  # doctest: +NORMALIZE_WHITESPACE
         array([[1., 0., 0., 0.],
                [0., 1., 0., 0.],
@@ -61,13 +61,17 @@ def __init__(self, matrix=None, reference=None):
         if matrix is not None:
             matrix = np.array(matrix)
             if matrix.ndim != 2:
-                raise TypeError('Affine should be 2D.')
+                raise TypeError("Affine should be 2D.")
             elif matrix.shape[0] != matrix.shape[1]:
-                raise TypeError('Matrix is not square.')
+                raise TypeError("Matrix is not square.")
             self._matrix = matrix
 
-            if np.any(self._matrix[3, :] != (0, 0, 0, 1)):
-                raise ValueError("Matrix does not represent a valid transform.")
+            if not np.allclose(self._matrix[3, :], (0, 0, 0, 1)):
+                raise ValueError("""The last row of a homogeneus matrix \
+should be (0, 0, 0, 1), got %s.""" % self._matrix[3, :])
+
+            # Normalize last row
+            self._matrix[3, :] = (0, 0, 0, 1)
 
     def __eq__(self, other):
         """
@@ -83,7 +87,7 @@ def __eq__(self, other):
         """
         _eq = np.allclose(self.matrix, other.matrix, rtol=EQUALITY_TOL)
         if _eq and self._reference != other._reference:
-            warnings.warn('Affines are equal, but references do not match.')
+            warnings.warn("Affines are equal, but references do not match.")
         return _eq
 
     @property
@@ -125,16 +129,16 @@ def map(self, x, inverse=False):
 
     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)])
+        xform = x5_root.create_dataset("Transform", data=[self._matrix])
+        xform.attrs["Type"] = "affine"
+        x5_root.create_dataset("Inverse", data=[np.linalg.inv(self._matrix)])
 
         if self._reference:
-            self.reference._to_hdf5(x5_root.create_group('Reference'))
+            self.reference._to_hdf5(x5_root.create_group("Reference"))
 
-    def to_filename(self, filename, fmt='X5', moving=None):
+    def to_filename(self, filename, fmt="X5", moving=None):
         """Store the transform in BIDS-Transforms HDF5 file format (.x5)."""
-        if fmt.lower() in ['itk', 'ants', 'elastix']:
+        if fmt.lower() in ["itk", "ants", "elastix"]:
             itkobj = io.itk.ITKLinearTransform.from_ras(self.matrix)
             itkobj.to_filename(filename)
             return filename
@@ -145,45 +149,45 @@ def to_filename(self, filename, fmt='X5', moving=None):
         else:
             moving = self.reference
 
-        if fmt.lower() == 'afni':
+        if fmt.lower() == "afni":
             afniobj = io.afni.AFNILinearTransform.from_ras(
                 self.matrix, moving=moving, reference=self.reference)
             afniobj.to_filename(filename)
             return filename
 
-        if fmt.lower() == 'fsl':
+        if fmt.lower() == "fsl":
             fslobj = io.fsl.FSLLinearTransform.from_ras(
                 self.matrix, moving=moving, reference=self.reference
             )
             fslobj.to_filename(filename)
             return filename
 
-        if fmt.lower() == 'fs':
+        if fmt.lower() == "fs":
             # xform info
             lt = io.LinearTransform()
-            lt['sigma'] = 1.
-            lt['m_L'] = self.matrix
+            lt["sigma"] = 1.
+            lt["m_L"] = self.matrix
             # Just for reference, nitransforms does not write VOX2VOX
-            lt['src'] = io.VolumeGeometry.from_image(moving)
-            lt['dst'] = io.VolumeGeometry.from_image(self.reference)
+            lt["src"] = io.VolumeGeometry.from_image(moving)
+            lt["dst"] = io.VolumeGeometry.from_image(self.reference)
             # to make LTA file format
             lta = io.LinearTransformArray()
-            lta['type'] = 1  # RAS2RAS
-            lta['xforms'].append(lt)
+            lta["type"] = 1  # RAS2RAS
+            lta["xforms"].append(lt)
 
-            with open(filename, 'w') as f:
+            with open(filename, "w") as f:
                 f.write(lta.to_string())
             return filename
 
         raise NotImplementedError
 
     @classmethod
-    def from_filename(cls, filename, fmt='X5',
+    def from_filename(cls, filename, fmt="X5",
                       reference=None, moving=None):
         """Create an affine from a transform file."""
-        if fmt.lower() in ('itk', 'ants', 'elastix'):
+        if fmt.lower() in ("itk", "ants", "elastix"):
             _factory = io.itk.ITKLinearTransformArray
-        elif fmt.lower() in ('lta', 'fs'):
+        elif fmt.lower() in ("lta", "fs"):
             _factory = io.LinearTransformArray
         else:
             raise NotImplementedError
@@ -193,7 +197,7 @@ def from_filename(cls, filename, fmt='X5',
         if cls == Affine:
             if np.shape(matrix)[0] != 1:
                 raise TypeError(
-                    'Cannot load transform array "%s"' % filename)
+                    "Cannot load transform array '%s'" % filename)
             matrix = matrix[0]
         return cls(matrix, reference=reference)
 
@@ -297,9 +301,9 @@ def map(self, x, inverse=False):
             affine = np.linalg.inv(affine)
         return np.swapaxes(affine.dot(coords), 1, 2)
 
-    def to_filename(self, filename, fmt='X5', moving=None):
+    def to_filename(self, filename, fmt="X5", moving=None):
         """Store the transform in BIDS-Transforms HDF5 file format (.x5)."""
-        if fmt.lower() in ('itk', 'ants', 'elastix'):
+        if fmt.lower() in ("itk", "ants", "elastix"):
             itkobj = io.itk.ITKLinearTransformArray.from_ras(self.matrix)
             itkobj.to_filename(filename)
             return filename
@@ -310,41 +314,41 @@ def to_filename(self, filename, fmt='X5', moving=None):
         else:
             moving = self.reference
 
-        if fmt.lower() == 'afni':
+        if fmt.lower() == "afni":
             afniobj = io.afni.AFNILinearTransformArray.from_ras(
                 self.matrix, moving=moving, reference=self.reference)
             afniobj.to_filename(filename)
             return filename
 
-        if fmt.lower() == 'fsl':
+        if fmt.lower() == "fsl":
             fslobj = io.fsl.FSLLinearTransformArray.from_ras(
                 self.matrix, moving=moving, reference=self.reference
             )
             fslobj.to_filename(filename)
             return filename
 
-        if fmt.lower() in ('fs', 'lta'):
+        if fmt.lower() in ("fs", "lta"):
             # xform info
             # to make LTA file format
             lta = io.LinearTransformArray()
-            lta['type'] = 1  # RAS2RAS
+            lta["type"] = 1  # RAS2RAS
             for m in self.matrix:
                 lt = io.LinearTransform()
-                lt['sigma'] = 1.
-                lt['m_L'] = m
+                lt["sigma"] = 1.
+                lt["m_L"] = m
                 # Just for reference, nitransforms does not write VOX2VOX
-                lt['src'] = io.VolumeGeometry.from_image(moving)
-                lt['dst'] = io.VolumeGeometry.from_image(self.reference)
-                lta['xforms'].append(lt)
+                lt["src"] = io.VolumeGeometry.from_image(moving)
+                lt["dst"] = io.VolumeGeometry.from_image(self.reference)
+                lta["xforms"].append(lt)
 
-            with open(filename, 'w') as f:
+            with open(filename, "w") as f:
                 f.write(lta.to_string())
             return filename
 
         raise NotImplementedError
 
     def apply(self, spatialimage, reference=None,
-              order=3, mode='constant', cval=0.0, prefilter=True, output_dtype=None):
+              order=3, mode="constant", cval=0.0, prefilter=True, output_dtype=None):
         """
         Apply a transformation to an image, resampling on the reference spatial object.
 
@@ -359,11 +363,11 @@ def apply(self, spatialimage, reference=None,
         order : int, optional
             The order of the spline interpolation, default is 3.
             The order has to be in the range 0-5.
-        mode : {'constant', 'reflect', 'nearest', 'mirror', 'wrap'}, optional
+        mode : {"constant", "reflect", "nearest", "mirror", "wrap"}, optional
             Determines how the input image is extended when the resamplings overflows
-            a border. Default is 'constant'.
+            a border. Default is "constant".
         cval : float, optional
-            Constant value for ``mode='constant'``. Default is 0.0.
+            Constant value for ``mode="constant"``. Default is 0.0.
         prefilter: bool, optional
             Determines if the image's data array is prefiltered with
             a spline filter before interpolation. The default is ``True``,
@@ -436,17 +440,17 @@ def apply(self, spatialimage, reference=None,
         return resampled
 
 
-def load(filename, fmt='X5', reference=None, moving=None):
+def load(filename, fmt="X5", reference=None, moving=None):
     """
     Load a linear transform file.
 
     Examples
     --------
-    >>> xfm = load(datadir / 'affine-LAS.itk.tfm', fmt='itk')
+    >>> xfm = load(datadir / "affine-LAS.itk.tfm", fmt="itk")
     >>> isinstance(xfm, Affine)
     True
 
-    >>> xfm = load(datadir / 'itktflist.tfm', fmt='itk')
+    >>> xfm = load(datadir / "itktflist.tfm", fmt="itk")
     >>> isinstance(xfm, LinearTransformsMapping)
     True
 
diff --git a/nitransforms/tests/data/regressions/bbregister.fsl b/nitransforms/tests/data/regressions/bbregister.fsl
@@ -0,0 +1,4 @@
+   0.99970585   -0.00953967   -0.02228683  -15.19269657 
+  -0.00599674    0.79344094   -0.60861677   70.83446503 
+   0.02348929    0.60857153    0.79315072   15.39251423 
+   0.00000000    0.00000000    0.00000000    1.00000000 
diff --git a/nitransforms/tests/data/regressions/bbregister.lta b/nitransforms/tests/data/regressions/bbregister.lta
@@ -0,0 +1,29 @@
+type      = 0 # LINEAR_VOX_TO_VOX
+nxforms   = 1
+mean      = 0.0000 0.0000 0.0000
+sigma     = 1.0000
+1 4 4
+-3.124080896377563e+00 2.981145866215229e-02 8.914728462696075e-02 1.741926879882812e+02 
+-1.405486371368170e-02 1.859627604484558e+00 -1.825850725173950e+00 5.312585067749023e+01 
+5.505303665995598e-02 1.426339745521545e+00 2.379452466964722e+00 1.154438781738281e+01 
+0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 9.999998807907104e-01 
+src volume info
+valid = 1  # volume info valid
+filename = /home/oesteban/tmp/fmriprep-ds005/fprep-work/fmriprep_wf/single_subject_01_wf/func_preproc_task_mixedgamblestask_run_01_wf/bold_reg_wf/bbreg_wf/bbregister/uni_xform_masked.nii.gz
+volume = 64 64 34
+voxelsize = 3.125000000000000e+00 3.125000000000000e+00 4.000000000000000e+00
+xras   = -1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
+yras   = 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00
+zras   = 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00
+cras   = 1.000000000000000e+00 2.800000000000000e+01 -3.100000000000000e+01
+dst volume info
+valid = 1  # volume info valid
+filename = /oak/stanford/groups/russpold/data/openfmri/ds000005/sub-01/anat/sub-01_T1w.nii.gz
+volume = 160 192 192
+voxelsize = 1.000000000000000e+00 1.333333015441895e+00 1.333333015441895e+00
+xras   = 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
+yras   = 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00
+zras   = 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00
+cras   = -1.000000000000000e+00 -5.000030517578125e+00 -1.000030517578125e+00
+subject sub-01
+fscale 0.100000
diff --git a/nitransforms/tests/data/regressions/bbregister.tfm b/nitransforms/tests/data/regressions/bbregister.tfm
@@ -0,0 +1,5 @@
+#Insight Transform File V1.0
+#Transform 0
+Transform: MatrixOffsetTransformBase_double_3_3
+Parameters: 0.999706355364118 0.0059967416827344956 0.023489316177446647 0.009539669218827523 0.7934419194347854 -0.6085721629337524 -0.02228683774477952 0.6086173680484374 0.7931513449655613 -5.5388974668589555 -45.57406723091313 -48.80406464758718
+FixedParameters: 0 0 0
diff --git a/nitransforms/tests/data/regressions/robust_register.lta b/nitransforms/tests/data/regressions/robust_register.lta
@@ -1,30 +1,28 @@
-# transform file /home/oesteban/tmp/fmriprep-ds005/fprep-work/fmriprep_wf/single_subject_01_wf/anat_preproc_wf/surface_recon_wf/fsnative2t1w_xfm/T1_robustreg.lta
-# created by oesteban on Sat Mar 14 19:28:37 2020
-
 type      = 1 # LINEAR_RAS_TO_RAS
 nxforms   = 1
 mean      = 129.0000 157.0000 132.0000
 sigma     = 10000.0000
 1 4 4
-9.999999403953552e-01 -1.698292035143822e-04 1.542967074783519e-04 -1.678466796875000e-04 
-1.698438863968477e-04 9.999999403953552e-01 -9.513227996649221e-05 -1.318359375000000e-02 
--1.542805403005332e-04 9.515848068986088e-05 9.999999403953552e-01 -6.271362304687500e-03 
-0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 9.999999403953552e-01 
+1.000000000000000e+00 1.698439009487629e-04 -1.542805694043636e-04 1.691182987997308e-04 
+-1.698292180662975e-04 1.000000000000000e+00 9.515849524177611e-05 1.318416278809309e-02 
+1.542967220302671e-04 -9.513228724244982e-05 1.000000000000000e+00 6.270134821534157e-03 
+0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 
 src volume info
 valid = 1  # volume info valid
-filename = /oak/stanford/groups/russpold/data/openfmri/derivatives/ds000005/freesurfer-6.0.1/sub-01/mri/T1.mgz
-volume = 256 256 256
-voxelsize = 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00
-xras   = -9.999999403953552e-01 0.000000000000000e+00 0.000000000000000e+00
-yras   = 0.000000000000000e+00 0.000000000000000e+00 -9.999999403953552e-01
-zras   = 0.000000000000000e+00 9.999999403953552e-01 0.000000000000000e+00
-cras   = -9.999847412109375e-01 -5.000015258789062e+00 -1.000038146972656e+00
-dst volume info
-valid = 1  # volume info valid
 filename = /oak/stanford/groups/russpold/data/openfmri/ds000005/sub-01/anat/sub-01_T1w.nii.gz
 volume = 160 192 192
 voxelsize = 1.000000000000000e+00 1.333333015441895e+00 1.333333015441895e+00
 xras   = 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
 yras   = 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00
 zras   = 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00
 cras   = -1.000000000000000e+00 -5.000030517578125e+00 -1.000030517578125e+00
+dst volume info
+valid = 1  # volume info valid
+filename = /oak/stanford/groups/russpold/data/openfmri/derivatives/ds000005/freesurfer-6.0.1/sub-01/mri/T1.mgz
+volume = 256 256 256
+voxelsize = 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00
+xras   = -9.999999403953552e-01 0.000000000000000e+00 0.000000000000000e+00
+yras   = 0.000000000000000e+00 0.000000000000000e+00 -9.999999403953552e-01
+zras   = 0.000000000000000e+00 9.999999403953552e-01 0.000000000000000e+00
+cras   = -9.999847412109375e-01 -5.000015258789062e+00 -1.000038146972656e+00
+fscale 0.100000
diff --git a/nitransforms/tests/data/regressions/robust_register.tfm b/nitransforms/tests/data/regressions/robust_register.tfm
@@ -0,0 +1,5 @@
+#Insight Transform File V1.0
+#Transform 0
+Transform: AffineTransform_double_3_3
+Parameters: 1 -0.00016982921806629747 -0.00015429673658218235 0.00016984390094876289 1 9.5132294518407434e-05 0.00015428056940436363 -9.5158495241776109e-05 1 0.00016784669423941523 0.013183594681322575 -0.0062713632360100746
+FixedParameters: 0 0 0
diff --git a/nitransforms/tests/test_conversions.py b/nitransforms/tests/test_conversions.py
@@ -0,0 +1,21 @@
+# 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
+from .. import linear as _l
+
+
+def test_conversions0(data_path):
+    """Check conversions between formats."""
+    lta = _l.load(data_path / "regressions" / "robust_register.lta", fmt="lta")
+    itk = _l.load(data_path / "regressions" / "robust_register.tfm", fmt="itk")
+
+    assert np.allclose(lta.matrix, itk.matrix)
+
+
+def test_conversions1(data_path):
+    """Check conversions between formats."""
+    lta = _l.load(data_path / "regressions" / "bbregister.lta", fmt="lta")
+    itk = _l.load(data_path / "regressions" / "bbregister.tfm", fmt="itk")
+
+    assert np.allclose(lta.matrix, itk.matrix)
