diff --git a/nitransforms/io/afni.py b/nitransforms/io/afni.py index 71186e05..086309b2 100644 --- a/nitransforms/io/afni.py +++ b/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. diff --git a/nitransforms/tests/test_io.py b/nitransforms/tests/test_io.py index 4379fb31..a2b9eaaf 100644 --- a/nitransforms/tests/test_io.py +++ b/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 diff --git a/nitransforms/tests/test_linear.py b/nitransforms/tests/test_linear.py index fa069eec..4462b212 100644 --- a/nitransforms/tests/test_linear.py +++ b/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) diff --git a/nitransforms/tests/test_nonlinear.py b/nitransforms/tests/test_nonlinear.py index 0662bd06..60ed9e9f 100644 --- a/nitransforms/tests/test_nonlinear.py +++ b/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]