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

nitransforms/io/afni.py

@@ -30,8 +30,7 @@ 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 _ras2afni(afni, moving=moving, reference=reference,
-                         inverse=True)
+        return _afni2ras(afni, moving=moving, reference=reference)
 
     def to_string(self, banner=True):
         """Convert to a string directly writeable to file."""
@@ -142,9 +141,9 @@ def _is_oblique(affine, thres=OBLIQUITY_THRESHOLD_DEG):
     return (obliquity(affine).min() * 180 / pi) > thres
 
 
-def _afni_to_oblique(oblique, plumb, inverse=False):
+def _afni_warpdrive_for(oblique, plumb, offset=True, inv=False):
     """
-    Calculate AFNI's matrices to (de)oblique affines.
+    Calculate AFNI's ``WARPDRIVE_MATVEC_FOR_000000`` (de)obliquing affine.
 
     Parameters
     ----------
@@ -163,37 +162,77 @@ def _afni_to_oblique(oblique, plumb, inverse=False):
     """
     R = np.linalg.inv(plumb[:3, :3]).dot(oblique[:3, :3])
     origin = oblique[:3, 3] - R.dot(oblique[:3, 3])
-    return from_matvec(R, origin) * AFNI_SIGNS
+    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, inverse=False):
+def _ras2afni(ras, moving=None, reference=None):
     """
-    Convert from AFNI to RAS+.
+    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 = LPS.copy()
-    post = LPS.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))
-        pre = pre.dot(_afni_to_oblique(M, plumb))
-        if inverse is True:
-            pre = np.linalg.inv(pre)
+        # Prepend the MATVEC_INV - AFNI will append MATVEC_FOR
+        pre = _afni_warpdrive_for(M, plumb, offset=False, inv=True)
 
     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))
-        post = post.dot(_afni_to_oblique(M, plumb))
-        if inverse is False:
-            post = np.linalg.inv(post)
+        # Append the MATVEC_FOR - AFNI will append MATVEC_INV
+        post = _afni_warpdrive_for(M, plumb, offset=False)
+
+    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 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``).
+
+    """
+    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(post.dot(ras.dot(pre)), 0, 1).T
+    return np.swapaxes(LPS.dot(afni_ras.dot(LPS)), 0, 1).T

nitransforms/tests/test_linear.py

@@ -57,9 +57,7 @@ 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'])
 def test_apply_linear_transform(
         tmpdir,