fix: finishing up this PR

With the generous support from Paul Taylor <@mrneont> and his answers in
afni/afni#353, I have managed to get all tests to PASS.

I am not resolving the problem of oblique datasets for displacements
fields this time.

Resolves: #45.
Resolves: #150.
Continues: #157.

Author: oesteban

nitransforms/io/afni.py

@@ -11,6 +11,7 @@
     DisplacementsField,
     LinearParameters,
     TransformFileError,
+    _ensure_image,
 )
 
 LPS = np.diag([-1, -1, 1, 1])
@@ -38,6 +39,15 @@ def to_string(self, banner=True):
     def from_ras(cls, ras, moving=None, reference=None):
         """Create an AFNI affine from a nitransform's RAS+ matrix."""
         # swapaxes is necessary, as axis 0 encodes series of transforms
+
+        reference = _ensure_image(reference)
+        if reference is not None and _is_oblique(reference.affine):
+            ras = ras @ _cardinal_rotation(reference.affine, False)
+
+        moving = _ensure_image(moving)
+        if moving is not None and _is_oblique(moving.affine):
+            ras = _cardinal_rotation(moving.affine, True) @ ras
+
         parameters = np.swapaxes(LPS @ ras @ LPS, 0, 1)
 
         tf = cls()
@@ -71,7 +81,16 @@ def from_string(cls, string):
     def to_ras(self, moving=None, reference=None):
         """Return a nitransforms internal RAS+ matrix."""
         # swapaxes is necessary, as axis 0 encodes series of transforms
-        return LPS @ np.swapaxes(self.structarr["parameters"].T, 0, 1) @ LPS
+        retval = LPS @ np.swapaxes(self.structarr["parameters"].T, 0, 1) @ LPS
+        reference = _ensure_image(reference)
+        if reference is not None and _is_oblique(reference.affine):
+            retval = retval @ _cardinal_rotation(reference.affine, True)
+
+        moving = _ensure_image(moving)
+        if moving is not None and _is_oblique(moving.affine):
+            retval = _cardinal_rotation(moving.affine, False) @ retval
+
+        return retval
 
 
 class AFNILinearTransformArray(BaseLinearTransformList):
@@ -244,6 +263,28 @@ def _dicom_real_to_card(oblique):
     return retval
 
 
+def _cardinal_rotation(oblique, real_to_card=True):
+    """
+    Calculate the rotation matrix to undo AFNI's deoblique operation.
+
+    Parameters
+    ----------
+    oblique : 4x4 numpy.array
+        affine that may not be aligned to the cardinal axes ("IJK_DICOM_REAL" for AFNI).
+
+    Returns
+    -------
+    plumb : 4x4 numpy.array
+        affine aligned to the cardinal axes ("IJK_DICOM_CARD" for AFNI).
+
+    """
+    card = _dicom_real_to_card(oblique)
+    return (
+        card @ np.linalg.inv(oblique) if real_to_card else
+        oblique @ np.linalg.inv(card)
+    )
+
+
 def _afni_warpdrive(oblique, forward=True, ras=False):
     """
     Calculate AFNI's ``WARPDRIVE_MATVEC_FOR_000000`` (de)obliquing affine.

nitransforms/tests/test_io.py

@@ -471,6 +471,30 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y,
     card_aff = afni._dicom_real_to_card(img.affine)
     assert np.allclose(card_aff, nb.load("deob_3drefit.nii.gz").affine)
 
+    # Check that nitransforms can emulate 3drefit -deoblique
+    nt3drefit = Affine(
+        afni._cardinal_rotation(img.affine, False),
+        reference="deob_3drefit.nii.gz",
+    ).apply("orig.nii.gz")
+
+    diff = (
+        np.asanyarray(img.dataobj, dtype="uint8")
+        - np.asanyarray(nt3drefit.dataobj, dtype="uint8")
+    )
+    assert np.sqrt((diff[10:-10, 10:-10, 10:-10] ** 2).mean()) < 0.1
+
+    # Check that nitransforms can revert 3drefit -deoblique
+    nt_undo3drefit = Affine(
+        afni._cardinal_rotation(img.affine, True),
+        reference="orig.nii.gz",
+    ).apply("deob_3drefit.nii.gz")
+
+    diff = (
+        np.asanyarray(img.dataobj, dtype="uint8")
+        - np.asanyarray(nt_undo3drefit.dataobj, dtype="uint8")
+    )
+    assert np.sqrt((diff[10:-10, 10:-10, 10:-10] ** 2).mean()) < 0.1
+
     # Check the target grid by 3dWarp and the affine & size interpolated by NiTransforms
     cmd = f"3dWarp -verb -deoblique -NN -prefix {tmpdir}/deob.nii.gz {tmpdir}/orig.nii.gz"
     assert check_call([cmd], shell=True) == 0

nitransforms/tests/test_linear.py

@@ -149,7 +149,14 @@ def test_linear_save(tmpdir, data_path, get_testdata, image_orientation, sw_tool
         # Account for the fact that FS defines LTA transforms reversed
         T = np.linalg.inv(T)
 
-    xfm = nitl.Affine(T)
+    xfm = (
+        nitl.Affine(T) if (sw_tool, image_orientation) != ("afni", "oblique") else
+        # AFNI is special when moving or reference are oblique - let io do the magic
+        nitl.Affine(nitl.io.afni.AFNILinearTransform.from_ras(T).to_ras(
+            reference=img,
+            moving=img,
+        ))
+    )
     xfm.reference = img
 
     ext = ""
@@ -190,7 +197,15 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient
 
     # Write out transform file (software-dependent)
     xfm_fname = "M.%s%s" % (sw_tool, ext)
-    xfm.to_filename(xfm_fname, fmt=sw_tool)
+    # Change reference dataset for AFNI & oblique
+    if (sw_tool, image_orientation) == ("afni", "oblique"):
+        nitl.io.afni.AFNILinearTransform.from_ras(
+            T,
+            moving=img,
+            reference=img,
+        ).to_filename(xfm_fname)
+    else:
+        xfm.to_filename(xfm_fname, fmt=sw_tool)
 
     cmd = APPLY_LINEAR_CMD[sw_tool](
         transform=os.path.abspath(xfm_fname),
@@ -209,19 +224,11 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient
     assert exit_code == 0
     sw_moved_mask = nb.load("resampled_brainmask.nii.gz")
 
-    # Change reference dataset for AFNI & oblique
-    if (sw_tool, image_orientation) == ("afni", "oblique"):
-        xfm.reference = "resampled_brainmask.nii.gz"
-
     nt_moved_mask = xfm.apply(msk, order=0)
     nt_moved_mask.set_data_dtype(msk.get_data_dtype())
-    nt_moved_mask.to_filename("nt_resampled_brainmask.nii.gz")
+    nt_moved_mask.to_filename("ntmask.nii.gz")
     diff = np.asanyarray(sw_moved_mask.dataobj) - np.asanyarray(nt_moved_mask.dataobj)
 
-    nt_moved_mask.__class__(
-        diff, sw_moved_mask.affine, sw_moved_mask.header
-    ).to_filename("diff.nii.gz")
-
     assert np.sqrt((diff ** 2).mean()) < RMSE_TOL
     brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool)
 

nitransforms/tests/test_nonlinear.py

@@ -112,6 +112,9 @@ def test_displacements_field1(
     axis,
 ):
     """Check a translation-only field on one or more axes, different image orientations."""
+    if (image_orientation, sw_tool) == ("oblique", "afni"):
+        pytest.skip("AFNI obliques are not yet implemented for displacements fields")
+
     os.chdir(str(tmp_path))
     nii = get_testdata[image_orientation]
     msk = get_testmask[image_orientation]