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)