diff --git a/unit_tests/avwutils/fix_orient/feedsRun b/unit_tests/avwutils/fix_orient/feedsRun new file mode 100755 index 0000000000000000000000000000000000000000..a56d814bc001d75c5eabcd87764a4dc3281d8e7e --- /dev/null +++ b/unit_tests/avwutils/fix_orient/feedsRun @@ -0,0 +1,160 @@ +#!/usr/bin/env python +"""Test the standard recipe for fixing a broken orientation, from "The labels +in FSLView are wrong or missing and no conversion to NIfTI can fix this", on +this page: + +https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/Orientation%20Explained + + If this really cannot be fixed at the conversion stage, then it is + possible to fix by the follow (only to be attempted after reading the + following section on background information on NIfTI orientation): + fslorient -deleteorient imagename ; fslswapdim imagename a b c imagename ; + fslorient -setqformcode 1 imagename ; where a b c need to be chosen, + usually by trial and error, so that the image after this is in the same + orientation (as viewed in FSLView) as the MNI152 images. However, unless + you have some way of telling left from right in the image (e.g. by a + vitamin capsule or clearly known structural asymmetry) this method can + easily result in an incorrect left-right labelling and so should not be + attempted. + +The recipe involves: + + 1. resetting the image affines + + 2. transposing the data array until it is in a + RAS orientation + + 3. Enabling the qform affine which, after step + 1, should be a diagonal matrix with pixdims + as scaling parameters. + +Specifically, the recipe is as follows: + + # 1. clear the sform/qform affines + fslorient -deleteorient image + + # 2. rotate the data into RAS orientation (exact + # command depends on original data orientation) + fslswapdim image x z y image + + # 3. reset the qform code. The qform should be + # a diagonal matrix, and the image data array + # should be in RAS orientation, so after this + # step the image should be labelled correctly. + fslorient -setqformcode 1 image + + +This recipe makes some key assumptions: + + - fslorient -deleteorient (step 1) deleting/ + clearing the affines - the qform must be + reset to be a scaling matrix. + + - fslswapdim (step 2) not touching the affines + when the s/qform codes are 0 + +These are required so that, in step 3, the qform is a diagonal matrix with +correct scaling paramaters. + +Some bugs were introduced in FSL 6.0.0 (fixed in FSL 6.0.6) which broke the +above assumptions: + + - fslorient -deleteorient was not clearing the original + affines + - fslswapdim was applying the rotations/flips to both the + data, and to the s/qforms, even if the s/qform codes were + set to 0. +""" + +import numpy as np +import nibabel as nib +import subprocess as sp +import shlex +import sys +import traceback + +def sprun(cmd): + sp.run(shlex.split(cmd), check=True) + + +# Test the original recipe for fixing the +# orientation of an image, as described above +def test_fix_orient(): + + data = np.random.random((10, 10, 10)) + orig_qform = np.array([[2, 0, 0, 10], + [0, 2, 0, -20], + [0, 0, 2, 50], + [0, 0, 0, 1]]) + image = nib.Nifti1Image(data, np.eye(4)) + image.set_sform(np.eye(4), 0) + image.set_qform(orig_qform, 1) + + image.to_filename('orig.nii.gz') + sprun('imcp orig cleared') + sprun('fslorient -deleteorient cleared') + sprun('fslswapdim cleared x z y swapped') + sprun('imcp swapped fixed') + sprun('fslorient -setqformcode 1 fixed') + + fixed = nib.load('fixed.nii.gz') + + new_qform, qcode = fixed.get_qform(coded=True) + _, scode = fixed.get_sform(coded=True) + + # new qform should be a scaling matrix (not + # necessarily equivalent to the original + assert qcode == 1 + assert scode == 0 + assert np.all(np.isclose(new_qform, np.diag([2, 2, 2, 1]))) + + +# An alternate strategy which involves +# explicitly clobbering the qform +def test_fix_orient_alternate(): + # original qform + data = np.random.random((10, 10, 10)) + orig_qform = np.array([[2, 0, 0, 10], + [0, 2, 0, -20], + [0, 0, 2, 50], + [0, 0, 0, 1]]) + image = nib.Nifti1Image(data, np.eye(4)) + image.set_sform(np.eye(4), 0) + image.set_qform(orig_qform, 1) + + image.to_filename('orig.nii.gz') + sprun('imcp orig cleared') + sprun('fslorient -deleteorient cleared') + sprun('fslswapdim cleared x z y swapped') + sprun('imcp swapped fixed') + sprun('fslorient -setqform 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 fixed') + sprun('fslorient -setqformcode 1 fixed') + + fixed = nib.load('fixed.nii.gz') + + new_qform, qcode = fixed.get_qform(coded=True) + _, scode = fixed.get_sform(coded=True) + + # new qform should be a scaling matrix (not + # necessarily equivalent to the original + assert qcode == 1 + assert scode == 0 + assert np.all(np.isclose(new_qform, np.diag([2, 2, 2, 1]))) + + + +if __name__ == '__main__': + tests = [test_fix_orient, + test_fix_orient_alternate] + + result = 0 + for test in tests: + try: + test() + print(f'\nTest {test.__name__} PASSED') + except Exception as e: + print(f'\nTest {test.__name__} FAILED') + traceback.print_exc() + result = 1 + + sys.exit(result)