add previous avwutils tests and new fslstats test

.gitignore

-.vscode
\ No newline at end of file
+.vscode
+pyfeeds-out
\ No newline at end of file

tests/avscale/feedsRun

+#!/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()

tests/fix_orient/feedsRun

+#!/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)

tests/fslchpixdim/feedsRun

+#!/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())

tests/fslcomplex/feedsRun

+#!/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)

tests/fslcreatehd/feedsRun

+#!/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())

tests/fslcreatehd/mni2mm.xml

+<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'
+/>

tests/fslstats/feedsRun

+#!/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()