ENH: AFNI - Enable conversion to RAS+ of affine transforms

nitransforms/io/afni.py

@@ -1,7 +1,8 @@
 """Read/write AFNI's transforms."""
 from math import pi
+import warnings
 import numpy as np
-from nibabel.affines import obliquity, voxel_sizes
+from nibabel.affines import obliquity, voxel_sizes, from_matvec
 
 from ..patched import shape_zoom_affine
 from .base import (
@@ -13,6 +14,8 @@
 
 LPS = np.diag([-1, -1, 1, 1])
 OBLIQUITY_THRESHOLD_DEG = 0.01
+B = np.ones((2, 2))
+AFNI_SIGNS = np.block([[B, -1.0 * B], [-1.0 * B, B]])
 
 
 class AFNILinearTransform(LinearParameters):
@@ -23,6 +26,12 @@ def __str__(self):
         param = self.structarr['parameters']
         return '\t'.join(['%g' % p for p in param[:3, :].reshape(-1)])
 
+    def to_ras(self, moving=None, reference=None):
+        """Convert to RAS+ coordinate system."""
+        afni = self.structarr['parameters'].copy()
+        # swapaxes is necessary, as axis 0 encodes series of transforms
+        return _afni2ras(afni, moving=moving, reference=reference)
+
     def to_string(self, banner=True):
         """Convert to a string directly writeable to file."""
         string = '%s\n' % self
@@ -33,29 +42,9 @@ def to_string(self, banner=True):
 
     @classmethod
     def from_ras(cls, ras, moving=None, reference=None):
-        """Create an ITK affine from a nitransform's RAS+ matrix."""
-        ras = ras.copy()
-        pre = LPS.copy()
-        post = LPS.copy()
-        if _is_oblique(reference.affine):
-            print('Reference affine axes are oblique.')
-            M = reference.affine
-            A = shape_zoom_affine(reference.shape,
-                                  voxel_sizes(M), x_flip=False, y_flip=False)
-            pre = M.dot(np.linalg.inv(A)).dot(LPS)
-
-        if _is_oblique(moving.affine):
-            print('Moving affine axes are oblique.')
-            M2 = moving.affine
-            A2 = shape_zoom_affine(moving.shape,
-                                   voxel_sizes(M2), x_flip=True, y_flip=True)
-            post = A2.dot(np.linalg.inv(M2))
-
-        # swapaxes is necessary, as axis 0 encodes series of transforms
-        parameters = np.swapaxes(post.dot(ras.dot(pre)), 0, 1)
-
+        """Create an AFNI affine from a nitransform's RAS+ matrix."""
         tf = cls()
-        tf.structarr['parameters'] = parameters.T
+        tf.structarr['parameters'] = _ras2afni(ras, moving, reference)
         return tf
 
     @classmethod
@@ -97,7 +86,7 @@ def to_string(self):
                 l.strip()
                 for l in xfm.to_string(banner=(i == 0)).splitlines()
                 if l.strip()]
-            strings += lines
+            strings += lines + ['']
         return '\n'.join(strings)
 
     @classmethod
@@ -150,3 +139,100 @@ def from_image(cls, imgobj):
 
 def _is_oblique(affine, thres=OBLIQUITY_THRESHOLD_DEG):
     return (obliquity(affine).min() * 180 / pi) > thres
+
+
+def _afni_warpdrive_for(oblique, plumb, offset=True, inv=False):
+    """
+    Calculate AFNI's ``WARPDRIVE_MATVEC_FOR_000000`` (de)obliquing affine.
+
+    Parameters
+    ----------
+    oblique : 4x4 numpy.array
+        affine that is not aligned to the cardinal axes.
+    plumb : 4x4 numpy.array
+        corresponding affine that is aligned to the cardinal axes.
+
+
+    Returns
+    -------
+    plumb_to_oblique : 4x4 numpy.array
+        the matrix that pre-pended to the plumb affine rotates it
+        to be oblique.
+
+    """
+    R = np.linalg.inv(plumb[:3, :3]).dot(oblique[:3, :3])
+    origin = oblique[:3, 3] - R.dot(oblique[:3, 3])
+    if offset is False:
+        origin = np.zeros(3)
+    matvec_inv = from_matvec(R, origin)  # * AFNI_SIGNS
+    if not inv:
+        return np.linalg.inv(matvec_inv)
+    return matvec_inv
+
+
+def _ras2afni(ras, moving=None, reference=None):
+    """
+    Convert from RAS+ to AFNI matrix.
+
+    inverse : bool
+        if ``False`` (default), return the matrix to rotate plumb
+        onto oblique (AFNI's ``WARPDRIVE_MATVEC_INV_000000``);
+        if ``True``, return the matrix to rotate oblique onto
+        plumb (AFNI's ``WARPDRIVE_MATVEC_FOR_000000``).
+
+    """
+    ras = ras.copy()
+    pre = np.eye(4)
+    post = np.eye(4)
+    if reference is not None and _is_oblique(reference.affine):
+        warnings.warn('Reference affine axes are oblique.')
+        M = reference.affine
+        plumb = shape_zoom_affine(reference.shape, voxel_sizes(M))
+        # Prepend the MATVEC_INV - AFNI will append MATVEC_FOR
+        pre = _afni_warpdrive_for(M, plumb, offset=False, inv=False)
+
+    if moving is not None and _is_oblique(moving.affine):
+        warnings.warn('Moving affine axes are oblique.')
+        M = moving.affine
+        plumb = shape_zoom_affine(moving.shape, voxel_sizes(M))
+        # Append the MATVEC_FOR - AFNI will append MATVEC_INV
+        post = _afni_warpdrive_for(M, plumb, offset=False, inv=True)
+
+    afni_ras = np.swapaxes(post.dot(ras.dot(pre)), 0, 1).T
+
+    # Combine oblique/deoblique matrices into RAS+ matrix
+    return np.swapaxes(LPS.dot(afni_ras.dot(LPS)), 0, 1).T
+
+
+def _afni2ras(afni, moving=None, reference=None):
+    """
+    Convert from AFNI to RAS+.
+
+    inverse : bool
+        if ``False`` (default), return the matrix to rotate plumb
+        onto oblique (AFNI's ``WARPDRIVE_MATVEC_INV_000000``);
+        if ``True``, return the matrix to rotate oblique onto
+        plumb (AFNI's ``WARPDRIVE_MATVEC_FOR_000000``).
+
+    """
+    afni = afni.copy()
+    pre = np.eye(4)
+    post = np.eye(4)
+    if reference is not None and _is_oblique(reference.affine):
+        warnings.warn('Reference affine axes are oblique.')
+        M = reference.affine
+        plumb = shape_zoom_affine(reference.shape, voxel_sizes(M))
+        # Append the MATVEC_FOR - AFNI would add it implicitly
+        pre = _afni_warpdrive_for(M, plumb, offset=False)
+
+    if moving is not None and _is_oblique(moving.affine):
+        warnings.warn('Moving affine axes are oblique.')
+        M = moving.affine
+        plumb = shape_zoom_affine(moving.shape, voxel_sizes(M))
+        # Prepend the MATVEC_INV - AFNI will add it implicitly
+        post = _afni_warpdrive_for(M, plumb, offset=False, inv=False)
+
+    afni_ras = np.swapaxes(post.dot(afni.dot(pre)), 0, 1).T
+
+    # Combine oblique/deoblique matrices into RAS+ matrix
+    return np.swapaxes(LPS.dot(afni_ras.dot(LPS)), 0, 1).T

nitransforms/tests/data/affine-oblique-catmatvec.afni

@@ -0,0 +1 @@
+     0.999998   -0.00206779  -0.000680893      -4.02301    0.00181873      0.621609      0.783326       1.51751    -0.0011965     -0.783325      0.621611      -9.66241

nitransforms/tests/data/affine-oblique-catmatvec2.afni

@@ -0,0 +1,2 @@
+     0.999998   -0.00206779  -0.000680893      -4.02301    0.00181873      0.621609      0.783326       1.51751    -0.0011965     -0.783325      0.621611      -9.66241
+     0.999998   -0.00206779  -0.000680895      -4.00307    0.00181873      0.621609      0.783326      -1.89802    -0.0011965     -0.783325      0.621611      -1.17171

nitransforms/tests/data/affine-oblique.afni

@@ -0,0 +1,2 @@
+# 3dvolreg matrices (DICOM-to-DICOM, row-by-row):
+0.999999	-0.000255644	0.00149077	4.00707	-0.00100885	0.62161	0.783326	1.72715	-0.00112693	-0.783327	0.621609	-9.6758

nitransforms/tests/test_io.py

@@ -135,15 +135,16 @@ def test_Linear_common(tmpdir, data_path, sw, image_orientation,
         xfm = factory.from_fileobj(f)
 
     # Test to_string
-    assert text == xfm.to_string()
+    assert text.strip() == xfm.to_string().strip()
 
     xfm.to_filename(fname)
     assert filecmp.cmp(fname, str((data_path / fname).resolve()))
 
     # Test from_ras
     RAS = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0])
     xfm = factory.from_ras(RAS, reference=reference, moving=moving)
-    assert np.allclose(xfm.to_ras(reference=reference, moving=moving), RAS)
+    recovered_ras = xfm.to_ras(reference=reference, moving=moving)
+    assert np.allclose(recovered_ras, RAS)
 
 
 @pytest.mark.parametrize('image_orientation', [

nitransforms/tests/test_linear.py

@@ -31,7 +31,7 @@
 }
 
 
-@pytest.mark.xfail(reason="Not fully implemented")
+@pytest.mark.xfail(raises=(FileNotFoundError, NotImplementedError))
 @pytest.mark.parametrize('image_orientation', ['RAS', 'LAS', 'LPS', 'oblique'])
 @pytest.mark.parametrize('sw_tool', ['itk', 'fsl', 'afni', 'fs'])
 def test_linear_save(tmpdir, data_path, get_testdata, image_orientation, sw_tool):
@@ -57,22 +57,28 @@ def test_linear_save(tmpdir, data_path, get_testdata, image_orientation, sw_tool
     assert_affines_by_filename(xfm_fname1, xfm_fname2)
 
 
-@pytest.mark.parametrize('image_orientation', [
-    'RAS', 'LAS', 'LPS',  # 'oblique',
-])
+@pytest.mark.parametrize('image_orientation', ['RAS', 'LAS', 'LPS', 'oblique'])
 @pytest.mark.parametrize('sw_tool', ['itk', 'fsl', 'afni'])
+@pytest.mark.parametrize('parameters', [
+    ((0.0, 0.0, 0.0), (0.0, 0.0, 0.0)),
+    ((0.0, 0.0, 0.0), (4.0, 2.0, -1.0)),
+    ((0.9, 0.001, 0.001), (0.0, 0.0, 0.0)),
+    ((0.9, 0.001, 0.001), (4.0, 2.0, -1.0)),
+])
 def test_apply_linear_transform(
         tmpdir,
         get_testdata,
         image_orientation,
-        sw_tool
+        sw_tool,
+        parameters
 ):
     """Check implementation of exporting affines to formats."""
     tmpdir.chdir()
 
     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])
+    T = from_matvec(euler2mat(*parameters[0]), parameters[1])
+
     xfm = ntl.Affine(T)
     xfm.reference = img
 
@@ -103,6 +109,50 @@ 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
 
+@pytest.mark.parametrize('cat', ['catmatvec', 'catmatvec2'])
+def test_afni_oblique(
+        tmpdir,
+        data_path,
+        get_testdata,
+        cat,
+):
+    """Check implementation of exporting affines to formats."""
+    tmpdir.chdir()
+
+    xfm_file0 = str(data_path / 'affine-RAS.afni')
+    xfm_file1 = str(data_path / 'affine-oblique-{}.afni'.format(cat))
+
+    get_testdata['RAS'].to_filename('RAS.nii.gz')
+    get_testdata['oblique'].to_filename('oblique.nii.gz')
+
+    cmd = """\
+3dAllineate -base {reference} -input {moving} -master {master} \
+-prefix resampled.nii.gz -1Dmatrix_apply {transform} -final NN""".format(
+        transform=os.path.abspath(xfm_file0),
+        reference=os.path.abspath('oblique.nii.gz'),
+        moving=os.path.abspath('oblique.nii.gz'),
+        master=os.path.abspath('RAS.nii.gz'))
+
+    # skip test if command is not available on host
+    exe = cmd.split(" ", 1)[0]
+    if not shutil.which(exe):
+        pytest.skip("Command {} not found on host".format(exe))
+
+    exit_code = check_call([cmd], shell=True)
+    assert exit_code == 0
+    shutil.move('resampled.nii.gz', 'resampled0.nii.gz')
+    sw_moved = nb.load('resampled0.nii.gz')
+
+    cmd = APPLY_LINEAR_CMD['afni'](
+        transform=os.path.abspath(xfm_file1),
+        reference=os.path.abspath('RAS.nii.gz'),
+        moving=os.path.abspath('RAS.nii.gz'))
+    check_call([cmd], shell=True)
+    nt_moved = nb.load('resampled.nii.gz')
+    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, data_path):
     """Test affine's operations."""