Skip to content
Snippets Groups Projects
Commit 5fa576c1 authored by Paul McCarthy's avatar Paul McCarthy :mountain_bicyclist:
Browse files

Remove avwutils unit tests - they have been added to the fsl/avwutils project

parent bd7bc0bd
No related branches found
No related tags found
1 merge request!53Remove avwutils unit tests - they have been added to the fsl/avwutils project
Pipeline #19257 skipped
#!/usr/bin/env fslpython
"""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 os
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_data = fixed.get_fdata()
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])))
assert np.all(np.isclose(new_data, data.transpose((0, 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 fixed1')
sprun('fslorient -setqform 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 fixed1')
sprun('imcp fixed1 fixed')
sprun('fslorient -setqformcode 1 fixed')
fixed = nib.load('fixed.nii.gz')
new_data = fixed.get_fdata()
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])))
assert np.all(np.isclose(new_data, data.transpose((0, 2, 1))))
if __name__ == '__main__':
tests = [test_fix_orient,
test_fix_orient_alternate]
result = 0
for test in tests:
try:
os.chdir(sys.argv[1])
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)
#!/usr/bin/env fslpython
import os
import sys
import subprocess as sp
import numpy as np
import nibabel as nib
def create_image(shape, pixdim):
pixdim = list(pixdim)
data = np.random.random(shape).astype(np.float32)
hdr = nib.Nifti1Header()
hdr.set_data_dtype(np.float32)
hdr.set_data_shape(shape)
hdr.set_zooms(pixdim[:len(shape)])
return nib.Nifti1Image(data, np.eye(4), hdr)
def check_image(origimg, changedimg, exppixdim):
gotshape = changedimg.shape
gotpixdim = list(changedimg.header.get_zooms())
gotdata = changedimg.get_fdata()
assert origimg.shape == gotshape
assert gotpixdim == exppixdim
assert np.all(origimg.get_fdata() == gotdata)
def run_tests():
# (image shape, command, exppixdims)
tests = [
((5, 5, 5), '2 2 2', (2, 2, 2)),
((5, 5, 5), '2 2 2 2', (2, 2, 2)),
((5, 5, 5, 5), '2 2 2', (2, 2, 2, 1)),
((5, 5, 5, 5), '2 2 2 2', (2, 2, 2, 2)),
]
for shape, cmd, exppixdim in tests:
imgfile = 'image.nii.gz'
img = create_image(shape, [1] * len(shape))
img.to_filename(imgfile)
try:
sp.run(['fslchpixdim', imgfile] + list(cmd.split()))
check_image(img, nib.load(imgfile), list(exppixdim))
finally:
os.remove(imgfile)
if __name__ == '__main__':
sys.exit(run_tests())
#!/usr/bin/env fslpython
"""Test that fslcomplex produces correct outputs. Specifically tests that the
output files have orientation info that matches that of the input files.
"""
import numpy as np
import nibabel as nib
import subprocess as sp
import os.path as op
import shlex
import sys
import traceback
from fsl.utils.tempdir import tempdir
from fsl.data.image import addExt
def sprun(cmd):
print(f'RUN {cmd}')
sp.run(shlex.split(cmd), check=True)
def imtest(path):
try:
addExt(path)
return True
except Exception:
return False
# radio=[True|False] -> whether or not
# to add a flip on the affine X axis
def create_test_input_data(radio, seed=1):
np.random.seed(seed)
real = np.random.randint(0, 100, (10, 10, 10)).astype(np.float32)
imag = np.random.randint(0, 100, (10, 10, 10)).astype(np.float32)
affine = np.diag([3, 3, 3, 1])
if radio:
affine[0, 0] *= -1
affine[:3, 3] = [37, 20, 30]
real = np.flip(real, 0)
imag = np.flip(imag, 0)
else:
affine[:3, 3] = [10, 20, 30]
comp = real + imag * 1j
abso = np.abs(comp)
phase = np.arctan2(comp.imag, comp.real)
comp4d = np.concatenate((comp.reshape(10, 10, 10, 1),
comp.reshape(10, 10, 10, 1)), 3)
real = nib.Nifti1Image(real, affine)
imag = nib.Nifti1Image(imag, affine)
comp = nib.Nifti1Image(comp, affine)
comp4d = nib.Nifti1Image(comp4d, affine)
abso = nib.Nifti1Image(abso, affine)
phase = nib.Nifti1Image(phase, affine)
real .set_qform(*real .get_sform(coded=True))
imag .set_qform(*comp .get_sform(coded=True))
comp .set_qform(*comp .get_sform(coded=True))
comp4d.set_qform(*comp4d.get_sform(coded=True))
abso .set_qform(*abso .get_sform(coded=True))
phase .set_qform(*phase .get_sform(coded=True))
real .to_filename('real_in.nii.gz')
imag .to_filename('imag_in.nii.gz')
comp .to_filename('complex_in.nii.gz')
comp4d.to_filename('complex4d_in.nii.gz')
abso .to_filename('abs_in.nii.gz')
phase .to_filename('phase_in.nii.gz')
def compare_images(file1, file2):
img1 = nib.load(addExt(file1))
img2 = nib.load(addExt(file2))
data1 = np.asanyarray(img1.dataobj)
data2 = np.asanyarray(img2.dataobj)
assert img1.header['sform_code'] == img2.header['sform_code']
assert img1.header['qform_code'] == img2.header['qform_code']
assert np.all(np.isclose(img1.get_sform(), img2.get_sform()))
assert np.all(np.isclose(img1.get_qform(), img2.get_qform()))
assert np.all(np.isclose(img1.header['pixdim'], img2.header['pixdim']))
assert np.all(np.isclose(data1, data2))
def test_fslcomplex_call(args, infiles, outfiles, radio):
with tempdir():
create_test_input_data(radio)
sprun(f'fslcomplex {args} {infiles} {outfiles}')
for outfile in outfiles.split(' '):
if not imtest(outfile):
continue
prefix = outfile.split('_')[0]
infile = f'{prefix}_in'
print(f'Comparing {outfile} against {infile}')
compare_images(infile, outfile)
if __name__ == '__main__':
tests = [
('-realabs', 'complex_in', 'abs_out'),
('-realphase', 'complex_in', 'phase_out'),
('-realpolar', 'complex_in', 'abs_out phase_out'),
('-realcartesian', 'complex_in', 'real_out imag_out'),
('-complex', 'real_in imag_in', 'complex_out'),
('-complexpolar', 'abs_in phase_in', 'complex_out'),
('-copyonly', 'complex_in', 'complex_out'),
('-complexsplit', 'complex4d_in', 'complex_out 0 0'),
('-complexsplit', 'complex4d_in', 'complex_out 1 1'),
('-complexmerge', 'complex_in complex_in', 'complex4d_out'),
]
result = 0
for radio in [True, False]:
for args, infiles, outfiles in tests:
try:
test_fslcomplex_call(args, infiles, outfiles, radio)
print(f'\nTest fslcomplex {args} {infiles} {outfiles} [radio: {radio}] PASSED')
except Exception as e:
print(f'\nTest fslcomplex {args} {infiles} {outfiles} [radio: {radio}] FAILED')
traceback.print_exc()
result = 1
sys.exit(result)
#!/usr/bin/env fslpython
import os
import os.path as op
import sys
import subprocess as sp
import numpy as np
import nibabel as nib
THISDIR = op.dirname(op.abspath(__file__))
# 2=char, 4=short, 8=int, 16=float, 64=double
DTYPE_MAPPING = {
2 : np.uint8,
4 : np.int16,
8 : np.int32,
16 : np.float32,
64 : np.float64
}
MNI152_2MM_AFFINE = np.array([[-2, 0, 0, 90],
[ 0, 2, 0, -126],
[ 0, 0, 2, -72],
[ 0, 0, 0, 1]])
def validate_result(imgfile, expshape, exppixdim, exporigin, expdtype):
img = nib.load(imgfile)
hdr = img.header
expshape = list(expshape)
exppixdim = list(exppixdim)
exporigin = list(exporigin)
if (len(expshape) == 3) or \
expshape[3] in (0, 1):
expndims = 3
else:
expndims = 4
expaffine = np.diag(exppixdim[:3] + [1])
expaffine[0, 0] *= -1
expaffine[:3, 3] = exporigin
expaffine[0, 3] *= exppixdim[0]
expaffine[1, 3] *= -exppixdim[1]
expaffine[2, 3] *= -exppixdim[2]
assert len(img.shape) == expndims
assert list(img.shape) == expshape[ :expndims]
assert list(hdr.get_zooms()) == exppixdim[:expndims]
assert img.get_data_dtype() == expdtype
assert np.all(np.isclose(img.affine, expaffine))
def create_image(shape, pixdim, origin, dtype):
pixdim = list(pixdim)
affine = np.diag(pixdim[:3] + [1])
affine[:3, 3] = origin
data = np.random.randint(1, 100, shape).astype(dtype)
hdr = nib.Nifti1Header()
hdr.set_data_dtype(dtype)
hdr.set_data_shape(shape)
hdr.set_zooms(pixdim[:len(shape)])
return nib.Nifti1Image(data, affine, hdr)
def test_new_file():
tests = [
# 4th dim of size 0 or 1 should
# result in a 3D image (this
# is coded in validate_result)
' 5 5 5 0 2 2 2 0 0 0 0 2',
' 5 5 5 0 2 2 2 1 0 0 0 2',
' 5 5 5 1 2 2 2 1 0 0 0 2',
' 5 5 5 1 2 2 2 1 0 0 0 4',
' 5 5 5 1 2 2 2 1 0 0 0 8',
' 5 5 5 1 2 2 2 1 0 0 0 16',
' 5 5 5 1 2 2 2 1 1 2 3 64',
' 5 5 5 1 0.5 1.5 1.25 1 0 0 0 2',
' 5 5 5 1 0.5 1.5 1.25 1.5 0 0 0 2',
' 5 5 5 1 0.5 1.5 1.25 0.5 0 0 0 2',
'30 30 30 5 5 10 3 5 10 20 10 2',
]
for test in tests:
args = list(test.split())
imgfile = 'image.nii.gz'
try:
os.remove(imgfile)
except:
pass
print(['fslcreatehd'] + args + [imgfile])
sp.run(['fslcreatehd'] + args + [imgfile])
args = [float(a) for a in args]
expshape = args[:4]
exppixdim = args[4:8]
exporigin = args[8:11]
expdtype = DTYPE_MAPPING[args[11]]
try:
validate_result(
imgfile, expshape, exppixdim, exporigin, expdtype)
finally:
os.remove(imgfile)
def test_existing_file():
# each test is a tuple containing:
# - dtype_of_existing_image
# - shape_of_existing_image
# - exp_shape (set to None if image should not be overwritten)
# - fslcreatehd_args
tests = [
# 4th dim of 0 or 1 should result in a 3D image
(2, (5, 5, 5), None, '5 5 5 0 2 2 2 1 0 0 0 2'),
(2, (5, 5, 5), None, '5 5 5 0 2 2 2 1 1 2 3 2'),
(2, (5, 5, 5), None, '5 5 5 1 2 2 2 2 0 0 0 2'),
# dtype arg should be ignored for existing images
(2, (5, 5, 5), None, '5 5 5 0 2 2 2 1 0 0 0 4'),
(4, (5, 5, 5), None, '5 5 5 0 2 2 2 1 0 0 0 2'),
# data should be overwritten if any
# of the first 3 dims are different,
# or if a 4th dim is specified
(2, (5, 5, 5), (5, 3, 5), '5 3 5 0 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5), (5, 3, 5), '5 3 5 0 2 2 2 2 0 0 0 4'),
(2, (5, 5, 5), (5, 3, 5), '5 3 5 1 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5), (5, 3, 5), '5 3 5 1 2 2 2 2 1 2 3 2'),
(2, (5, 5, 5), (5, 3, 5, 5), '5 3 5 5 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5), (5, 3, 5, 5), '5 3 5 5 2 2 2 2 1 2 3 2'),
# 4D - same rules apply
(2, (5, 5, 5, 5), None, '5 5 5 5 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5, 5), None, '5 5 5 5 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5, 5), None, '5 5 5 5 2 2 2 2 0 0 0 4'),
(4, (5, 5, 5, 5), None, '5 5 5 5 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5, 5), None, '5 5 5 5 2 2 2 2 1 2 3 2'),
# data overwritten if nelements
# change
(2, (5, 5, 5, 5), (5, 5, 5), '5 5 5 0 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5, 5), (5, 5, 5), '5 5 5 1 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5, 5), (5, 3, 5), '5 3 5 0 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5, 5), (5, 3, 5), '5 3 5 1 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5, 5), (5, 3, 5), '5 3 5 1 2 2 2 2 1 2 3 2'),
(2, (5, 5, 5, 5), (5, 5, 5, 2), '5 5 5 2 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5, 5), (5, 3, 5, 5), '5 3 5 5 2 2 2 2 0 0 0 2'),
(2, (5, 5, 5, 5), (5, 3, 5, 5), '5 3 5 5 2 2 2 2 1 2 3 2'),
]
for (dtype, shape, expshape, args) in tests:
imgfile = 'image.nii.gz'
dtype = DTYPE_MAPPING[dtype]
args = list(args.split())
img = create_image(shape, (1, 1, 1, 1), (0, 0, 0), dtype)
img.to_filename(imgfile)
sp.run(['fslcreatehd'] + args + [imgfile])
checkdata = expshape is None
if expshape is None:
expshape = shape
args = [float(a) for a in args]
exppixdim = args[4:8]
exporigin = args[8:11]
try:
validate_result(imgfile, expshape, exppixdim, exporigin, dtype)
if checkdata:
data = np.asanyarray(img.dataobj)
datacopy = np.asanyarray(nib.load(imgfile).dataobj)
assert np.all(data == datacopy)
finally:
os.remove(imgfile)
def test_new_file_xml():
xmlfile = op.join(THISDIR, 'mni2mm.xml')
sp.run(['fslcreatehd', xmlfile, 'mni.nii.gz'])
img = nib.load('mni.nii.gz')
assert img.shape == (91, 109, 91)
assert img.header.get_zooms()[:3] == (2.0, 2.0, 2.0)
assert img.header['intent_code'] == 10
assert img.header['intent_p1'] == 20
assert img.header['intent_p2'] == 30
assert img.header['intent_p3'] == 40
def test_existing_file_xml_same_shape():
xmlfile = op.join(THISDIR, 'mni2mm.xml')
img = create_image((91, 109, 91), (3, 3, 3), (5, 5, 5), np.float32)
img.to_filename('image.nii.gz')
sp.run(['fslcreatehd', xmlfile, 'image.nii.gz'])
result = nib.load('image.nii.gz')
assert result.shape == (91, 109, 91)
assert result.header.get_zooms()[:3] == (2.0, 2.0, 2.0)
# bug in FSL<=6.0.4 caused intent codes to not be set
assert result.header['intent_code'] == 10
assert result.header['intent_p1'] == 20
assert result.header['intent_p2'] == 30
assert result.header['intent_p3'] == 40
assert np.all(result.affine == MNI152_2MM_AFFINE)
# data should be preserved
assert np.all(result.get_fdata() == img.get_fdata())
def test_existing_file_xml_different_shape():
xmlfile = op.join(THISDIR, 'mni2mm.xml')
img = create_image((20, 20, 20), (3, 3, 3), (5, 5, 5), np.float32)
img.to_filename('image.nii.gz')
sp.run(['fslcreatehd', xmlfile, 'image.nii.gz'])
result = nib.load('image.nii.gz')
assert result.shape == (91, 109, 91)
assert result.header.get_zooms()[:3] == (2.0, 2.0, 2.0)
assert result.header['intent_code'] == 10
assert result.header['intent_p1'] == 20
assert result.header['intent_p2'] == 30
assert result.header['intent_p3'] == 40
assert np.all(result.affine == MNI152_2MM_AFFINE)
# data should be cleared
assert np.all(result.get_fdata() == 0)
def main():
outdir = sys.argv[1]
os.chdir(outdir)
test_new_file()
test_existing_file()
test_new_file_xml()
test_existing_file_xml_same_shape()
test_existing_file_xml_different_shape()
if __name__ == '__main__':
sys.exit(main())
<nifti_image
image_offset = '352'
ndim = '3'
nx = '91'
ny = '109'
nz = '91'
nt = '1'
dx = '2'
dy = '2'
dz = '2'
dt = '1'
datatype = '4'
nvox = '902629'
nbyper = '2'
scl_slope = '1'
scl_inter = '0'
intent_code = '10'
intent_p1 = '20'
intent_p2 = '30'
intent_p3 = '40'
intent_name = ''
toffset = '0'
xyz_units = '2'
time_units = '8'
freq_dim = '0'
phase_dim = '0'
slice_dim = '0'
descrip = 'FSL5.0'
aux_file = ''
qform_code = '4'
qfac = '-1'
quatern_b = '0'
quatern_c = '1'
quatern_d = '0'
qoffset_x = '90'
qoffset_y = '-126'
qoffset_z = '-72'
sform_code = '4'
sto_xyz_matrix = '-2 0 0 90 0 2 0 -126 0 0 2 -72 0 0 0 1'
slice_code = '0'
slice_start = '0'
scl_end = '0'
scl_duration = '0'
/>
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment