From 08e56c4632edf175221ec9c263dc1b77efb99e92 Mon Sep 17 00:00:00 2001
From: Taylor Hanayik <hanayik@gmail.com>
Date: Fri, 16 Jun 2023 14:08:31 +0100
Subject: [PATCH] add previous avwutils tests and new fslstats test
---
.gitignore | 3 +-
tests/avscale/feedsRun | 116 +++++++++++++++
tests/fix_orient/feedsRun | 166 ++++++++++++++++++++++
tests/fslchpixdim/feedsRun | 58 ++++++++
tests/fslcomplex/feedsRun | 132 +++++++++++++++++
tests/fslcreatehd/feedsRun | 265 +++++++++++++++++++++++++++++++++++
tests/fslcreatehd/mni2mm.xml | 44 ++++++
tests/fslstats/feedsRun | 180 ++++++++++++++++++++++++
8 files changed, 963 insertions(+), 1 deletion(-)
create mode 100755 tests/avscale/feedsRun
create mode 100755 tests/fix_orient/feedsRun
create mode 100755 tests/fslchpixdim/feedsRun
create mode 100755 tests/fslcomplex/feedsRun
create mode 100755 tests/fslcreatehd/feedsRun
create mode 100644 tests/fslcreatehd/mni2mm.xml
create mode 100755 tests/fslstats/feedsRun
diff --git a/.gitignore b/.gitignore
index 600d2d3..49b91c2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,2 @@
-.vscode
\ No newline at end of file
+.vscode
+pyfeeds-out
\ No newline at end of file
diff --git a/tests/avscale/feedsRun b/tests/avscale/feedsRun
new file mode 100755
index 0000000..32062a0
--- /dev/null
+++ b/tests/avscale/feedsRun
@@ -0,0 +1,116 @@
+#!/usr/bin/env fslpython
+
+# test avscale with different input affines
+
+
+import sys
+import os
+
+import fsl.utils.run as run
+
+import numpy as np
+
+
+PI = np.pi
+PION2 = np.pi / 2
+
+
+# Each test has the form (input-affine, expected-output)
+# where expected-output comprises:
+# - scales
+# - translations
+# - rotations
+# - skews
+# - L/R orientation (preserved or swapped)
+tests = [
+
+ ([[1, 0, 0, 0],
+ [0, 1, 0, 0],
+ [0, 0, 1, 0],
+ [0, 0, 0, 1]], [(1, 1, 1), (0, 0, 0), (0, 0, 0), (0, 0, 0), 'preserved']),
+ ([[2, 0, 0, 0],
+ [0, 2, 0, 0],
+ [0, 0, 2, 0],
+ [0, 0, 0, 1]], [(2, 2, 2), (0, 0, 0), (0, 0, 0), (0, 0, 0), 'preserved']),
+ ([[-2, 0, 0, 0],
+ [ 0, 2, 0, 0],
+ [ 0, 0, 2, 0],
+ [ 0, 0, 0, 1]], [(2, 2, 2), (0, 0, 0), (0, 0, PI), (0, 0, 0), 'swapped']),
+ ([[0, 0, 1, 0],
+ [0, 1, 0, 0],
+ [1, 0, 0, 0],
+ [0, 0, 0, 1]], [(1, 1, 1), (0, 0, 0), (0, -PION2, 0), (0, 0, 0), 'swapped']),
+ ([[ 0, 0, -1, 0],
+ [ 0, -1, 0, 0],
+ [-1, 0, 0, 0],
+ [ 0, 0, 0, 1]], [(1, 1, 1), (0, 0, 0), (-PI, PION2, 0), (0, 0, 0), 'preserved']),
+ ([[0, 1, 0, 0],
+ [1, 0, 0, 0],
+ [0, 0, 1, 0],
+ [0, 0, 0, 1]], [(1, 1, 1), (0, 0, 0), (0, 0, PION2), (0, 0, 0), 'swapped']),
+ ([[1, 0, 0, 0],
+ [0, 0, 1, 0],
+ [0, 1, 0, 0],
+ [0, 0, 0, 1]], [(1, 1, 1), (0, 0, 0), (PION2, 0, 0), (0, 0, 0), 'swapped']),
+ ([[-2, 0, 0, 90],
+ [ 0, 2, 0, -126],
+ [ 0, 0, 2, -72],
+ [ 0, 0, 0, 1]], [(2, 2, 2), (90, -126, -72), (0, 0, PI), (0, 0, 0), 'swapped']),
+
+]
+
+def read_avscale_output(output):
+ lines = output.split('\n')
+ # scales
+ # translations
+ # rotations
+ # skews
+ # l/r orientation (preserved or swapped)
+ scales = [l for l in lines if l.startswith('Scales ')] [0]
+ translations = [l for l in lines if l.startswith('Translations ')] [0]
+ rotations = [l for l in lines if l.startswith('Rotation Angles')][0]
+ skews = [l for l in lines if l.startswith('Skews ')] [0]
+ orient = [l for l in lines if l.startswith('Left-Right ')] [0]
+
+ scales = [float(s) for s in scales .split()[-3:]]
+ translations = [float(s) for s in translations.split()[-3:]]
+ rotations = [float(s) for s in rotations .split()[-3:]]
+ skews = [float(s) for s in skews .split()[-3:]]
+ orient = orient.split()[-1]
+
+ return {
+ 'scales' : np.array(scales),
+ 'translations' : np.array(translations),
+ 'rotations' : np.array(rotations),
+ 'skews' : np.array(skews),
+ 'orient' : orient.lower(),
+ }
+
+
+def test_avscale():
+
+ for affine, expected in tests:
+ affine = np.array(affine)
+ np.savetxt('affine.mat', affine, '%0.6f')
+ print(affine)
+ output = run.runfsl('avscale --allparams affine.mat')
+ output = read_avscale_output(output)
+
+ scales, translations, rotations, skews, orient = expected
+
+ try:
+ assert np.all(np.isclose(scales, output['scales']))
+ assert np.all(np.isclose(translations, output['translations']))
+ assert np.all(np.isclose(rotations, output['rotations']))
+ assert np.all(np.isclose(skews, output['skews']))
+ assert orient == output['orient']
+ except AssertionError:
+ print(affine)
+ print(output)
+ raise
+
+
+if __name__ == '__main__':
+ if len(sys.argv) > 1:
+ os.chdir(sys.argv[1])
+ test_avscale()
diff --git a/tests/fix_orient/feedsRun b/tests/fix_orient/feedsRun
new file mode 100755
index 0000000..0e0e575
--- /dev/null
+++ b/tests/fix_orient/feedsRun
@@ -0,0 +1,166 @@
+#!/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)
diff --git a/tests/fslchpixdim/feedsRun b/tests/fslchpixdim/feedsRun
new file mode 100755
index 0000000..56bcee4
--- /dev/null
+++ b/tests/fslchpixdim/feedsRun
@@ -0,0 +1,58 @@
+#!/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())
diff --git a/tests/fslcomplex/feedsRun b/tests/fslcomplex/feedsRun
new file mode 100755
index 0000000..a84cc1d
--- /dev/null
+++ b/tests/fslcomplex/feedsRun
@@ -0,0 +1,132 @@
+#!/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)
diff --git a/tests/fslcreatehd/feedsRun b/tests/fslcreatehd/feedsRun
new file mode 100755
index 0000000..339553d
--- /dev/null
+++ b/tests/fslcreatehd/feedsRun
@@ -0,0 +1,265 @@
+#!/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())
diff --git a/tests/fslcreatehd/mni2mm.xml b/tests/fslcreatehd/mni2mm.xml
new file mode 100644
index 0000000..06b8043
--- /dev/null
+++ b/tests/fslcreatehd/mni2mm.xml
@@ -0,0 +1,44 @@
+<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'
+/>
diff --git a/tests/fslstats/feedsRun b/tests/fslstats/feedsRun
new file mode 100755
index 0000000..c6d1305
--- /dev/null
+++ b/tests/fslstats/feedsRun
@@ -0,0 +1,180 @@
+#!/usr/bin/env fslpython
+
+# test fslstats with a bunch of different options
+#
+# Usage: fslstats [preoptions] <input> [options]
+#
+# preoption -t will give a separate output line for each 3D volume of a 4D timeseries
+# preoption -K < indexMask > will generate seperate n submasks from indexMask, for indexvalues 1..n where n is the maximum index value in indexMask, and generate statistics for each submask
+# Note - options are applied in order, e.g. -M -l 10 -M will report the non-zero mean, apply a threshold and then report the new nonzero mean
+
+# -l <lthresh> : set lower threshold
+# -u <uthresh> : set upper threshold
+# -r : output <robust min intensity> <robust max intensity>
+# -R : output <min intensity> <max intensity>
+# -e : output mean entropy ; mean(-i*ln(i))
+# -E : output mean entropy (of nonzero voxels)
+# -v : output <voxels> <volume>
+# -V : output <voxels> <volume> (for nonzero voxels)
+# -m : output mean
+# -M : output mean (for nonzero voxels)
+# -s : output standard deviation
+# -S : output standard deviation (for nonzero voxels)
+# -w : output smallest ROI <xmin> <xsize> <ymin> <ysize> <zmin> <zsize> <tmin> <tsize> containing nonzero voxels
+# -x : output co-ordinates of maximum voxel
+# -X : output co-ordinates of minimum voxel
+# -c : output centre-of-gravity (cog) in mm coordinates
+# -C : output centre-of-gravity (cog) in voxel coordinates
+# -p <n> : output nth percentile (n between 0 and 100)
+# -P <n> : output nth percentile (for nonzero voxels)
+# -a : use absolute values of all image intensities
+# -n : treat NaN or Inf as zero for subsequent stats
+# -k <mask> : use the specified image (filename) for masking - overrides lower and upper thresholds
+# -d <image> : take the difference between the base image and the image specified here
+# -h <nbins> : output a histogram (for the thresholded/masked voxels only) with nbins
+# -H <nbins> <min> <max> : output a histogram (for the thresholded/masked voxels only) with nbins and histogram limits of min and max
+
+# Note - thresholds are not inclusive ie lthresh<allowed<uthresh
+
+import sys
+import os
+import fsl.utils.run as run
+import numpy as np
+import nibabel as nib
+
+# set the random seed so that the tests are reproducible
+np.random.seed(0)
+
+options = [
+ {"option": "-r", "expected": "0.018533 0.978824"},
+ {"option": "-R", "expected": "0.000546 0.999809"},
+ {"option": "-e", "expected": "0.906103"},
+ {"option": "-E", "expected": "0.906103"},
+ {"option": "-v", "expected": "1000 1000.000000"},
+ {"option": "-V", "expected": "1000 1000.000000"},
+ {"option": "-m", "expected": "0.495922"},
+ {"option": "-M", "expected": "0.495922"},
+ {"option": "-s", "expected": "0.290744"},
+ {"option": "-S", "expected": "0.290744"},
+ {"option": "-w", "expected": "0 10 0 10 0 10 0 1"},
+ {"option": "-x", "expected": "9 7 4"},
+ {"option": "-X", "expected": "8 2 1"},
+ {"option": "-c", "expected": "4.498706 4.545873 4.613188"},
+ {"option": "-C", "expected": "4.498706 4.545873 4.613188"},
+ {"option": "-p 50", "expected": "0.482584"},
+ {"option": "-P 50", "expected": "0.482584"},
+]
+
+# create a list of tests and expected results
+tests = []
+for option in options:
+ tests.append(
+ {
+ "preoptions": "",
+ "mask": "",
+ "options": option["option"],
+ "expected": option["expected"],
+ }
+ )
+
+tests_with_preoptions = [
+ {
+ "preoptions": "-K",
+ "mask": "mask.nii.gz",
+ "options": "-m",
+ "expected": "0.948930 \n0.947671 \n1.003258 \n1.010696",
+ },
+ {
+ "preoptions": "-t -K",
+ "mask": "mask.nii.gz",
+ "options": "-m",
+ "expected": "0.459736 \n0.476035 \n0.504080 \n0.549485 \n0.489194 \n0.471636 \n0.499178 \n0.461211",
+ },
+ {
+ "preoptions": "-t",
+ "mask": "",
+ "options": "-m",
+ "expected": "0.496675 \n0.487950",
+ },
+]
+
+# taken from fslchpixdim test
+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 create_mask(input_img, n_rois=4):
+ '''
+ Create a mask image from the input image.
+ The mask image will have n_rois different values with random sizes randomly placed in the image.
+ '''
+ data = input_img.get_fdata()
+ mask = np.zeros_like(data)
+ for i in range(n_rois + 1):
+ roi_size = np.random.randint(1, data.size)
+ roi_value = i
+ roi = np.random.choice(data.size, roi_size, replace=False)
+ mask.flat[roi] = roi_value
+ return nib.Nifti1Image(mask, input_img.affine, input_img.header)
+
+
+def test_fslstats():
+ imgfile = 'test.nii.gz'
+ shape = (10,10,10)
+ img = create_image(shape, [1] * len(shape))
+ img.to_filename(imgfile)
+ mask = create_mask(img)
+ mask.to_filename('mask.nii.gz')
+
+ # run the tests witoout preoptions
+ for test in tests:
+ cmd = f"fslstats {test['preoptions']} {test['mask']} {imgfile} {test['options']}"
+ # remove any double spaces from empty test options
+ cmd = " ".join(cmd.split())
+ print(cmd)
+ output = run.runfsl(cmd)
+ # trim any trailing whitespace from the output (fslstats adds a newline at the end)
+ output = output.rstrip()
+ try:
+ assert output == test["expected"]
+ except AssertionError:
+ print(cmd)
+ print(output)
+ print(test["expected"])
+ raise
+
+ # run the tests with preoptions
+ imgfile = 'test_4d.nii.gz'
+ shape = (10,10,10, 2)
+ img = create_image(shape, [1] * len(shape))
+ img.to_filename(imgfile)
+ for test in tests_with_preoptions:
+ cmd = f"fslstats {test['preoptions']} {test['mask']} {imgfile} {test['options']}"
+ # remove any double spaces from empty test options
+ cmd = " ".join(cmd.split())
+ print(cmd, flush=True)
+ output = run.runfsl(cmd)
+ # trim any trailing whitespace from the output (fslstats adds a newline at the end)
+ output = output.rstrip()
+ try:
+ assert output == test["expected"]
+ except AssertionError:
+ print('Failed test')
+ print(cmd)
+ print('Output')
+ print(output)
+ print('Expected')
+ print(test["expected"])
+ raise
+
+
+if __name__ == "__main__":
+ if len(sys.argv) > 1:
+ os.chdir(sys.argv[1])
+ test_fslstats()
--
GitLab