Skip to content
Snippets Groups Projects
feedsRun 6.61 KiB
#!/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()