fix: implement AFNI's dicom real to dicom card and test

Author: oesteban

nitransforms/io/afni.py

@@ -214,16 +214,44 @@ def _afni_deobliqued_grid(oblique, shape):
     return plumb, nshape
 
 
-def _afni_warpdrive(oblique, cardinal, forward=True, ras=False):
+def _dicom_real_to_card(oblique):
+    """
+    Calculate the corresponding "DICOM cardinal" for "DICOM real" (AFNI jargon).
+
+    Implements the internal "deobliquing" operation of ``3drefit`` and other tools, which
+    just *drop* the obliquity from the input affine.
+
+    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).
+
+    """
+    # Origin is kept from input
+    retval = np.eye(4)
+    retval[:3, 3] = oblique[:3, 3]
+
+    # Calculate director cosines and project to closest canonical
+    cosines = oblique[:3, :3] / np.abs(oblique[:3, :3]).max(0)
+    cosines[np.abs(cosines) < 1.0] = 0
+    # Once director cosines are calculated, scale by voxel sizes
+    retval[:3, :3] = np.round(voxel_sizes(oblique), decimals=4) * cosines
+    return retval
+
+
+def _afni_warpdrive(oblique, forward=True, ras=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.
-    cardinal : 4x4 numpy.array
-        corresponding affine that is aligned to the cardinal axes.
     forward : :obj:`bool`
         Returns the forward transformation if True, i.e.,
         the matrix to convert an oblique affine into an AFNI's plumb (if ``True``)
@@ -237,18 +265,8 @@ def _afni_warpdrive(oblique, cardinal, forward=True, ras=False):
         AFNI's *warpdrive* forward or inverse matrix.
 
     """
-    # Rotate the oblique affine to align with imaging axes
-    # Calculate director cosines and project to closest canonical
-
-    # plumb_r = oblique[:3, :3] / np.abs(oblique[:3, :3]).max(0)
-    # plumb_r[np.abs(plumb_r) < 1.0] = 0
-    # # Scale by min voxel size (AFNI's default)
-    # plumb_r *= vs.min()
-    # plumb = np.eye(4)
-    # plumb[:3, :3] = plumb_r
-
     ijk_to_dicom_real = np.diag(LPS) * oblique
-    ijk_to_dicom = cardinal
+    ijk_to_dicom = _dicom_real_to_card(oblique)
     R = np.linalg.inv(ijk_to_dicom) @ ijk_to_dicom_real
     return np.linalg.inv(R) if forward else R
 

nitransforms/tests/test_io.py

@@ -12,6 +12,7 @@
 from nibabel.eulerangles import euler2mat
 from nibabel.affines import from_matvec
 from scipy.io import loadmat
+from nitransforms.linear import Affine
 from ..io import (
     afni,
     fsl,
@@ -446,59 +447,36 @@ def test_regressions(file_type, test_file, data_path):
 @pytest.mark.parametrize("dir_y", (-1, 1))
 @pytest.mark.parametrize("dir_z", (1, -1))
 @pytest.mark.parametrize("swapaxes", [
-    None,  # (0, 1), (1, 2), (0, 2),
+    None, (0, 1), (1, 2), (0, 2),
 ])
 def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, dir_z):
     tmpdir.chdir()
-    img = nb.load(testdata_path / "someones_anatomy.nii.gz")
-    shape = np.array(img.shape[:3])
-    hdr = img.header.copy()
-    aff = img.affine.copy()
-    data = np.asanyarray(img.dataobj, dtype="uint8")
-
-    directions = (dir_x, dir_y, dir_z)
-    if directions != (1, 1, 1):
-        last_ijk = np.hstack((shape - 1, 1.0))
-        last_xyz = aff @ last_ijk
-        aff = np.diag((dir_x, dir_y, dir_z, 1)) @ aff
-
-        for ax in range(3):
-            if (directions[ax] == -1):
-                aff[ax, 3] = last_xyz[ax]
-                data = np.flip(data, ax)
-
-        hdr.set_qform(aff, code=1)
-        hdr.set_sform(aff, code=1)
-        img.__class__(data, aff, hdr).to_filename("flips.nii.gz")
-
-    if swapaxes is not None:
-        data = np.swapaxes(data, swapaxes[0], swapaxes[1])
-        aff[list(reversed(swapaxes)), :] = aff[(swapaxes), :]
-
-        hdr.set_qform(aff, code=1)
-        hdr.set_sform(aff, code=1)
-        img.__class__(data, aff, hdr).to_filename("swaps.nii.gz")
+    img, R = _generate_reoriented(
+        testdata_path / "someones_anatomy.nii.gz",
+        (dir_x, dir_y, dir_z),
+        swapaxes,
+        parameters
+    )
+    img.to_filename("orig.nii.gz")
 
-    R = from_matvec(euler2mat(**parameters), [0.0, 0.0, 0.0])
+    # Run AFNI's 3drefit -deoblique
+    if not shutil.which("3drefit"):
+        pytest.skip("Command 3drefit not found on host")
 
-    # img_center = aff @ np.hstack((shape * 0.5, 1.0))
-    # R[:3, 3] += (img_center - (R @ aff @ np.hstack((shape * 0.5, 1.0))))[:3]
-    newaff = R @ aff
-    hdr.set_qform(newaff, code=1)
-    hdr.set_sform(newaff, code=1)
-    img = img.__class__(data, newaff, hdr)
-    img.to_filename("oblique.nii.gz")
+    shutil.copy(f"{tmpdir}/orig.nii.gz", f"{tmpdir}/deob_3drefit.nii.gz")
+    cmd = f"3drefit -deoblique {tmpdir}/deob_3drefit.nii.gz"
+    assert check_call([cmd], shell=True) == 0
 
-    # Run AFNI's 3dDeoblique
-    if not shutil.which("3dWarp"):
-        pytest.skip("Command 3dWarp not found on host")
+    # Check that nitransforms can make out the deoblique affine:
+    card_aff = afni._dicom_real_to_card(img.affine)
+    assert np.allclose(card_aff, nb.load("deob_3drefit.nii.gz").affine)
 
-    cmd = f"3dWarp -verb -deoblique -NN -prefix {tmpdir}/deob.nii.gz {tmpdir}/oblique.nii.gz"
+    # 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
 
-    # Check the target grid by 3dWarp and the affine & size interpolated by NiTransforms
-    deobaff, deobshape = afni._afni_deobliqued_grid(newaff, shape)
     deobnii = nb.load("deob.nii.gz")
+    deobaff, deobshape = afni._afni_deobliqued_grid(img.affine, img.shape)
 
     assert np.all(deobshape == deobnii.shape[:3])
     assert np.allclose(deobaff, deobnii.affine)
@@ -509,13 +487,21 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y,
         deobaff,
         deobnii.header
     )).apply(img, order=0)
-    ntdeobnii.to_filename("ntdeob.nii.gz")
+
+    # Generate an internal box to exclude border effects
+    box = np.zeros(img.shape, dtype="uint8")
+    box[10:-10, 10:-10, 10:-10] = 1
+    ntdeobmask = Affine(np.eye(4), reference=ntdeobnii).apply(
+        nb.Nifti1Image(box, img.affine, img.header),
+        order=0,
+    )
+    mask = np.asanyarray(ntdeobmask.dataobj, dtype=bool)
+
     diff = (
         np.asanyarray(deobnii.dataobj, dtype="uint8")
         - np.asanyarray(ntdeobnii.dataobj, dtype="uint8")
     )
-    deobnii.__class__(diff, deobnii.affine, deobnii.header).to_filename("diff.nii.gz")
-    assert np.sqrt((diff[20:-20, 20:-20, 20:-20] ** 2).mean()) < 0.1
+    assert np.sqrt((diff[mask] ** 2).mean()) < 0.1
 
     # Confirm AFNI's rotation of axis is consistent with the one we introduced
     afni_warpdrive_inv = afni._afni_header(
@@ -526,9 +512,34 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y,
     assert np.allclose(afni_warpdrive_inv[:3, :3], R[:3, :3])
 
     # Check nitransforms' estimation of warpdrive with header
-    # Still haven't gotten my head around orientation, this test should not fail
-    # nt_warpdrive_inv = afni._afni_warpdrive(newaff, deobaff, forward=False)
-    # assert not np.allclose(
-    #     afni_warpdrive_inv[:3, :3],
-    #     nt_warpdrive_inv[:3, :3],
-    # )
+    nt_warpdrive_inv = afni._afni_warpdrive(img.affine, forward=False)
+    assert not np.allclose(afni_warpdrive_inv, nt_warpdrive_inv)
+
+
+def _generate_reoriented(path, directions, swapaxes, parameters):
+    img = nb.load(path)
+    shape = np.array(img.shape[:3])
+    hdr = img.header.copy()
+    aff = img.affine.copy()
+    data = np.asanyarray(img.dataobj, dtype="uint8")
+
+    if directions != (1, 1, 1):
+        last_ijk = np.hstack((shape - 1, 1.0))
+        last_xyz = aff @ last_ijk
+        aff = np.diag((*directions, 1)) @ aff
+
+        for ax in range(3):
+            if (directions[ax] == -1):
+                aff[ax, 3] = last_xyz[ax]
+                data = np.flip(data, ax)
+
+    if swapaxes is not None:
+        data = np.swapaxes(data, swapaxes[0], swapaxes[1])
+        aff[list(reversed(swapaxes)), :] = aff[(swapaxes), :]
+
+    R = from_matvec(euler2mat(**parameters), [0.0, 0.0, 0.0])
+
+    newaff = R @ aff
+    hdr.set_qform(newaff, code=1)
+    hdr.set_sform(newaff, code=1)
+    return img.__class__(data, newaff, hdr), R

nitransforms/tests/test_linear.py

@@ -218,7 +218,9 @@ def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orient
     nt_moved_mask.to_filename("nt_resampled_brainmask.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")
+    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)