quasil.in

#!/bin/bash

# QuASIL: QUASAR Bayesian Arterial SpIn Labeling parameter estimation
#
# Michael Chappell, IBME & FMRIB Image Analysis Group
#
# Copyright (c) 2011-2012 University of Oxford
#
# SHCOPYRIGHT

# Make script use local copies of helper scripts/programs in the same
# directory, if present. This allows for multiple versions of the scripts
# to be used, possibly with bundled dependencies
if [ -z "${FSLDEVDIR}" ]; then
    FSLPATH="${FSLDIR}/bin"
else
    FSLPATH="${FSLDEVDIR}/bin:${FSLDIR}/bin"
fi
PATH=`dirname $0`:${FSLPATH}:${PATH}

abspath() {                                               
    # Return an absolute path if the input is relative
    cd "$(dirname "$1")"
    printf "%s/%s\n" "$(pwd)" "$(basename "$1")"
}

Usage() {
    echo "QUASAR Bayesian Inference for Arterial Spin Labelling MRI"
    echo ""
    echo "Usage (optional parameters in {}):"
    echo " -i         : specify data file"
    echo " {-o}       : specify output directory"
    echo " {-m}       : specify brain mask file"
    echo ""
    echo " Extended options:"
    echo " --t1b      : Set the value for T1 of arterial blood {default: 1.6 s}"
    echo " --disp     : include bolus dispersion in the model (gamma kernel)"
    echo " --infertau : estimate bolus duration from data"
    echo " --mfree    : Do model-free (SVD deconvolution) analysis"
    echo " --corrcal  : Correct partial volume effects on the edge of calibration image M0a"
    echo " --alpha    : Inversion efficiency {default: 0.91}"
    echo ""
    echo " Partial volume effects correction options:"
    echo " --pvcorr   : Set partial volume effect correction on. You should provide high resolution partial volume estimates and a structural image."
    #echo " --pvgm     : GM PV Estimates"
    #echo " --pvwm     : WM PV Estimates"
    echo " --fslanat  : Name of the directory containing the output from fsl_anat"
    echo " --t1wm     : T1 for WM {default: 1.1 s}"
    echo ""
    echo " Sequence parameters:"
    echo " --slicedt  : Set the increase in TI with slice {default: 0.035 s}"
    echo " --fa       : Flip angle for LL readout {default: 35 degrees}"
    echo " --lfa      : Lower flip angle for final phase of data {default: 11.7 degrees}"
    echo " --tis      : comma separated list of TI values"
    echo "            {default: 0.04,0.34,0.64,0.94,1.24,1.54,1.84,2.14,2.44,2.74,3.04,3.34,3.64}"
    echo ""
}

Version() {
echo "@GIT_SHA1@ @GIT_DATE@"
exit 0
}

# deal with options
if [ -z $1 ]; then
    Usage
    exit 1
fi

until [ -z $1 ]; do

# look at this option and determine if has an argument specified by an =
option=`echo $1 | sed s/=.*//`
arg="" #specifies if an argument is to be read from next item on command line (=1 is is when = is used)
if [ $option = $1 ]; then
# no argument to this command has been found with it (i.e. after an =)
# if there is an argument it will be the next option
    argument=$2
else
    arg=1
    argument=`echo $1 | sed s/.*=//`
fi
takeargs=0;boolarg="";isbool="";

    case $option in
    -o) outflag=1 outdir=$argument
        shift;;
    -i) inflag=1 infile=$argument #input/data file
        shift;;
    -m) mask=$argument
        shift;;
    --t1b) t1b=$argument
        shift;;
    --t1) t1=$argument
        shift;;
    --t1wm) t1wm=$argument
        shift;;
    --slicedt) slicedt=$argument
        shift;;
    --fa) fa=$argument
        shift;;
    --lfa) lfa=$argument
        shift;;
    --disp) disp=1
        ;;
    --mfree) mfree=1
        ;;
    --edgecorr) isbool=1;
		boolarg=edgecorr;
		;;
    --tis) tis=$argument
        shift;;
    --iform) iform=$argument # to choose the input form of the data
        shift;;
    --tau) tau=$argument
        shift;;
    --corrcal) corrcal=1
        ;;
    --alpha) alpha=$argument
        shift;;
    --infertau) infertau=1
        ;;

    --pvcorr) pvcorr=1
        ;;
    #--pvgm) pvgm_highres=$argument
        #shift;;
    #--pvwm) pvwm_highres=$argument
        #shift;;
    #--s) strim=$argument # structure images for registration and partial volume estimation
        #shift;;
    --fslanat) fslanat=$argument # Directory containing the output from fsl_anat
        shift;;

    --ccmds) calibcmds=$argument
        shift;;
    --debug) debug=1 #debugging option
        ;;
    --version) Version
        ;;
    *)  Usage
        echo "Error! Unrecognised option on command line: $1"
        echo ""
        exit 1;;
    esac

    # sort out a shift required by a command line option that takes arguments
    if [ -z $arg ]; then
	# an argument has been supplied on the command NOT using an =
	if [ $takeargs -eq 1 ]; then
	    shift;
	fi
    fi
    
    if [ ! -z $isbool ]; then
	    # this is an (explicit) boolean setting
	if [ ! -z $arg ]; then
	    # an argument has been supplied on the command using an =
	    # set the variable based on the argument
	    case $argument in
		on) eval $boolarg=1
		    ;;
		off) eval $boolarg=""
		     ;;
		1) eval $boolarg=1
		   ;;
		0) eval $boolarg=""
		   ;;
		*)  Usage
		    echo "Error! Unrecognised setting for boolean option: $1"
		    echo ""
		    exit 1;;
	    esac
	else
	    # no argument has been suppled with this command (NOTE that you cannot supply an arugment to a bool option without an =)
	    # this sets the variable to true
	    eval $boolarg=1;
	fi
    fi

    shift
done

#### --- Procedural ---
asl_file=asl_file
fabber=fabber_asl
asl_mfree=asl_mfree ###~/cproject/asl_mfree/asl_mfree

#### --- Housekeeping ---
# set the output directory here if not specified
if [ -z $outflag ]; then
    echo "Ouput being placed in basil subdirectory of input directory"
    outdir=$indir/quasil;
fi

# Start by looking for the output directory (and create if need be)
count=0
while [ -d $outdir ]; do
    outdir=$outdir"+"
    count=`expr $count + 1`

    if [ $count -gt 20 ]; then
    echo "Error: $outdir too many existing output directories (i.e. shall not add another +)"
    exit
    fi
done
echo "Creating output directory: $outdir"
mkdir $outdir;

# save the starting directory
stdir=`pwd`

# Full path of output directory
outdir=$(abspath $outdir)

# make a temp directory to work in
tmpbase=`tmpnam`
tempdir=${tmpbase}_quasil
mkdir $tempdir

# deal with the TIs
if [ -z $tis ]; then
# default QUASAR list of TIs
tis="0.04,0.34,0.64,0.94,1.24,1.54,1.84,2.14,2.44,2.74,3.04,3.34,3.64"
fi

count=0
tislist=""
thetis=`echo $tis | sed 's:,: :g'`
for ti in $thetis; do
    count=`expr ${count} + 1`
    tislist=`echo $tislist --ti${count}=$ti`
done
# echo "TIs list: $tislist" >> $log
ntis=$count;

if [ -z $iform ]; then
    iform="q"
fi

# parameters
#bolus duration - default 0.64 s
if [ -z $tau ]; then
tau=0.64;
fi

#T1b
if [ -z $t1b ]; then
t1b=1.6;
fi

#T1 - this si the prior value, since T1 will be estimated from the data
if [ -z $t1 ]; then
t1=1.3;
fi

#T1WM
if [ -z $t1wm ]; then
    t1wm=1.1;
fi

# sequence parameters
# slicedt
if [ -z $slicedt ]; then
    slicedt=0.035;
fi
# Flip angles
if [ -z $fa ]; then
    fa=35;
fi
if [ -z $lfa ]; then
    lfa=11.7;
fi
if [ -z $alpha ]; then
    alpha=0.91;
fi

#### --- Pre-processing ---
echo "Pre-processing"
imcp $infile $tempdir/data

if [ ! -z $fslanat ]; then
    cp -R $fslanat $tempdir
fi

cd $tempdir
if [ $iform = "q" ]; then
# input is one big file 13x84 volumes containing raw data (TC pairs) grouped as phases(7) - repeats(6) - tis(13)
# (NB nesting order is from left to right - so that phases are together for one repeat at one TI in this case)
# need to get it into right form for fabber: tis(13) - phases(7), mean over repeats, both subtracted and raw data
    
    # first break out all the TIs
    $asl_file --data=data --ntis=$ntis --ibf=tis --iaf=tc --split=ti
    # now we have 13 tis each with 84 volumes
    
    # Within each TI: separate the phases 
    for ((i=0; i<$ntis; i++)); do
    mkdir ti$i
    tifile=`ls ti$i.nii.gz ti0$i.nii.gz ti00$i.nii.gz 2>/dev/null`
    echo $tifile
    $asl_file --data=$tifile --ntis=7 --ibf=rpt --iaf=tc --split=ti$i/phs
        # NB using asl_file to split the phases (pseudo TIs)
        # leaves TC pairs together
    done

    #now assemble the multiTI files
    phslist=""
    for ((j=0; j<7; j++)); do
    #within each phase
    filelist=""
    for ((i=0; i<$ntis; i++)); do
            #within each TI
        filelist=$filelist" ti$i/phs00$j"
        done
    fslmerge -t aslraw_ph$j $filelist
        #take mean within TI
    $asl_file --data=aslraw_ph$j --ntis=$ntis --ibf=tis --iaf=tc --mean=aslraw_ph$j
    phslist=$phslist" aslraw_ph$j"
     done
    fslmerge -t aslraw $phslist
# data is now in 'f' form
elif [ $iform = "f" ]; then
# data is (already) in 'f' form: one file with 13x7 volumes raw data (TC pairs) grouped as tis(13) - phases(7)
    immv data aslraw
fi

# TC difference
$asl_file --data=aslraw --ntis=$ntis --ibf=tis --iaf=tc --diff --out=asldata
# discard the final (low flip angle) phase from the differenced data
# we do not (currently) use this for the main analysis
nkeep=`expr $ntis \* 6`
fslroi asldata asldata 0 $nkeep
# extract control images
$asl_file --data=aslraw --ntis=$ntis --ibf=tis --iaf=tc --spairs --out=aslraw
immv aslraw_even aslcontrol

if [ -z $mask ]; then
# auto generate mask
    fslmaths aslcontrol -Tmean aslmean
    bet aslmean mask -m
else
    cd "$stdir"
    imcp $mask $tempdir/mask
    cd $tempdir
fi

# copy mask to output for future reference
cd "$stdir" 
imcp $tempdir/mask $outdir/mask
cd $tempdir

#### --- Calibration ---
if [ -z "$calibcmds" ]; then
#voxelwise M0 calibration
    echo "#QUASAR analysis calibration options" > calib_options.txt
    echo "--mask=mask" >> calib_options.txt
    echo "--method=spatialvb --noise=white --param-spatial-priors=MN+" >> calib_options.txt
    echo "--model=satrecov" >> calib_options.txt
    echo "--repeats=1" >> calib_options.txt
    echo "--phases=6" >> calib_options.txt #NB 6 (normal) phases plus 1 LFA phase
    echo $tislist >> calib_options.txt
    echo "--t1=$t1 --FA=$fa --LFA=$lfa" >> calib_options.txt
    echo "--slicedt=$slicedt" >> calib_options.txt
    echo "--link-to-latest" >> calib_options.txt # Here we create a shortcut to the latest results directory
    $fabber --data=aslcontrol --data-order=singlefile --output=calib -@ calib_options.txt

    if [ ! -z $corrcal ]; then
        echo "Correct partial volume effects on the edge of M0 image"
        # First make a copy of the original M0 image
        ${FSLDIR}/bin/imcp calib/mean_M0t calib/mean_M0t_uncorr
        ${FSLDIR}/bin/imcp calib/mean_M0t_uncorr $outdir/M0t_uncorr
        # Use a median filter to correct the artefact
        fslmaths calib/mean_M0t -fmedian calib/mean_M0t_median
        # Erode the edge voxels
        fslmaths calib/mean_M0t_median -ero calib/mean_M0t_ero
        # Extrapolate back the eroded voxels
        $asl_file --data=calib/mean_M0t_ero --ntis=1 --mask=mask --extrapolate --neighbour=5 --out=calib/mean_M0t_corr
        # Rename the corrected M0 image
        ${FSLDIR}/bin/immv calib/mean_M0t_corr calib/mean_M0t
    fi

    # deal with outputs
    ${FSLDIR}/bin/immv calib/mean_T1t calib/T1t
    ${FSLDIR}/bin/immv calib/mean_g calib/g
    ${FSLDIR}/bin/immv calib/mean_M0t calib/M0t

    #fslmaths instruction for calibration (for execution whilst back in starting dir)
    #cinstr=" -div $tempdir/calib/M0t -mul 0.9 " # partition coefficient 0.9
    cinstr=" -div $tempdir/calib/M0t -div 0.9 -div $alpha " # partition coefficient 0.9

    echo $cinstr

    #save calibration results to output directory for reference
    ${FSLDIR}/bin/imcp calib/M0t $outdir/M0t
    ${FSLDIR}/bin/imcp calib/T1t $outdir/T1t

else
    # we have some commands to pass to asl_calib
    cd $stdir #NB run this in the original starting directory 
    asl_calib -c $tempdir/aslcontrol $calibcmds --mode satrecov -o $tempdir/calib --bmask $tempdir/mask --tis $tis --fa $fa --lfa $lfa --nphases

    #fslmaths instruction for calibration (for execution whilst back in starting dir)
    cinstr=" -div `cat $tempdir/calib/M0a.txt` " 

    #return to working in temporary directory
    cd $tempdir

    #save calibration results to output directory for reference
    cp calib/M0a.txt $outdir/M0a.txt
    imcp calib/mean_T1t $outdir/T1t
fi


# Partial Volume Correction options
if [ ! -z $pvcorr ]; then
    # If the output from fsl_anat is provided
    if [ ! -z $fslanat ]; then

    #fslanat=$(abspath $fslanat)
    echo "Copy user specified PVE to temp directory (current directory)"
    echo "Do ASL and Structural image registration to obtain transformation matrix"
    struct=$fslanat/T1
    struct_brain=$fslanat/T1_brain
    # Brain extraction
    # This step may be optional if fsl_anat has done the same job
    # But we perform it here just in case
    ${FSLDIR}/bin/bet $struct $struct_brain
    # Registration
    ${FSLDIR}/bin/asl_reg -i calib/T1t -o reg_dir -s $struct --sbet $struct_brain -c calib/M0t -m mask
    # Transform high resolution PV estimates to low (ASL) resolution
    echo "Downsample to low (ASL) resolution space"
    pvgm_highres=$fslanat/T1_fast_pve_1
    pvwm_highres=$fslanat/T1_fast_pve_2
    ${FSLDIR}/bin/applywarp --ref=calib/T1t --in=$pvgm_highres --out=pvgm_lowres --premat=reg_dir/struct2asl.mat --super --interp=spline --superlevel=4
    ${FSLDIR}/bin/applywarp --ref=calib/T1t --in=$pvwm_highres --out=pvwm_lowres --premat=reg_dir/struct2asl.mat --super --interp=spline --superlevel=4
    pvgm=pvgm_lowres
    pvwm=pvwm_lowres
    
    # Use the estimated T1 image from calibration image to obtain partial volume estimates
    else
        echo "High resolution T1 image not provided. Using estiamted T1 image for partial volume correction..."

        # Upsample the T1 image estimated from calibration and use standard brain as reference
        echo "Upsample the T1 image estimated from calibration and use standard brain as reference"
        standard_brain=$FSLDIR/data/linearMNI/MNI152lin_T1_2mm_brain
        ${FSLDIR}/bin/flirt -in calib/T1t -out T1highres -ref $standard_brain -dof 6 -omat regmat.mat
        #${FSLDIR}/bin/flirt -in calib/T1t -out T1highres -applyisoxfm 1 -ref calib/T1t
        #${FSLDIR}/bin/flirt -in calib/M0t -out PDhighres -applyisoxfm 1 -ref calib/M0t

        # Compute Standard to ASL space transformation matrix
        echo "Compute Standard to ASL space transformation matrix"
        ${FSLDIR}/bin/convert_xfm -omat regmat_inv.mat -inverse regmat.mat

        # Do segmentation
        echo "Tissue segmentation"
        ${FSLDIR}/bin/fast -N -p T1highres

        # Downsample to ASL space
        echo "Downsample to low (ASL) space"
        ${FSLDIR}/bin/applywarp --ref=calib/T1t --in=T1highres_pve_0 --out=pvcsf_lowres --premat=regmat_inv.mat --super --interp=spline --superlevel=4
        ${FSLDIR}/bin/applywarp --ref=calib/T1t --in=T1highres_pve_1 --out=pvwm_lowres --premat=regmat_inv.mat --super --interp=spline --superlevel=4
        ${FSLDIR}/bin/applywarp --ref=calib/T1t --in=T1highres_pve_2 --out=pvgm_lowres --premat=regmat_inv.mat --super --interp=spline --superlevel=4

        #${FSLDIR}/bin/applywarp -i T1highres_pve_2 -r calib/T1t -o pvgm_lowres -s --interp=trilinear
        #${FSLDIR}/bin/applywarp -i T1highres_pve_1 -r calib/T1t -o pvwm_lowres -s --interp=trilinear

        pvgm=pvgm_lowres
        pvwm=pvwm_lowres

    fi

    # Threshold PV estimates - remove voxels below 10%
    ${FSLDIR}/bin/fslmaths $pvgm -thr 0.1 -min 1 $pvgm
    ${FSLDIR}/bin/fslmaths $pvwm -thr 0.1 -min 1 $pvwm
    ${FSLDIR}/bin/fslmaths $pvgm -bin gm_mask
    ${FSLDIR}/bin/fslmaths $pvwm -bin wm_mask

    #save the PVE and masks for future reference
    ${FSLDIR}/bin/imcp $pvgm $outdir/pvgm_lowres
    ${FSLDIR}/bin/imcp $pvwm $outdir/pvwm_lowres
    ${FSLDIR}/bin/imcp gm_mask $outdir/gm_mask
    ${FSLDIR}/bin/imcp wm_mask $outdir/wm_mask

fi



### --- Analysis ---
if [ -z $mfree ]; then
# --- [Model Based] ---
echo "Begin model-based analysis"

    echo "#QUASAR analysis options" > options.txt
    echo "--mask=mask" >> options.txt
    echo "--method=spatialvb" >> options.txt
    echo "--noise=white" >> options.txt
    echo "--model=quasar" >> options.txt
    echo "--inferart" >> options.txt
    echo "--repeats=1" >> options.txt
    echo $tislist >> options.txt
    echo "--t1=$t1 --t1b=$t1b --t1wm=$t1wm --tau=$tau --fa=$fa " >> options.txt
    echo "--slicedt=$slicedt" >> options.txt
    echo "--infert1 ">>options.txt
    echo "--artdir" >> options.txt

# use calibration information within inference
    
    echo "--usecalib ">>options.txt

    # Save model fitting results and residue
    echo "--save-model-fit" >> options.txt
    echo "--save-residuals" >> options.txt

    # Here we create a shortcut to the latest results directory
    echo "--link-to-latest" >> options.txt

    if [ -z $fixt1 ]; then
    t1sp=I
    if [ ! -z $infertau ]; then
        echo "--image-prior6=calib_latest/T1t " >> options.txt
    else
        echo "--image-prior5=calib_latest/T1t " >> options.txt
    fi
    else
    t1sp=N
    fi
    
    if [ -z $pvcorr ]; then
        if [ ! -z $infertau ]; then
           #infer bolus duration
            echo "--infertau --tauboff" >> options.txt #Note we have only a single tau for both arterial and tissue signal (both also share the same dispersion properties)
            echo "--image-prior13=calib_latest/g" >> options.txt
            echo "--param-spatial-priors=MNNAN${t1sp}NNNNNNI" >> options.txt
        else
            # spatial prior list without bolus duration
           echo "--image-prior12=calib_latest/g" >> options.txt
           echo "--param-spatial-priors=MNAN${t1sp}NNNNNNI" >> options.txt
    fi
    else
        # PV correction
        echo "--inferwm " >> options.txt
        echo "--usepve " >> options.txt
        echo "--max-iterations=200" >> options.txt # Maximum 1000 iterations
        if [ ! -z $infertau ]; then
            #infer bolus duration
            echo "--infertau --tauboff" >> options.txt #Note we have only a single tau for both arterial and tissue signal (both also share the same dispersion properties)
            if [ -z $fixt1 ]; then
                echo "--image-prior11=calib_latest/T1t " >> options.txt # WM starts with the same measured T1 value as GM
            fi
            echo "--image-prior12=pvgm_lowres " >> options.txt
            echo "--image-prior13=pvwm_lowres " >> options.txt
            echo "--image-prior19=calib_latest/g" >> options.txt
            echo "--param-spatial-priors=MNNAN${t1sp}NMNN${t1sp}IINNNNNI" >> options.txt
        else
        # spatial prior list without bolus duration
            if [ -z $fixt1 ]; then
                echo "--image-prior9=calib_latest/T1t " >> options.txt # WM starts with the same measured T1 value as GM
            fi
        echo "--image-prior10=pvgm_lowres " >> options.txt
        echo "--image-prior11=pvwm_lowres " >> options.txt
        echo "--image-prior17=calib_latest/g" >> options.txt
        echo "--param-spatial-priors=MNAN${t1sp}NMN${t1sp}IINNNNNI" >> options.txt
        fi
    fi
    

    if [ ! -z $disp ]; then
        # include dispersion in the model (default is gamma)
        if [ ! -z $kern ]; then
            kern="gamma"
        fi
    else
        kern="none"
    fi

    $fabber --data=asldata --data-order=singlefile --disp=$kern --output=full -@ options.txt

   #copy results to output directory
    cd "$stdir"
    if [ -z $pvcorr ]; then
    fslmaths $tempdir/full/mean_ftiss $cinstr -mul 6000 $outdir/perfusion
    else
    fslmaths $tempdir/full/mean_ftiss $cinstr -mul 6000 -mas $tempdir/gm_mask $outdir/perfusion_gm
    fi
    fslmaths $tempdir/full/mean_ftiss $outdir/perfusion_raw
    imcp $tempdir/full/mean_delttiss $outdir/arrival
    fslmaths $tempdir/full/mean_fblood $cinstr $outdir/aCBV
    if [ ! -z $infertau ]; then
    imcp $tempdir/full/mean_tautiss $outdir/bolus_duration
    fi
    if [ ! -z $pvcorr ]; then
    fslmaths $tempdir/full/mean_fwm $cinstr -mul 6000 -mas $tempdir/wm_mask $outdir/perfusion_wm
    fslmaths $tempdir/full/mean_fwm $outdir/perfusion_wm_raw 
    imcp $tempdir/full/mean_deltwm $outdir/arrival_wm
    fi
else
# --- [Model Free] ---
echo "Begin model-free analysis"

# need to separate tissue and arterial signals
# first split up the differenced data into the separate phases (treating as TIs)
$asl_file --data=asldata --ntis=6 --ibf=tis --iaf=diff --split=asldata_ph
fslmaths asldata_ph002 -add asldata_ph005 -mul 0.5 asl_nocrush
fslmaths asldata_ph000 -add asldata_ph001 -add asldata_ph003 -add asldata_ph004 -mul 0.25 asl_tissue
fslmaths asl_nocrush -sub asl_tissue asl_blood

# fit GVF for the AIF
    echo "#QUASAR analysis AIF options" > aifoptions.txt
    echo "--data-order=singlefile" >> aifoptions.txt
    echo "--mask=mask" >> aifoptions.txt
    echo "--method=spatialvb" >> aifoptions.txt
    echo "--noise=white" >> aifoptions.txt
    echo "--model=quasar" >> aifoptions.txt
    echo $tislist >> aifoptions.txt
    echo "--t1=$t1 --t1b=$t1b --tau=$tau --fa=$fa " >> aifoptions.txt
    echo "--slicedt=$slicedt" >> aifoptions.txt
    echo "--repeats=1" >> aifoptions.txt
    echo "--infert1 ">> aifoptions.txt
    echo "--inferart --tissoff" >> aifoptions.txt
    echo "--onephase" >> aifoptions.txt
    echo "--artdir" >> aifoptions.txt
# use calibration information within inference
    echo " --usecalib ">> aifoptions.txt
    echo "--image-prior10=calib/g" >> aifoptions.txt
    echo "--param-spatial-priors=MNNNNNNNNI" >> aifoptions.txt
    # Save model fitting results and residue
    echo "--save-model-fit" >> aifoptions.txt
    echo "--save-residuals" >> aifoptions.txt
    # Here we create a shortcut to the latest results directory
    echo "--link-to-latest" >> aifoptions.txt
    
    $fabber --data=asl_blood --disp=gvf --output=aif -@ aifoptions.txt

    # need aBV image (in absolute units) - to determine what voxels contain viable aif
    fslmaths aif_latest/mean_fblood $cinstr aBV
    
    # need aif shapes (scale aifs by the aBV)
    fslmaths aif_latest/modelfit -div aif_latest/mean_fblood aifs

    #smooth data (a little) before model-free analysis
    fslmaths asl_tissue -s 2.1 asl_tissue

    # do deconvolution

    $asl_mfree --data=asl_tissue --mask=mask --out=modfree --aif=aifs --dt=0.3 --metric=aBV --mthresh=0.012 --tcorrect --t1=1.6 --fa=$fa

    #copy results to output directory
    cd "$stdir"
    fslmaths $tempdir/modfree_magntiude $cinstr -mul 6000 -div $tau $outdir/perfusion 
    # note that in the calibration we have to account for the scaling of the AIF by the bolus duration
    # this is still required (even though we have tau in the model-fitting for the AIF) becuase we normalise the AIF above before deconvolution
    imcp $tempdir/aBV $outdir/aCBV
fi


# clearup
cd "$stdir" # make sure we are back where we started
if [ -z $debug ]; then
    echo "Tidying up"
    rm -r $tempdir
else
mv $tempdir .
fi

echo "QUASIL done"