#!/bin/bash

# OXFORD_ASL: Converts ASL images in to perfusion maps
#
# Michael Chappell, FMRIB Image Analysis Group & IBME
#
# Copyright (c) 2008-2016 University of Oxford
#
#   Part of FSL - FMRIB's Software Library
#   http://www.fmrib.ox.ac.uk/fsl
#   fsl@fmrib.ox.ac.uk
#
#   Developed at FMRIB (Oxford Centre for Functional Magnetic Resonance
#   Imaging of the Brain), Department of Clinical Neurology, Oxford
#   University, Oxford, UK
#
#
#   LICENCE
#
#   FMRIB Software Library, Release 6.0 (c) 2018, The University of
#   Oxford (the "Software")
#
#   The Software remains the property of the Oxford University Innovation
#   ("the University").
#
#   The Software is distributed "AS IS" under this Licence solely for
#   non-commercial use in the hope that it will be useful, but in order
#   that the University as a charitable foundation protects its assets for
#   the benefit of its educational and research purposes, the University
#   makes clear that no condition is made or to be implied, nor is any
#   warranty given or to be implied, as to the accuracy of the Software,
#   or that it will be suitable for any particular purpose or for use
#   under any specific conditions. Furthermore, the University disclaims
#   all responsibility for the use which is made of the Software. It
#   further disclaims any liability for the outcomes arising from using
#   the Software.
#
#   The Licensee agrees to indemnify the University and hold the
#   University harmless from and against any and all claims, damages and
#   liabilities asserted by third parties (including claims for
#   negligence) which arise directly or indirectly from the use of the
#   Software or the sale of any products based on the Software.
#
#   No part of the Software may be reproduced, modified, transmitted or
#   transferred in any form or by any means, electronic or mechanical,
#   without the express permission of the University. The permission of
#   the University is not required if the said reproduction, modification,
#   transmission or transference is done without financial return, the
#   conditions of this Licence are imposed upon the receiver of the
#   product, and all original and amended source code is included in any
#   transmitted product. You may be held legally responsible for any
#   copyright infringement that is caused or encouraged by your failure to
#   abide by these terms and conditions.
#
#   You are not permitted under this Licence to use this Software
#   commercially. Use for which any financial return is received shall be
#   defined as commercial use, and includes (1) integration of all or part
#   of the source code or the Software into a product for sale or license
#   by or on behalf of Licensee to third parties or (2) use of the
#   Software or any derivative of it for research with the final aim of
#   developing software products for sale or license to a third party or
#   (3) use of the Software or any derivative of it for research with the
#   final aim of developing non-software products for sale or license to a
#   third party, or (4) use of the Software to provide any service to an
#   external organisation for which payment is received. If you are
#   interested in using the Software commercially, please contact Oxford
#   University Innovation ("OUI"), the technology transfer company of the
#   University, to negotiate a licence. Contact details are:
#   fsl@innovation.ox.ac.uk quoting Reference Project 9564, FSL.
export LC_ALL=C

# 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}

# Standard options for applying transformations
APPLYXFM_OPTS="-interp trilinear -paddingsize 1"
APPLYWARP_OPTS="--interp=trilinear --paddingsize=1 --super --superlevel=a"

# Save args for reference
ARGS="$@"

Usage() {
    echo "OXFORD_ASL"
    echo "Calculate perfusion maps from ASL data"
    echo ""
    echo "Usage: oxford_asl -i <asl_data> -o <output_dir_name> [options]"
    echo ""
    echo ""
    echo " Main options:"
    echo " -m         : mask (in native space of ASL data) - {default: automatically generated}"
    echo " --spatial  : Use adaptive spatial smoothing on perfusion {default: on}"
    echo " --wp       : analysis that conforms to the 'white paper' (Alsop et al. 2014)"
    echo " --mc       : apply motion correction using mcflirt"
    echo " --more     : see extended options and usage notes"
    echo ""
    echo " Acquisition/Data specific"
    echo " --iaf       : input ASl format: diff,tc,ct {default: diff}"
    echo " --ibf       : input block format (for multi-TI): rpt,tis {default: rpt}"
    echo " --plds      : comma separated list of post-labelling delay times in second, e.g. --plds 0.2,0.4,0.6"
    echo " --tiimg     : 4D image containing voxelwise TI values"
    echo " --casl      : ASL acquisition is  pseudo cASL (pcASL) rather than pASL"
    echo " --artsupp   : Arterial suppression (vascular crushing) was used"
    echo " --bolus     : Bolus duration - {default: 1}"
    echo " --bat       : Bolus arrival time - {default: 0.7 (pASL); 1.3 (cASL)}"
    echo " --t1        : Tissue T1 value - {default: 1.3}"
    echo " --t1b       : Blood T1 value - {default: 1.65}"
    echo " --slicedt   : Timing difference between slices - {default: 0}"
    echo " --sliceband : Number of slices per band in a multi-band setup - {default: n/a}."
    echo " --fa        : Flip angle for Look-Locker readout correction."
    echo ""
    echo "Structural image (optional) (see also Registration)"
    echo " --fslanat=<dir> : An fsl_anat directory from structural image"
    echo " Alternatively:"
    echo " -s              : structural image (whole head)"
    echo " --sbrain        : structural image (already BETed)"
    echo " {--fastsrc=<image_stub>} : images from a FAST segmenation - if not set FAST will be run on structural"
    echo " Options in all cases:"
    echo " --senscorr      : use bias field (from segmentation) for sensitivity correction"
    echo ""
    echo " Calibration"
    echo " --M0        : (single) precomputed M0 value (e.g. from having run a separate calibration)"
    echo " --alpha     : Inversion efficiency - {default: 0.98 (pASL); 0.85 (cASL)}"
    echo "  To supply calibration data:"
    echo "  -c          : M0 calibration image (proton density or mean control image)"
    echo "  --tr        : TR of calibration data - {default: 3.2 s}"
    echo "  --cmethod   : Calibration method: "
    echo "                single - default if structural image is supplied"
    echo "                 M0 value will be calculated within automatically created CSF mask"
    echo "                voxel  - default if no structral image is supplied"
    echo "                 voxelwise M0 values derrived from calibration image"
    echo "                both   - Do single and voxelwise calibaration methods, outputing results"
    echo "                 in separate subdirectories of output directories"
    echo ""
}

Usage_extended() {
    echo "Extended options (all optional):"
    echo " Analysis"
    echo " --artoff     : Do not infer arterial signal - (same as --artsupp)"
    echo " --artonly    : Only infer the arterial signal"
    echo " --fixbolus   : Bolus duration is fixed, e.g. by QUIPSSII or CASL (otheriwse it will be estimated)"
    echo "                 {default: on (cASL); off (pASL)"
    echo " --fixbat     : Fix the bolus arrival time (to value specified by --bat)"
    echo " --batim <img>: Voxelwise BAT (ATT) in seconds estimates in native (ASL) space"
    echo " --infert1    : Incorporate uncertainty in T1 values in analysis"
    echo " --t1im <img> : Voxelwise T1 estimates in native (ASL) space"
    echo " --fulldata   : Never average multiple measurements at each TI"
    echo " --noiseprior : Use an informative prior for the noise estimation"
    echo "   --noisesd  : Set a custom noise std. dev. for the nosie prior"
    echo "   --noisesd  : Set a custom noise std. dev. for the nosie prior"
    echo " --qc-output  : Include extended quality control output"

    echo " Registration (requires structural image)"
    echo " --asl2struc     : transformation matrix from asl data to structural image"
    echo "                   (skips registration)"
    echo " --regfrom       : image to use a basis for inital registration (already BETed)"
    echo "                   (must be same resolution as ASL data)"
    echo " --regfrom-method: If --regfrom is not specified, how to choose the regfrom image."
    echo "                   Choices are 'pwi' (mean differenced ASL data), 'calib' for "
    echo "                   calibration image, 'mean' for mean ASL data. If no regfrom image"
    echo "                   or method specified, the calibration image is used for subtracted"
    echo "                   ASL data, and the mean ASL image for label/control pairs. The "
    echo "                   image is brain extracted before use in registration."
    echo " --reg-init-bbr  : Use BBR in the initial registration. This makes sense mostly when"
    echo "                   using --regfrom-method=pwi"
    echo " --regfrom-bet-thresh: Fractional threshold to use when brain-extracting images for registration (default 0.2)"
#    echo " -t          : structural to standard space transformation matrix"
#    echo "               (requires structural image)"
#    echo " --structout : (also) save the maps in structural space (for use with -t)"
#    echo " -S          : standard brain image - {default: MNI152_T1_2mm}"
    echo " -r          : low resolution structural image (already BETed)"

    echo " Calibration"
    echo "  --cgain     : Relative gain between calibration and ASL image - {default: 1}"
    echo "   extended options for voxel"
    echo "    --t1t      : Assumed T1 of tissue (only used if TR<5 s) {default: 1.3 s}"
    echo "   extended options for single"
    echo "    --tissref  : Specify tissue reference type: csf, wm, gm or none {default: csf}"
    echo "    --csf      : Manually specify csf mask for calibration"
    echo "    --t1csf    : T1 of CSF (s) - {for default see asl_calib}"
    echo "    --te       : Echo time for the readout (to correct for T2/T2* effect in calibration)"
    echo "    --t2star   : Correct for T2* rather than T2"
    echo "    --t2csf    : T2 of CSF (ms) - {for default see asl_calib}"
    echo "    --t2bl     : T2 of blood (ms) - {for default see asl_calib}"
    echo "    --cref     : Supply a reference image for sensitivity correction"

    echo " Distortion correction using fieldmap (see epi_reg):"
    echo "  requires structural image to be provided"
    echo " --fmap=<image>         : fieldmap image (in rad/s)"
    echo " --fmapmag=<image>      : fieldmap magnitude image - wholehead extracted"
    echo " --fmapmagbrain=<image> : fieldmap magnitude image - brain extracted"
    echo " --echospacing=<val>    : Effective EPI echo spacing (sometimes called dwell time) - in seconds"
    echo " --pedir=<dir>          : Phase encoding direction, dir = x/y/z/-x/-y/-z. This is relative to the axes"
    echo "                          in the image voxel array which usually (but not always) corresponds to RAS."
    echo "                          If in doubt, compare output with and without distortion correction."
    echo " {--nofmapreg}          : do not perform registration of fmap to T1 (use if fmap already in T1-space) "
    echo " Distortion correction using phase-encode-reversed calibration image (TOPUP):"
    echo " --cblip=<image>        : phase-encode-reversed (blipped) calibration image"
    echo " --echospacing=<val>    : Effective EPI echo spacing (sometimes called dwell time) - in seconds"
    echo " --pedir=<dir>          : phase encoding direction, dir = x/y/z/-x/-y/-z"
    echo " Gradient distortion correction:"
    echo " --gdcwarp=<image>      : Warp image for gradient distortion correction - will be combined with"
    echo "                          fieldmap or TOPUP distortion correction"
    echo ""
 
    echo "Partial volume correction"
    echo " --pvcorr    : Do partial volume correction"
    echo "  PV estimates will be taken from:"
    echo "   fsl_anat dir (--fslanat), if supplied"
    echo "   exising fast segmentation (--fastsrc), if supplied"
    echo "   FAST segmenation of structural (if using -s and --sbet)"
    echo "   User supplied PV estimates (--pvgm, --pvwm)"   
    echo "   --pvgm    : Partial volume estimates for GM"
    echo "   --pvwm    : Partial volume estimates for WM"

    echo " Epochs "
    echo " --elen      : Length of each epoch in TIs"
    echo " --eol       : Overlap of each epoch in TIs- {deafult: 0}"
    echo ""
    echo " Miscellaneous "
    echo " --model-options=<file>  : File containing additional model options to be passed to BASIL/Fabber"
    echo "                           This is an advanced setting for model-specific options ONLY"
    echo "                           Do NOT include any standard oxford_asl options here"
    echo " --region-analysis       : Run region analysis using ROIs (by default from Harvard-Oxford standard atlas)"
    echo "                           This requires either an FSL_ANAT directory input or a structural->standard space transformation"
    echo " --region-analysis-save-rois : Save ROIs used for region analysis"
    echo " --region-analysis-atlas : Label image containing ROIs to use for region analysis"
    echo " --region-analysis-atlas-labels : Text file containing names of regions in the label image"
    echo " --region-analysis-psf   : Text file containing Z-blurring point spread function for blurred region analysis"
    echo " --gm-thresh=<val>       : Probability threshold for a voxel to be considered 'pure' GM (default 0.8)"
    echo " --wm-thresh=<val>       : Probability threshold for a voxel to be considered 'pure' WM (default 0.9)"
    echo ""
    echo " Notes:"
    echo " Input data - Any combination of label (tag) control images is accepted, see asl_file for more detailed"
    echo "              usage of the input options"
    echo " Analysis  -  If only one TI is specified then the following options are set:"
    echo "              --fixbat --fixbolus --artoff --fulldata"
    echo "              White Paper mode - to conform to the 'white paper' the follow options are also set:"
    echo "              --t1=1.65 --exch=simple --cmethod=voxel"
    echo " Registration - Results are saved in native space by default in case you need to revisit this."
    echo "                By default the perfusion image is used with BBR for the final registration"
    echo "                An inital registration is performed using mean of the ASL data"
    echo "                An alternative base image for (inital) registration can be supplied with --regfrom"
    echo "                Final (BBR) registration can be turned off using --finalonly"
    echo "                For custom registration: use asl_reg with native space results."
    echo " Calibration - Performed using asl_calib using CSF as a reference ('longtr' mode) where structural image is supplied"
    echo "               Voxelwise calibration is performed in 'white paper' mode or in the absence of structural image"
    echo "               For custom calibration: do not set M0 image and then use asl_calib separately"
    echo " Masking - for processing purposes a brain mask is applied to the data, this will be"
    echo "            derrived from (in order of preference):"
    echo "            > Brain extracted structural and (inital) registration"
    echo "            > mean ASL data"
    echo "            > 'regfrom' image if supplied"

}

Version() {
echo "v4.0.29-5-g09a5cba5-dirty Thu Jul 13 14:20:21 2023"
exit 0
}

Output() {
# Function that deals with creating the right output files for a given parameter
    param=$1 #parameter as know by basil (and related scripts)
    shift
    parname=$1 #parameter name as called in the output directory
    shift
    subdir_in=$1 # a subdirectory name in which to find and place results
    shift
    subdir_out=$1 # Subdirectory to place results

Log "Output: $param $parname $subdir_in $subdir_out"
if [ ! -z $nativeout ]; then
    if [ ! -d  $outdir/native_space/$subdir_out ]; then mkdir -p $outdir/native_space/$subdir_out; fi
    imcp $tempdir/$subdir_in/$param $outdir/native_space/$subdir_out/$parname
fi
if [ ! -z $structout ]; then
    if [ ! -d  $outdir/struct_space/$subdir_out ]; then mkdir -p $outdir/struct_space/$subdir_out; fi
    flirt -in $tempdir/$subdir_in/$param -applyxfm -init $tempdir/asl2struct.mat -ref $tempdir/struc -out $outdir/struct_space/$subdir_out/$parname $APPLYXFM_OPTS
fi
if [ ! -z $stdout ]; then
    if [ ! -d  $outdir/std_space/$subdir_out ]; then mkdir -p $outdir/std_space/$subdir_out; fi
    if [ ! -z $warp ]; then
        applywarp --in=$tempdir/$subdir_in/$param --out=$outdir/std_space/$subdir_out/$parname --ref=$std --premat=$tempdir/asl2struct.mat --warp=$warp $APPLYWARP_OPTS
    elif [ ! -z $trans ]; then
        flirt -in $tempdir/$subdir_in/$param -applyxfm -init $tempdir/asl2std.mat -ref $std -out $outdir/std_space/$subdir_out/$parname $APPLYXFM_OPTS
    else
        Log "ERROR: cannot do standard space out - neither warp or transformation matrix found"
    fi
fi

}

OutputMasked() {
    # Function to output images having been masked
    # currently we only do this in native space

    parname=$1
    subdir_out=$2
    maskname=$3

    if [ ! -z $nativeout ]; then
    fslmaths $outdir/native_space/$subdir_out/$parname -mas $tempdir/${maskname}mask $outdir/native_space/$subdir_out/${parname}_masked
    fi
}

Calibrate() {

    param=$1 #parameter as know by basil (and related scripts)
    parname=$2 #parameter name as called in the output directory
    Moval=$3 # Either a number or an image in native space for the calibration scaling
    multiplier=$4 # multiplier that we might use to get parameter into physiological units
    subdir_in=$5 # subdirectory to use for finding the result
    subdir_out=$6 # subdirectory to use for saving the result

    fslmaths  $tempdir/$subdir_in/$param -div $Moval -mul $multiplier $tempdir/$subdir_in/${param}_calib
    Output ${param}_calib ${parname}_calib $subdir_in $subdir_out
}

Report() {
    if [ ! -z $pvexist ]; then
    # generate text reports on parameters - the parameter must have been output first for this to work
    #NB we only do this is the PVE are available (and thus the reqd masks will exist)
    parname=$1
    subdir_out=$2
    masktype=$3

    repval=`fslstats $outdir/native_space/$subdir_out/$parname -k $tempdir/${masktype}mask_pure -m`
    echo $repval > $outdir/native_space/$subdir_out/${parname}_${masktype}_mean.txt
    Log "Mean $parname in $masktype is $repval"

    if [ -f "$tempdir/${masktype}mask_pure_cort."* ]; then
        repval=`fslstats $outdir/native_space/$subdir_out/$parname -k $tempdir/${masktype}mask_pure_cort -m`
        echo $repval > $outdir/native_space/$subdir_out/${parname}_cortical_${masktype}_mean.txt
        Log "Mean $parname in cortical $masktype is $repval"
    fi
    # This nonsense exists because we want the WM mask minus subcortial structures to be
    # referred to as 'cerebral' rather than 'cortical'
    if [ -f "$tempdir/${masktype}mask_pure_cereb."* ]; then
        repval=`fslstats $outdir/native_space/$subdir_out/$parname -k $tempdir/${masktype}mask_pure_cereb -m`
        echo $repval > $outdir/native_space/$subdir_out/${parname}_cerebral_${masktype}_mean.txt
        Log "Mean $parname in cerebral $masktype is $repval"
    fi

    fi
}

Normalise() {
    if [ ! -z $pvexist ]; then
    # also output the perfusion images normalised by the mean mask value - the parameter must have been output first for this to work
    #NB we only do this is the PVE are available (and thus the reqd masks will exist)
    parname=$1
    subdir_out=$2
    masktype=$3
    
    # get normalization from reported value in the output directory
    normval=`cat $outdir/native_space/$subdir_out/${parname}_${masktype}_mean.txt`
    
    if [ ! -z $stdout ]; then
        fslmaths $outdir/std_space/$subdir_out/$parname -div $normval $outdir/std_space/$subdir_out/${parname}_norm
    fi
    if [ ! -z $nativeout ]; then
        fslmaths $outdir/native_space/$subdir_out/$parname -div $normval $outdir/native_space/$subdir_out/${parname}_norm
    fi
    if [ ! -z $structout ]; then
        fslmaths $outdir/struct_space/$subdir_out/$parname -div $normval $outdir/struct_space/$subdir_out/${parname}_norm
    fi

    fi
}

Normalise_var() {
    if [ ! -z $pvexist ]; then
    # normalisaiton for a variance image
    parname=$1
    subdir_out=$2
    masktype=$3

    # get normalization from reported value in the output directory
    normval=`cat $outdir/native_space/$subdir_out/${parname}_${masktype}_mean.txt`
    #need to square the value as we are outputting a variance
    normval=`echo "$normval * $normval" | bc`
    
    
    
    if [ ! -z $stdout ]; then
    fslmaths $outdir/std_space/$subdir_out/${parname}_var -div $normval $outdir/std_space/$subdir_out/${parname}_var_norm
    fi
    if [ ! -z $nativeout ]; then
    fslmaths $outdir/native_space/$subdir_out/${parname}_var -div $normval $outdir/native_space/$subdir_out/${parname}_var_norm
    fi
    if [ ! -z $structout ]; then
    fslmaths $outdir/struct_space/$subdir_out/${parname}_var -div $normval $outdir/struct_space/$subdir_out/${parname}_var_norm
    fi

    fi
}


Registration() {
    Log "Performing registration"
    regbase=$1 #the i/p to the function is the image to use for registration
    transopt=$2 # other options to pass to asl_reg
    distout=$3 # we want to do distortion correction and save in the subdir distout
    
    extraoptions=""
    if [ ! -z $lowstrucflag ]; then
    extraoptions=$extraoptions"-r $tempdir/lowstruc"
    fi
    if [ ! -z $debug ]; then
    extraoptions=$extraoptions" --debug"
    fi
    
    #if [ ! -z $reginit ]; then
    #    extraoptions=$extraoptions" -c $reginit"
    #fi
    
    if [ -z $distout ]; then
    # normal registration
    $asl_reg -i $regbase -o $tempdir -s $tempdir/struc --sbet $tempdir/struc_bet $transopt $extraoptions

    if [ ! -z $trans ]; then
        # compute the transformation needed to standard space if we have the relvant structural to standard transform
        convert_xfm -omat $tempdir/asl2std.mat -concat $trans $tempdir/asl2struct.mat
    fi
        
    else
    # registration for distortion correction
    fmapregstr=""
    if [ ! -z $nofmapreg ]; then
        fmapregstr="--nofmapreg"
    fi
    $asl_reg -i $regbase -o $tempdir/$distout -s $tempdir/struc --sbet $tempdir/struc_bet $transopt $extraoptions --fmap=$tempdir/fmap --fmapmag=$tempdir/fmapmag --fmapmagbrain=$tempdir/fmapmagbrain --pedir=$pedir --echospacing=$echospacing $fmapregstr
    fi


}

Calibration() {
 Log "Calculating M0a - calling ASL_CALIB"
    extraoptions=""
    if [ ! -z $debug ]; then
    extraoptions=$extraoptions" --debug"
    fi

    #if [ ! -z $cref ]; then
    # pass calibration reference image to asl_calib
    #extraoptions=$extraoptions" --cref $tempdir/cref"
    if [ ! -z $senscorr ]; then
        # use a sensitivity iamge from elsewhere
        Log "Sensitivity image $outdir/native_space/sensitivity being loaded into asl_calib"
        extraoptions=$extraoptions" --isen $outdir/native_space/sensitivity"
    fi

    if [ -z $te ]; then
    #by default assume TE is zero
    te=0
    fi

    if [ ! -z $t2star ]; then
    # tell asl_calib to correct for T2* rather than T2
    extraoptions=$extraoptions" --t2star"
    fi

    if [ ! -z $tissref ]; then
    # Specify reference tissue type
    extraoptions=$extraoptions" --tissref $tissref"
    fi

    if [ ! -z $t1csf ]; then
    # supply the T1 of csf
    extraoptions=$extraoptions" --t1r $t1csf"
    fi

    if [ ! -z $t2csf ]; then
    # Supply the T2(*) of CSF
    extraoptions=$extraoptions" --t2r $t2csf"
    fi

    if [ ! -z $t2bl ]; then
    # Supply the T2(*) of blood
    extraoptions=$extraoptions" --t2b $t2bl"
    fi

    if [ ! -z $debug ]; then
    #run asl_calib in debug mode
    extraoptions=$extraoptions" --debug"
    fi

    # setup the main options that we will pass to aslcalib regardless of whether we are auot generating reference mask
    maincaliboptions="--cgain $cgain --te $te --tr $tr"

    if [ -z $csfflag ]; then
    # call asl_calib in normal (auto csf) mode

    # use low res structural for auto generation of csf mask if availible
    # otherwise just use structural image
#    if [ -z $lowstrucflag ]; then
#        usestruc=$tempdir/struc_bet
#        usetrans=$tempdir/asl2struct.mat
#    else
#        usestruc=$tempdir/lowstruc_bet
#        usetrans=$tempdir/asl2lowstruct.mat
#    fi

    usestruc=$tempdir/struc_bet
    usetrans=$tempdir/asl2struct.mat

    if [ ! -z $fasthasrun ]; then
            # we have already run FAST so we can pass the PVE to asl_calib (to save running FAST again)
            if [ -z $tissref ]; then
              extraoptions=$extraoptions" --refpve $tempdir/pvcsf_struct"
            else
              extraoptions=$extraoptions" --refpve $tempdir/pv${tissref}_struct"
        fi
        fi

    $asl_calib -c $tempdir/calib -s $usestruc -t $usetrans -o $outdir/calib --bmask $tempdir/mask --osen $outdir/native_space/sensitivity $maincaliboptions $extraoptions 

    else
    # a manual csf mask has been supplied
    $asl_calib -c $tempdir/calib -m $csf -o $outdir/calib --bmask $tempdir/mask --osen $outdir/native_space/sensitivity $maincaliboptions $extraoptions
    fi
}

Dobasil() {
# inputs: datafile tempdir_subdir initmvn

if [ -z $fast ]; then
    fast=""
else 
    fast="--fast $fast"
fi

Log "Run time basil options:"
Log "$basil_options"
Log "---"

if [ ! -z $3 ]; then
    # we are being supplied with an intial MVN - pass to BASIL
    initmvn="--init $3"
    Log "Initial MVN for BASIL is: $3"
else
    initmvn=""
fi

# Run Basil with dilated mask to ensure all required voxels are fitted
$basil -i $1 -o $2/basil -m $tempdir/mask_dil -@ $2/basil_options.txt $fast $initmvn $basil_options


# work out which is the final step from BASIL
finalstep=`ls -d $2/basil/step? | sed -n '$ p'`
Log "Using BASIL step $finalstep"

# extract images from BASIL results (and throw away values below zero)
fslmaths ${finalstep}/mean_ftiss -thr 0 $2/ftiss
if [ ! -z $senscorr ]; then
# sensitivity correction
    fslmaths $2/ftiss -div $outdir/native_space/sensitivity $2/ftiss
fi

if [ -z $fixbat ]; then
    fslmaths ${finalstep}/mean_delttiss -thr 0 $2/delttiss
fi
if [ -z $artoff ]; then
    fslmaths ${finalstep}/mean_fblood -thr 0 $2/fblood
    if [ ! -z $senscorr ]; then
# sensitivity correction
    fslmaths $2/fblood -div $outdir/native_space/sensitivity $2/fblood
    fi
fi

#Partial volume correction - sort out basil results when PV corrected
if [ ! -z $basil_pvcorr ]; then
    fslmaths ${finalstep}/mean_fwm -thr 0 $2/ftisswm
    if [ ! -z $senscorr ]; then
        # sensitivity correction
    fslmaths $2/ftisswm -div $outdir/native_space/sensitivity $2/ftisswm
    fi
    if [ -z $fixbat ]; then
    fslmaths ${finalstep}/mean_deltwm -thr 0 $2/deltwm
    fi

fi

if [ ! -z $varout ]; then
#get varainces out of finalMVN
    fabber_var -d ${finalstep} -m $tempdir/mask_dil
# do correction of negative values
    fslmaths ${finalstep}/var_ftiss -bin -add 1 -uthr 1 -mul 1e12 -add ${finalstep}/var_ftiss $2/var_ftiss
    if [ ! -z $senscorr ]; then
# sensitivity correction
    fslmaths $2/var_ftiss -div $outdir/native_space/sensitivity -div $outdir/native_space/sensitivity $2/var_ftiss
fi
    if [ -z $fixbat ]; then
    fslmaths ${finalstep}/var_delttiss -bin -add 1 -uthr 1 -mul 1e12 -add ${finalstep}/var_delttiss $2/var_delttiss
    fi

    #Partial volume correction - output variances of PVC results
    if [ ! -z $basil_pvcorr ]; then
        fslmaths ${finalstep}/var_fwm -bin -add 1 -uthr 1 -mul 1e12 -add ${finalstep}/var_fwm $2/var_ftisswm
        if [ ! -z $senscorr ]; then
            # sensitivity correction
            fslmaths $2/var_ftisswm -div $outdir/native_space/sensitivity -div $outdir/native_space/sensitivity $2/var_ftisswm
        fi
        if [ -z $fixbat ]; then
            fslmaths ${finalstep}/var_deltwm -bin -add 1 -uthr 1 -mul 1e12 -add ${finalstep}/var_deltwm $2/var_deltwm
        fi
    fi
fi

#copy the final MVN to the temp directory for future use
imcp ${finalstep}/finalMVN $2/finalMVN
cp ${finalstep}/paramnames.txt $2/paramnames.txt


}

Dooutput() {

# Do all the outputs - using the supplied subdirectiory of the results
subdir_in=$1
subdir_out=$2
cmethod_output=$3

# perfusion
Output ftiss perfusion $subdir_in $subdir_out
Report perfusion $subdir_out gm
Normalise perfusion $subdir_out gm

# arrival
if [ -z $fixbat ]; then
    Output delttiss arrival $subdir_in $subdir_out
    Report arrival $subdir_out gm
fi
# aBV
if [ -z $artoff ]; then
    Output fblood aCBV $subdir_in $subdir_out
fi

# white matter values
if [ $subdir_in = "pvcorr" ]; then
    Output ftisswm perfusion_wm $subdir_in $subdir_out
    Report perfusion_wm $subdir_out wm
    Normalise perfusion_wm $subdir_out wm
    if [ -z $fixbat ]; then
        Output deltwm arrival_wm $subdir_in $subdir_out
        Report arrival_wm $subdir_out wm
    fi
    
else 
    Report perfusion $subdir_out wm
    if [ -z $fixbat ]; then
        Report arrival $subdir_out wm
    fi
fi

# Masked results (PVcorr)
if [ $subdir_in = "pvcorr" ]; then
    OutputMasked perfusion $subdir_out gm
    OutputMasked perfusion_wm $subdir_out wm
    if [ -z $fixbat ]; then
        OutputMasked arrival $subdir_out gm
        OutputMasked arrival_wm $subdir_out wm
    fi
fi

# Optionally provide variance results
if [ ! -z $varout ]; then
    Output var_ftiss perfusion_var $subdir_in $subdir_out
    Normalise_var perfusion $subdir_out gm
    if [ -z $fixbat ]; then
        Output var_delttiss arrival_var $subdir_in $subdir_out
    fi
        
    # PVcorr
    if [ $subdir_in = "pvcorr" ]; then
        Output var_ftisswm perfusion_wm_var $subdir_in $subdir_out
        Normalise_var perfusion_wm $subdir_out wm
        if [ -z $fixbat ]; then
            Output var_deltwm arrival_wm_var $subdir_in $subdir_out
        fi
    fi
fi

# calibrated results
if [ ! -z $calibflag ]; then
    if [ $cmethod_output = 'single' ]; then
        malpha=`echo "$Mo * $alpha" | bc` #include the inversion efficiency when we do the final calibration
    elif [ $cmethod_output = 'voxel' ]; then
        # Sensitivity map is re-calculated in PVC mode, so 'undo' the old correction and re-do the new one
        if [ ! -z $senscorr ]; then
            if [ $subdir_in = "pvcorr" ]; then
                fslmaths $outdir/calib/M0 -mul $outdir/native_space/sensitivity_non_pvcorr $outdir/calib/M0
            fi
            fslmaths $outdir/calib/M0 -div $outdir/native_space/sensitivity $outdir/calib/M0
        fi
        fslmaths $outdir/calib/M0 -mul $alpha $tempdir/malpha
        malpha=$tempdir/malpha
    fi

    Calibrate ftiss perfusion $malpha 6000 $subdir_in $subdir_out
    Report perfusion_calib $subdir_out gm
   
    if [ $subdir_in = "pvcorr" ]; then
        OutputMasked perfusion_calib $subdir_out gm
        Calibrate ftisswm perfusion_wm $malpha 6000 $subdir_in $subdir_out
        Report perfusion_wm_calib $subdir_out wm
        OutputMasked perfusion_wm_calib $subdir_out wm
    else
        Report perfusion_calib $subdir_out wm
    fi

    if [ ! -z $varout ]; then
        if [ $cmethod_output = 'single' ]; then
            Mosq=`echo "$Mo * $Mo * $alpha * $alpha" | bc` #include the inversion efficiency when we do the final calibration
        elif [ $cmethod_output = 'voxel' ]; then
            fslmaths $outdir/calib/M0 -mul $outdir/calib/M0 -mul $alpha -mul $alpha $tempdir/mosq
            Mosq=$tempdir/mosq
        fi
        
        Calibrate var_ftiss perfusion_var $Mosq 36000000 $subdir_in $subdir_out

        if [ $subdir_in = "pvcorr" ]; then
            Calibrate var_ftisswm perfusion_wm_var $Mosq 36000000 $subdir_in $subdir_out
        fi
    fi

    if [ -z $artoff ];then
        # output aCBV as a percentage
        Calibrate fblood aCBV $malpha 100 $subdir_in $subdir_out
    fi
fi

# advanced output
if [ ! -z $advout ]; then
    if [ ! -d  $outdir/advanced/$subdir_out ]; then mkdir $outdir/advanced/$subdir_out; fi
    imcp $tempdir/$subdir_in/finalMVN $outdir/advanced/$subdir_out/finalMVN
   cp $tempdir/$subdir_in/paramnames.txt $outdir/advanced/$subdir_out/paramnames.txt
fi

}

Log() {
    # Output text to log and to standard output stream
    echo $1
    if [ ! -z $log ]; then
        echo $1 >> $log
    fi
}

ApplyWarps() {
    ### Apply correction warps to the data, including a motion correction transformation if
    ### applicable, and a nonlinear warp field

    corrwarp=$1
    jacobian=$2

    Log "  Applying combined warps to original data:"
    Log "   - Warp: $corrwarp"
    Log "   - Jacobian: $jacobian"

    Log "  Applying warps to ASL data"
    # Add motion correction if required
    appremat=""
    if [ ! -z $moco ]; then
        moco_xfm=$3
        Log "   - Moco: $moco_xfm"
        appremat="--premat=$moco_xfm"
    fi
    applywarp -i $tempdir/asldata_orig -r $tempdir/nativeref -o $tempdir/asldata $appremat -w $corrwarp --rel $APPLYWARP_OPTS

    # Correct for quantitative signal magnitude as volume has been locally scaled by the Jacobian
    fslmaths $tempdir/asldata -mul $jacobian $tempdir/asldata
    
    # Repeat post correction as may be used as regfrom
    fslmaths $tempdir/asldata -Tmean $tempdir/meanasl

    # Apply warp to calibration images
    if [ ! -z $calib ]; then
        Log "  Applying warps to calibration data"
        appremat=""
        if [ ! -z $moco ]; then
            calib_xfm=$4
            Log "   - Calib->ASL: $calib_xfm"
            appremat="--premat=$calib_xfm"
        fi

        applywarp -i $tempdir/calib_orig -r $tempdir/nativeref -o $tempdir/calib $appremat -w $corrwarp --rel $APPLYWARP_OPTS
        fslmaths $tempdir/calib -mul $jacobian $tempdir/calib
        if [ ! -z $cref ]; then
            applywarp -i $tempdir/cref_orig -r $tempdir/nativeref -o $tempdir/cref $appremat -w $corrwarp --rel $APPLYWARP_OPTS
            fslmaths $tempdir/cref -mul $jacobian $tempdir/cref
        fi
        if [ ! -z $cblip ]; then
            applywarp -i $tempdir/cblip_orig -r $tempdir/nativeref -o $tempdir/cblip $appremat -w $corrwarp --rel $APPLYWARP_OPTS
            fslmaths $tempdir/cblip -mul $jacobian $tempdir/cblip
        fi
    fi
}

# deal with options

if [ -z $1 ]; then
    Usage
    exit 1
elif [ $1 = "--more" ]; then
    Usage
    Usage_extended
    exit 1
fi

#basil=basil
#asl_reg=asl_reg
#asl_calib=asl_calib

basil=basil
asl_reg=asl_reg
asl_calib=asl_calib

# defaults that (boolean) command line options can overide
spatial=1; # we always use spatial priors, use --spatial=off to turn off.
finalreg=1; # we always try to refine the registration using the perfusion image (Registration 2)
fixbolus=undef; # the default for this is determined by other command line parameters.
edgecorr=1; # we do edge correction with the voxelwise calibration method by default
regfrom_bet_thresh=0.2; # Extraction threshold fraction to use with BET - by default conservative to avoid erosion of brain

# Declare arrays for additional ROI atlases
declare -a region_analysis_atlases
declare -a region_analysis_atlases_labels

# parse command line here
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
        takeargs=1;;
    -i) inflag=1 infile=$argument #input/data file
        takeargs=1;;
    --iaf) iaf=$argument # input asl format (asl_file syntax)
           takeargs=1;;
    --ibf) ibf=$argument # input block format (asl_file syntax)
           takeargs=1;;
    --rpts) rpts=$argument # the number of repeats at each TI (when --ibf=tis) - to be passed to asl_file
        takeargs=1;;
    -c) calibflag=1 calib=$argument #calibration image
        takeargs=1;;
    --wp) whitepaper=1; # operate in 'white paper' mode
        ;;
    -s) strucflag=1  #NOTE: use strucflag to determin if structural information (in any form) has been provided
        struc=$argument # strucutral image (not BETed)
        takeargs=1;;
    --sbrain) strucbet=$argument # user supplied BETed structural
          takeargs=1;;
    --fslanat) strucflag=1;
           fslanat=$argument # user supplied fslanat dir (overrides the structural image inputs)
           takeargs=1;;
    --fastsrc) strucflag=1;
           fastsrc=$argument # user supplied FAST results (this is the stub of the filenames)
           takeargs=1;;
    -r) lowstrucflag=1 lowstruc=$argument #low resolution structural image
        takeargs=1;;
    -t) transflag=1;
        trans=$argument;
        stdout=1;
        takeargs=1;;
    -S) stdflag=1 std=$argument
        takeargs=1;;
    --warp) transflag=1;
        warp=$argument; #structural to standard warp
        stdout=1;
        takargs=1;;
    -m) mask=$argument
        takeargs=1;;
    --mc) isbool=1; # do motion correction using mcflirt
          boolarg=moco;
          ;;
    --tis) tis=$argument
        takeargs=1;;
    --plds) plds=$argument
        takeargs=1;;
    --tiimg) tiimg=$argument
        takeargs=1;;
    --bolus) boluset=1 boluslen=$argument
        takeargs=1;;
    --bat) bat=$argument
        takeargs=1;;
    --batsd) batsd=$argument
        takeargs=1;;
    --slicedt) slicedt=$argument
        takeargs=1;;
    --sliceband) sliceband=$argument
        takeargs=1;;
    --fa) fa=$argument
        takeargs=1;;
    --t1) t1set=$argument # the T1 of tissue to be used in kinetic model
        takeargs=1;;
    --t1b) t1bset=$argument
           takeargs=1;;
    --t1im) t1im=$argument # A T1 (of tissue) image
        takeargs=1;;
    --batim) batim=$argument # A BAT (seconds) image
        takeargs=1;;
    --noiseprior) noiseprior=1
        ;;
    --noisesd) noisesd=$argument
        takeargs=1;;
    --cmethod) cmethod=$argument
        takeargs=1;;
    --tissref) tissref=$argument
        takeargs=1;;
    --edgecorr) isbool=1;
            boolarg=edgecorr;
            ;;
    --csf) csfflag=1 csf=$argument
        takeargs=1;;
    --cref) cref=$argument; senscorr=1;
        takeargs=1;;
    --isen) isen=$argument; senscorr=1;
        takeargs=1;;
    --senscorr) needseg=1; senscorr=1; #disbale 
        ;;
    --M0) M0=$argument; calibflag=1;
        takeargs=1;;
    --t1t) t1tset=$argument; #the T1 of tissue to be used in calibration
        takeargs=1;;
    --tr) tr=$argument
        takeargs=1;;
    --te) te=$argument #supply the echo time for calibration correction for T2
        takeargs=1;;
    --t2star) t2star=1 #do calibration with T2* rather than T2
        ;;
    --t1csf) t1csf=$argument #custom T1 for CSF
        takeargs=1;;
    --t2csf) t2csf=$argument #custom T2 for CSF
        takeargs=1;;
    --t2bl) t2bl=$argument #custom T2 of blood
        takeargs=1;;
    --regfrom) regfromflag=1 regfrom=$argument
        takeargs=1;;
    --regfrom-method) regfrom_method=$argument
        takeargs=1;;
    --reg-init-bbr) reg_init_bbr=1
        ;;
    --regfrom-bet-thresh) regfrom_bet_thresh=$argument
        takeargs=1;;
    --finalreg) isbool=1; # To turn off the final registration step that uses the perfusion image itself (with BBR) - helpful is the perfusion image is poor and thus not a good basis for registration
        boolarg=finalreg;
        ;;
    --asl2struc) asl2struc=$argument
        takeargs=1;;
    --cgain) cgain=$argument
        takeargs=1;;
    --alpha) alpha=$argument
        takeargs=1;;
    --zblur) zblur=1
        ;;
    --structout) structout=1
        ;;
    --advout) advout=1
        ;;
    --spatial) #spatial=1
        isbool=1;
        boolarg=spatial;
        ;;
    --infert1) infert1=1
        ;;
    --artoff) artoff=1
        ;;
    --artonly) artonly=1
        ;;
    --artsupp) artoff=1 #this does same job as --artoff, but is explicitly linked to vascular crushers in the data
        ;;
    --fixbat) fixbat=1
        ;;
    --fixbolus) isbool=1
        boolarg=fixbolus;
        ;;
    --casl) casl=1
        ;;
    --disp) disp=$argument
        takeargs=1;;
    --exch) exch=$argument
        takeargs=1;;
    --pvcorr) pvcorr=1
        ;;
    --pvgm) pvgm=$argument
        takeargs=1;;
    --pvwm) pvwm=$argument
        takeargs=1;;
    --fulldata) fulldata=1
        ;;
    --elen) epoch=1 elen=$argument
        takeargs=1;;
    --eol) eol=$argument
           takeargs=1;;
    --fast) fast=2 #do a one step analysis (NB generally a good idea to have spatial on for this option)
        ;;

    # fieldmap/distorition correction related arguments
    --fmap) fmap=$argument
        takeargs=1;;
    --fmapmag) fmapmag=$argument
        takeargs=1;;
    --fmapmagbrain) fmapmagbrain=$argument
        takeargs=1;;
    --echospacing) echospacing=$argument
        takeargs=1;;
    --pedir) pedir=$argument
        takeargs=1;;
    --nofmapreg) isbool=1
        boolarg=nofmapreg
        ;;
    --gdcwarp) gdcwarp=$argument #User-supplied gradient distortion correction warp
        takeargs=1;; 
    --cblip) cblip=$argument #blip reversed calibration image
        takeargs=1;; 
    --model-options) model_options=$argument
        takeargs=1;;
    --region-analysis) region_analysis=1
        ;;
    --region-analysis-save-rois) region_analysis_save_rois=1
        ;;
    --region-analysis-atlas) region_analysis_atlases+=( "$argument" )
        takeargs=1;;
    --region-analysis-atlas-labels) region_analysis_atlases_labels+=( "$argument" )
        takeargs=1;;
    --region-analysis-psf) region_analysis_psf+=( "$argument" )
        takeargs=1;;
    --gm-thresh) gm_thresh=$argument
        takeargs=1;;
    --wm-thresh) wm_thresh=$argument
        takeargs=1;;
    --qc-output) qc_output=1
        ;;
    --verbose) verbose=$argument
        takeargs=1;;
    --debug) debug=1
        ;; 
    --devel) devel=1
        ;;
    --version) Version
        ;;
    *)  #Usage
        echo "Error! Unrecognised option on command line: $option"
        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 to move on to next parameter
    shift
done

nativeout=1 #we always keep the native space images!!
varout=1 # we always output the variance images

if [ -z $verbose ]; then
    verbose=1
fi

# deal with the temporary directory
tmpbase=`tmpnam`
tempdir=${tmpbase}_ox_asl
mkdir $tempdir

echo "#FABBER options created by Oxford_asl" > $tempdir/basil_options.txt

# deal with default output format
if [ ! -z $transflag ]; then
   #if transformation matrix included then output in std space
    stdout=1;
    if [ -z $strucflag ]; then
    echo "ERROR: Structural image is required along with transformation matrix to output results in standard space"
    exit 1
    fi
fi

if [ ! -z $strucflag ]; then
    # else-if strucutral image included output in structural space
    structout=1;
fi

# set the output directory here if not specified
if [ -z $outflag ]; then
    Log "Ouput being placed in input directory"
    outdir=`pwd`;
fi

# Start by looking for the output directory (and create if need be)
if [ ! -d $outdir ]; then
  mkdir $outdir;
  Log "Created output directory $outdir"
fi

log=$outdir/logfile
rm "$log"

Log "oxford_asl $ARGS"
Log "OXFORD_ASL - running"
Log "Version: v4.0.29-5-g09a5cba5-dirty Thu Jul 13 14:20:21 2023"

#check required inputs are present
if [ -z $inflag ]; then
    echo "ERROR: no input file specified"
    exit 1
else
    if [ `imtest $infile` -eq 0 ]; then
    echo "ERROR: $infile is not an image/has not been found"
    exit 1
    fi
fi
Log "Input file: $infile"

if [ ! -z $strucflag ]; then
    if [ ! -z $struc ]; then
    if [ `imtest $struc` -eq 0 ]; then
        echo "ERROR: $struc is not an image/has not been found"
        exit 1
    fi
    Log "Structural image: $struc"
    fi
    
fi




if [ ! -z $lowstruc ]; then
    if [ `imtest $lowstruc` -eq 0 ]; then
    echo "ERROR: $lowstruc is not an image/has not been found"
    exit 1
    fi
    Log "Low res. structural image: $lowstruc"
fi


# Setup option outputs - main subdirectories and anything that would be common to all epochs
if [ ! -z $nativeout ] && [ ! -d $outdir/native_space ]; then
    Log "Saving results in natve (ASL aquisition) space to $outdir/native_space"
    mkdir $outdir/native_space
fi
if [ ! -z $structout ] && [ ! -d $outdir/struct_space ]; then
    Log "Saving results in structural space to $outdir/struct_space"
    mkdir $outdir/struct_space
fi
if [ ! -z $advout ] && [ ! -d $outdir/advanced ]; then
Log "Saving advanced outputs"
    mkdir $outdir/advanced
fi

### Command line parameter interactions
# white paper mode - this overrieds defaults, but can be overwritten (below) by command line specification of individual parameters
if [ ! -z $whitepaper ]; then
    Log "Analysis in white paper mode"
    t1set=1.65
    cmethod=voxel
    bat=0 # the white paper model ignores BAT
    # note that other things related to the model will be set by single TI mode below
fi

# bolus duration inference
# if we are doing CASL then fix the bolus duration, except where the user has explicitly told us otherwise
if [ $fixbolus = "undef" ]; then
    # fixbolus is to take its default value
    if [ ! -z $casl ]; then
    fixbolus=1;
    else
    fixbolus="";
    fi
fi

### End of Command line parameter interactions

# general pre-processing
Log "Pre-processing"
if [ ! -z $fslanat ]; then
    # we are being supplied with an fslanat directory
    # copy over the structural and brain extracted strucutral so that we can use them
    if ls $fslanat/T1_biascorr.* 1> /dev/null 2>&1; then
    Log "Using structural images (bias corrected) from fsl_anat: $fslanat"
    imcp $fslanat/T1_biascorr $tempdir/struc
    imcp $fslanat/T1_biascorr_brain $tempdir/struc_bet
    else
    Log "Using structural images from fsl_anat: $fslanat"
    imcp $fslanat/T1 $tempdir/struc
    imcp $fslanat/T1_biascorr_brain $tempdir/struc_bet
    fi
    chmod u+w $tempdir/struc.* $tempdir/struc_bet.*
    
    # use the registration to standard space - if it is there
    if [ -e $fslanat/T1_to_MNI_nonlin_coeff.nii.gz ]; then
        stdout=1; transflag=1;
    warp=$fslanat/T1_to_MNI_nonlin_coeff.nii.gz
    elif [ -e $fslanat/T1_to_MNI_lin.mat ]; then
        stdout=1; transflag=1;
    trans=$fslanat/T1_to_MNI_lin.mat
    fi
    
elif [ ! -z $struc ]; then
    fslmaths $struc $tempdir/struc
    Log "Structural image is: $struc"
    if [ -z $strucbet ]; then
       #bet the structural for calibration and registration
    bet $struc $tempdir/struc_bet
    Log "BET on structural image"
    else
    fslmaths $strucbet $tempdir/struc_bet
    fi
fi

if [ ! -z $transflag ]; then
    # deal with Standard brain image
    if [ -z $stdflag ]; then
    std=${FSLDIR}/data/standard/MNI152_T1_2mm
    fi
    Log "Standard brain is: $std"

    if [ ! -z $warp ]; then
    Log "Structural to standard transformation warp: $warp"
    elif [ ! -z $trans ]; then
    Log "Structural to standard transformation matrix: $trans"
    fi
fi

if [ ! -z $lowstrucflag ]; then
    fslmaths $lowstruc $tempdir/lowstruc
## bet the low res. struc (if present)
#    bet $lowstruc $tempdir/lowstruc_bet
#    Log "BET on low res. structural image"
fi

# standard pre-processing of calibration image
if [ ! -z $calib ]; then
    tsize=`fslinfo $calib | grep "^dim4" | sed 's:dim4[ ]*::'`
    if [ $tsize -gt 1 ]; then
        #cut - throw away first volume
    Log "Removing first volume of calibration time series - to ensure data is in steady state"
    tsize=`expr $tsize - 1`
    fslroi $calib $tempdir/calib 1 $tsize

    if [ ! -z $moco ]; then
        #motion correction
        mcflirt -in $tempdir/calib -o $tempdir/calib
    fi
    # take the mean
    fslmaths $tempdir/calib -Tmean $tempdir/calib
    else
    fslmaths $calib $tempdir/calib
    fi
    # Save original so we can reapply motion correction/warps to original data
    imcp $tempdir/calib $tempdir/calib_orig
fi

# same thing for cref image
if [ ! -z $cref ]; then
    tsize=`fslinfo $cref | grep "^dim4" | sed 's:dim4[ ]*::'`
    if [ $tsize -gt 1 ]; then
        #cut - throw away first volume
    Log "Removing first volume of calibration reference time series - to ensure data is in steady state"
    tsize=`expr $tsize - 1`
    fslroi $cref $tempdir/cref 1 $tsize

    if [ ! -z $moco ]; then
        #motion correction
        mcflirt -in $tempdir/cref -o $tempdir/cref
    fi
    # take the mean
    fslmaths $tempdir/cref -Tmean $tempdir/cref
    else
    fslmaths $cref $tempdir/cref
    fi
    # Save original so we can reapply motion correction/warps to original data
    imcp $tempdir/cref $tempdir/cref_orig
fi

# same thing for cblip image
if [ ! -z $cblip ]; then
    tsize=`fslinfo $cblip | grep "^dim4" | sed 's:dim4[ ]*::'`
    if [ $tsize -gt 1 ]; then
        #cut - throw away first volume
    Log "Removing first volume of blipped calibration time series - to ensure data is in steady state"
    tsize=`expr $tsize - 1`
    fslroi $cblip $tempdir/cblip 1 $tsize

    if [ ! -z $moco ]; then
        #motion correction
        mcflirt -in $tempdir/cblip -o $tempdir/cblip
    fi
    # take the mean
    fslmaths $tempdir/cblip -Tmean $tempdir/cblip
    else
    fslmaths $cblip $tempdir/cblip
    fi
    # NOTE cblip is still wholehead, we dont need it masked
    # Save original so we can reapply motion correction/warps to original data
    imcp $tempdir/cblip $tempdir/cblip_orig
fi

# read in ASL data
imcp $infile $tempdir/asldata # this is the MAIN data that we will reflect any corrections applied
chmod u+w $tempdir/asldata.*
# take a copy that will not be subject to any subsequent corrections
imcp $tempdir/asldata $tempdir/asldata_orig

# Take the mean of the ASL data to use as a reference image in applywarp. 
# This is required because in FSL6, applywarp duplicates the output for every
# volume of the reference image, so it is imperative that it is a 3D image!
# Note that in applywarp the only function of the reference image is to define
# the output space, so this never needs to be changed as our 'native' space is
# always the same space as the input ASL data
fslmaths $tempdir/asldata -Tmean $tempdir/nativeref

### Motion Correction (main)
# note motion correction within calibration data is done above
if [ ! -z $moco ]; then
    Log "Motion Correction"
    if [ ! -z $calib ]; then
    # we use the calibration image as our reference for motion correction - since this will be most consistent if the data has a range of different TIs and background suppression etc
    # this also removes motion effects between asldata and calibration image
    Log "Motion correction to calibration image"
    mcflirt -in $tempdir/asldata -out $tempdir/asldata -r $tempdir/calib -mats -plots

    # to reduce interpolation of the ASL data change the transformations so that we end up in the space of the central volume of asldata
    # Extract the middle transformation
    middlemat=`ls $tempdir/asldata.mat/MAT_* | awk '{ lines[NR]=$0; } END { print lines[int(NR/2)+1] }'`
    convert_xfm -omat $tempdir/calib2asl.mat -inverse $middlemat
    # convert all the transformations
    for mat in `ls $tempdir/asldata.mat/MAT_*`; do
        convert_xfm -omat $mat -concat $tempdir/calib2asl.mat $mat
    done
    # Concatenate the motion matrices for use with applywarp
    cat $tempdir/asldata.mat/MAT* > $tempdir/asldata.cat
    applywarp --ref=$tempdir/nativeref --in=$tempdir/asldata_orig --out=$tempdir/asldata --premat=$tempdir/asldata.cat $APPLYWARP_OPTS
    # Convert all calibration images to align with asldata
    flirt -in $tempdir/calib -out $tempdir/calib -ref $tempdir/nativeref -init $tempdir/calib2asl.mat -applyxfm $APPLYXFM_OPTS
    if [ ! -z $cref ]; then
        flirt -in $tempdir/cref -out $tempdir/cref -ref $tempdir/nativeref -init $tempdir/calib2asl.mat -applyxfm $APPLYXFM_OPTS
    fi
    if [ ! -z $cblip ]; then
        flirt -in $tempdir/cblip -out $tempdir/cblip -ref $tempdir/nativeref -init $tempdir/calib2asl.mat -applyxfm $APPLYXFM_OPTS
    fi
    else
    Log "Motion correction to middle volume of ASL data"
    mcflirt -in $tempdir/asldata -out $tempdir/asldata -mats -plots
    # Save concatenated motion matrices for later use with applywarp if required
    cat $tempdir/asldata.mat/MAT* > $tempdir/asldata.cat
    fi
fi

# Take mean of the asl data as we might need this in the next step, and subsequent steps
fslmaths $tempdir/asldata -Tmean $tempdir/meanasl

if [ ! -z $gdcwarp ]; then
    fslinfo $gdcwarp >/dev/null 2>/dev/null
    if [ $? -eq 0 ]; then
        ### User supplied distortion correction as (additional) warp
        ### We apply this now, but later if additional distortion correction warps are derived, we will
        ### combine them together and re-apply to raw data to avoid multiple interpolations

        Log "Initial application of GDC warp"
        mkdir $tempdir/distcorr

        # Convert the supplied warp image to extract the jacbbian (as parts)
        # (meanasl here is just a reference). Using relative warp convention
        convertwarp -r $tempdir/nativeref -o $tempdir/asldist_userwarp -w $gdcwarp --rel -j $tempdir/distcorr/jacobian_parts --constrainj

        # Calculation of the jacobian for the warp - method suggested in:
        # https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=FSL;d3fee1e5.0908
        fnirtfileutils -i $tempdir/asldist_userwarp -f spline -o $tempdir/asldist_userwarp_coef
        fnirtfileutils -i $tempdir/asldist_userwarp_coef -j $tempdir/distcorr/jacobian
        fslcpgeom $tempdir/asldist_userwarp $tempdir/distcorr/jacobian -d

        ApplyWarps $tempdir/asldist_userwarp $tempdir/distcorr/jacobian $tempdir/asldata.cat $tempdir/calib2asl.mat
    else
        Log "WARNING: GDC warp file $gdcwarp not found - will not be used"
        gdcwarp=""
    fi
fi

### Bolus duration(s)
if [ -z $boluset ]; then
    boluslen=1.8; # use the WP default
fi

taucount=0
tauslist=""
thetaus=`echo $boluslen | sed 's:,: :g'`
for tau in $thetaus; do
    taucount=`expr ${taucount} + 1`
    tauslist=`echo $tauslist --tau${taucount}=$tau`
done

if [ $taucount -eq 1 ]; then
    tauslist="--tau=$tau" #single univerisal bolus duration
    Log "Bolus duration: $tau"
else
    Log "Bolus duration list: $tauslist"
fi
### End of Bolus duration(s)

### deal with TIs (as specified on the command line)
if [ ! -z $tiimg ]; then
    # Cannot easily define alltis array - therefore cannot do epochwise analysis
    if [ ! -z $epoch ]; then
        echo "ERROR: cannot do epochwise analysis with a voxelwise TI image"
        exit 1
    fi
    if [ ! -z $tis ]; then
        echo "ERROR: Cannot specify a TI image as well as a list of TIs"
        exit 1
    fi
    if [ ! -z $slicedt ]; then
        echo "ERROR: Cannot specify slice timing correction with a voxelwise TI image"
        exit 1
    fi
    Log "Using voxelwise TI image: $tiimg"
    ntis=`fslinfo $tiimg | grep "^dim4" | sed 's:dim4[ ]*::' |tr -d [:blank:]`
    tislist=" --tiimg=$tiimg"
    Log "Number of TIs: $ntis"
else
    if [ ! -z $tis ] && [ ! -z $plds ]; then
        echo "ERROR: Cannot specify both TIs and PLDs"
        exit 1
    fi
    ticount=0
    tislist=""
    thetis=`echo $tis $plds | sed 's:,: :g'`
    tausarray=($thetaus)
    #echo $thetis
    for ti in $thetis; do
        ticount=`expr ${ticount} + 1`
        if [ ! -z $plds ]; then
            if [ $taucount -eq 1 ]; then
                # Single univerisal labelling duration
                ti=`echo $ti + $tau |bc`
            else
                tauidx=`expr ${ticount} - 1`
                ti=`echo $ti + ${tausarray[$tauidx]} |bc`
            fi
        fi
        tislist=`echo $tislist --ti${ticount}=$ti`
        alltis[`expr ${ticount} - 1`]=$ti
    done
    Log "Number of TIs in list: $ticount"
    Log "TIs list: $tislist"
    if [ $taucount -ne 1 ] && [ $taucount -ne $ticount ]; then
        echo "ERROR: Number of bolus durations specified does not match the number of TIs - this is not possible for multiple bolus duration processing"
        exit 1
    fi
    ntis=$ticount
fi
### End of: Deal with the TIs

### Label-control subtraction (we repeat subtraction after doing distortion correction - when applicable)
# We get data into correct block format for BASIL here
if [ -z $iaf ]; then
    # DEFAULT input ASL format is 'diff' (label-control difference data) - maintains backward compatibility
    iaf="diff"
fi
Log "Input ASL format is: $iaf"
if [ -z $ibf ]; then
    # DEFAULT input block format is 'rpt' (data is as typically acquired) - maintains backward compatibility
    ibf="rpt"
fi
Log "Input block format is: $ibf"

afrptstr=""
if [ $ibf == "tis" ]; then
    if [ ! -z $rpts ]; then
    # the repeats at each TI have been specified
    afrptstr="--rpts=$rpts"

    if [ ! -z $epoch ]; then
        echo "ERROR: cannot do epochwise analysis when there are not the same number of repeats for each TI"
        exit 1
    fi
    fi
fi

if [ $iaf = 'diff' ]; then
    # make sure the block format is correct for BASIL
    asl_file --data=$tempdir/asldata --ntis=$ntis --ibf=$ibf $afrptstr --iaf=$iaf --obf=tis --out=$tempdir/diffdata --mean=$tempdir/diffdata_mean
    Log "Label-control difference data provided: $infile"
else
    # create label-control difference data using asl_file - this gets it into the right block form for BASIL (blocks of TIs)
    asl_file --data=$tempdir/asldata --ntis=$ntis --ibf=$ibf $afrptstr --iaf=$iaf --obf=tis --diff --out=$tempdir/diffdata --mean=$tempdir/diffdata_mean
    
    # pull out (label &) control images (NB these will have suffix _even or _odd)  *TODO asl_file should now be able to do suffix of _label and _control, which would be more helpful
    # NOT enabled at the moment - since this is only specifically useful for extensing the calibration options - TODO 
    #asl_file --data=$tempdir/asldata --ntis=$ntis --ibf=$ibf --iaf=$iaf --obf=tis --split --out=$tempdir/asldata --mean=$tempdir/asldata_mean

fi

# Generate a perfusion-weighted image by taking the mean over all TIs of the differenced data
fslmaths $tempdir/diffdata_mean -Tmean $tempdir/pwi

### End of : Label-control subtraction

### Establish the number of repeats in data - query the diffdata file (that will contain all repeats)
tpoints=`fslinfo $tempdir/diffdata | grep "^dim4" | sed 's:dim4[ ]*::'`

Log "Number of inversion times: $ntis"
Log "Number of timepoints in data: $tpoints"

if [ -z $afrptstr ]; then
    repeats=`expr $tpoints / $ntis`
    Log "Number of repeats in data: $repeats"
else
    repeats=0 #variable number of repeats at each TI - have to handle specially
    ### deal with repeats (as specified on the command line)
    count=0
    rptslist=""
    therpts=`echo $rpts | sed 's:,: :g'`
    for rp in $therpts; do
    count=`expr ${count} + 1`
    rptslist=`echo $rptslist --rpt${count}=$rp`
    done
    Log "RPTs list: $rptslist"
fi
### End of: Establish number of repeats

# Single or Multi TI setup
if [ $ntis -lt 2 ]; then
# single TI data - dont average send to BASIL as-is (helps with noise estimation)
#    Log "Single TI data to be passed to BASIL"
    #datafile=$tempdir/diffdata
    # OPERATE in 'single TI mode'
    artoff=1; fixbolus=1; #fixbat=1; #fast=1; 
    Log "Operating in Single TI mode"
    singleti=1;
    
#elif [ $repeats -gt 1 ]; then
 #   if [ -z $fulldata ]; then
        # take the mean over the TIs for faster analysis
#    Log " Multi TI data, mean is being taken at each TI to pass to BASIL"
#    datafile=$tempdir/diffdata_mean
#    repeats=1
 #   else
#    Log "Multi TI data, all data passed to BASIL for analysis"
 #       datafile=$tempdir/diffdata
 #   fi
#else
    # if there is only one repeat it just gets passed stright to BASIL
#    echo " Multi-TI data (single measurment at each TI) to be passed to BASIL" >> $log
    #    datafile=$tempdir/diffdata
fi

# finish filling the alltis array - this is for epochwise analysis - this only works when there are the same number of repeats at each TI
# Epochwise analysis presume that the data was acquired in blocks of all TIs (even if that is not how is supplied to oxford_asl)
# Otherwise epoch analysis doens't make sense
# NB: we use the tislist and numerb of repeats for BASIl for the conventional perfusion image
for ((r=1; r<$repeats; r++)); do
    for ((i=0; i<$ntis; i++)); do
    idx=`expr $r \* $ntis + $i`
    alltis[$idx]=${alltis[$i]}
    done
done
echo ${alltis[*]}
Log `echo ${alltis[*]}`

### Segmentation of structural image - if we have a structural image we ALWAYS ensure we have a segmentation
if [ ! -z $fslanat ]; then
    Log "Using FSL_ANAT outputs at: $fslanat"
    # we are being supplied with an fslanat directory
    fasthasrun=1 #this means that we have PVE for calibration & PVC purposes
    
    # copy over the things we need and place them using the names used elsewhere
    imcp $fslanat/T1_fast_pve_0 $tempdir/pvcsf_struct #indicate that it is in structural space!
    imcp $fslanat/T1_fast_pve_1 $tempdir/pvgm_struct
    imcp $fslanat/T1_fast_pve_2 $tempdir/pvwm_struct
    chmod u+w $tempdir/pvcsf_struct.* $tempdir/pvgm_struct.* $tempdir/pvwm_struct.*

    if [ ! -z $fslanat/T1_fast_bias ]; then # test to check that there is a bias field here
       Log "Bias field extracted from $fslanat sucessfully"
       imcp $fslanat/T1_fast_bias $tempdir/biasfield_struct
       chmod u+w $tempdir/biasfield_struct.*
    else
    Log "No Bias field found in $fslanat"
    fi
    
elif [ ! -z $struc ]; then
    # do we have the results from FAST already? If not run it
    if [ -z $fastsrc ]; then
    Log "Segmenting the structural image"
    fast -B -b -o $tempdir/seg -p $tempdir/struc_bet
    fastsrc=$tempdir/seg
    else
    # FAST has been run externally
    Log "Using FAST outputs at: $fastsrc"
    fi
    
    # we are now sure we have FAST outputs
    fasthasrun=1

     # copy over the things we need and place them using the names used elsewhere
    imcp ${fastsrc}_pve_0 $tempdir/pvcsf_struct #indicate that it is in structural space!
    imcp ${fastsrc}_pve_1 $tempdir/pvgm_struct
    imcp ${fastsrc}_pve_2 $tempdir/pvwm_struct
    chmod u+w $tempdir/pvcsf_struct.* $tempdir/pvgm_struct.* $tempdir/pvwm_struct.*

    if [ ! -z ${fastsrc}_bias ]; then # test to see if there is a bias field in the FAST output
    Log "Bias field extracted from ${fastsrc} sucessfully"
    imcp ${fastsrc}_bias $tempdir/biasfield_struct
    chmod u+w $tempdir/biasfield_struct.*
    else
    Log "No Bias field found with ${fastsrc}"
    fi
fi
### End of: Segmentation

### Registration (1/2)
# Make sure we have some form of transformation between the ASL data and the structural (if that has been supplied)
# only 'initial' step in asl_reg is used here
register=0
if [ ! -z $strucflag ]; then # if structural image has not been suppled then skip the registration
    register=1
    if [ ! -z $asl2struc ]; then # we have been supplied with a transformation matrix - we do not need registration, but we do want to transform the results
    register=0
    Log "Using existing asl to structural transform: $asl2struc"
    cp $asl2struc $tempdir/asl2struct.mat
    convert_xfm -omat $tempdir/struct2asl.mat -inverse $tempdir/asl2struct.mat
    fi
fi

if [ -z $regfrom_method ]; then
    if [ $iaf = "diff" ] && [ ! -z $calib ]; then
        Log "Defaulting to brain extracted calibration image as initial registration reference"
        regfrom_method="calib"
    else
        Log "Defaulting brain extracted mean ASL image as initial registration reference"
        regfrom_method="mean"
    fi
fi

if [ -z $regfrom ]; then
    # No regfrom iamge supplied. If --regfrom-pwi was specified we will use the PWI, if
    # --regfrom-mean we will use the mean ASL image (subtracted or unsubtracted), if --regfrom-calib
    # we will use the calibration image.
    #
    # The default if nothing is user-specified is to use the calibration image for differenced
    # data (if available) or the mean ASL image if unsubtracted or there is no calibration imag
    if [ $regfrom_method == "pwi" ]; then
        Log "Using brain extracted PWI as initial registration reference"
        bet $tempdir/pwi $tempdir/pwi_brain
        regfrom=$tempdir/pwi_brain
    elif [ $regfrom_method == "calib" ]; then
        Log "Using brain extracted calibration image as initial registration reference"
        bet $tempdir/calib $tempdir/calib_brain
        regfrom=$tempdir/calib_brain
    elif [ $regfrom_method == "mean" ]; then
        Log "Using brain extracted mean ASL image as initial registration reference"
        bet $tempdir/meanasl $tempdir/meanasl_brain -f $regfrom_bet_thresh
        regfrom=$tempdir/meanasl_brain
    else
        Log "ERROR: Unrecognized regfrom_method: $regfrom_method"
        exit 1
    fi
else
    Log "Registration reference image specified by user: $regfrom"
fi

# Create the wmseg here in case it's needed for the initial registration
if [ ! -z $fasthasrun ]; then
    # create a tissseg (wmseg) image for BBR in asl_reg
    Log "Creating WMseg for BBR registration"
    fslmaths $tempdir/pvwm_struct -thr 0.5 -bin ${tempdir}/tissseg
fi

if [ $register -eq 1 ]; then
    # registration here using asl_reg (inital only)
    if [ ! -z $regfrom ]; then
        Log "Performing registration"
        Log "Using $regfrom as base for inital registration"

        if [ ! -z $lowstrucflag ]; then
            extraoptions=$extraoptions" -r $tempdir/lowstruc"
        fi

        if [ ! -z $reg_init_bbr ]; then
            Log "Using BBR for initial registration"
            fslmaths $regfrom -thr 0 -bin $tempdir/mask_reg_init
            extraoptions=$extraoptions" -m $tempdir/mask_reg_init --tissseg $tempdir/tissseg"
            regto=low2high_final
        else
            Log "Initial registration is FLIRT only"
            extraoptions=$extraoptions" --mainonly"
            regto=lowtohigh
        fi

        Registration $regfrom "$extraoptions"

        # Save the registered image
        imcp $tempdir/tmp_asl_reg/$regto $tempdir/asl2struc_initial
    fi
    convert_xfm -omat $tempdir/struct2asl.mat -inverse $tempdir/asl2struct.mat

    # Do a registration of the calibration image as this may be easier to visualise
    if [ ! -z $calib ]; then
        Log "Generating calibration image in structural space"
        flirt -in $tempdir/calib -ref $tempdir/struc -applyxfm -init $tempdir/asl2struct.mat -out $tempdir/calib2struc_initial $APPLYXFM_OPTS
    fi
fi

# Some useful preproc to do with FAST outputs once we have the initial registration
if [ ! -z $fasthasrun ]; then
    if [ ! -z $tempdir/biasfield_struct ]; then
        # transform the bias field and invert to use for sensitivity correction in calibration
        applywarp --ref=$tempdir/nativeref --in=$tempdir/biasfield_struct --out=$tempdir/biasfield --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
    fi

    if [ ! -z $senscorr ]; then
        if [ ! -z $tempdir/biasfield ]; then #make sure we have the biasfield (from above) before attempting this
            Log "Creating sensitivity map from biasfield"
            fslmaths $tempdir/biasfield -recip $outdir/native_space/sensitivity
            imcp $outdir/native_space/sensitivity $tempdir/sensitivity
            chmod u+w $tempdir/sensitivity.*
        fi
    fi
fi
### End of: Registration (1/2)

### Generate mask for ASL data
# sort out the mask for processing the data
if [ -z $mask ]; then
Log "Creating mask"
# preferred option is to use brain extracted structural
if [ ! -z $strucflag ]; then
    fslmaths $tempdir/struc_bet -bin $tempdir/struc_bet_mask
    flirt -in $tempdir/struc_bet_mask -ref $regfrom -applyxfm -init $tempdir/struct2asl.mat -out $tempdir/mask  $APPLYXFM_OPTS
    fslmaths $tempdir/mask -thr 0.25 -bin -fillh $tempdir/mask
    fslcpgeom $regfrom $tempdir/mask
# otherwise use the regfrom image (should already be BETed) - note that regfrom may have been set in the Registration (1/2) section.
elif [ ! -z $regfrom ]; then
    fslmaths $regfrom -bin $tempdir/mask
    Log "Mask generated from regfrom image: $regfrom"
    
# We are unlikey to use these options as regfrom has probably been set above - leave for future reference
# next option is to use betted version of mean M0 calib image as mask
elif [ ! -z $calib ]; then
    bet $tempdir/calib $tempdir/calib_bet
    fslmaths $tempdir/calib_bet -bin $tempdir/mask
    Log "Mask generated from calibration image (post BET)"
 # use the low resolution strucutral image to create mask (ahould already be BETed)
elif [ ! -z $lowstrucflag ]; then
 # resample
    flirt -in $tempdir/lowstruc -applyxfm -init $FSLDIR/etc/flirtsch/ident.mat -out $tempdir/mask -paddingsize 0.0 -interp trilinear -ref $tempdir/nativeref
 # make binary
    fslmaths $tempdir/mask -bin $tempdir/mask
    Log "Mask generated from low res. structural"
# otherwise just use mean time series
else
    bet $tempdir/meanasl $tempdir/meanasl -f $regfrom_bet_thresh
    fslmaths $tempdir/meanasl -bin $tempdir/mask
    Log "Mask generated from mean time series"
fi
else
# mask has been supplied
    fslmaths $mask -bin $tempdir/mask # just to be sure binarise the mask here
    Log "Using mask: $mask"
fi

# Generate dilated mask for BASIL analysis to ensure that when we do the final registration
# and re-generate the mask, (hopefully) all the voxels in the final mask have been through
# BASIL. We don't need to do this when the mask was user-specified as then it is unaffected
# by registration
if [ -z $mask ]; then
    Log "Generating dilated mask for initial BASIL runs"
    fslmaths $tempdir/mask -kernel boxv 5 -dilF $tempdir/mask_dil
else
    imcp $tempdir/mask $tempdir/mask_dil
fi

### End of: Generate mask

### Distortion Correction
# Do TOPUP if applicable
if [ ! -z $cblip ]; then
    if [ -z $calib ]; then
    Log "WARNING: Cannot do TOPUP on blip-reversed calibration image ($cblip) without correpsonding calibration image"
    elif [ -z $echospacing ] || [ -z $pedir ]; then
    Log "WARNING: Cannot do TOPUP on blip-reversed calibration image without echospacing (dwell time) and pahse encode direction"
    else
    Log "Distortion correction: running topup"
    
    # Create topup params. Note that TOPUP requires to total readout time which is the echo spacing multiplied
    # by the output dimension size - 1. We are interpreting the xyz axis as being the axes of the underlying
    # voxel array - normally these will correspond to RAS space but this is not guaranteed. However this seems
    # to be what other FSL tools use which take an encoding direction.
    case $pedir in
        x)
        dimsize=`fslinfo $cblip | grep "^dim1" | sed 's:dim1[ ]*::'`
        totaltime=`echo "$echospacing * ($dimsize-1)" | bc`
        echo "1 0 0 $totaltime" > $tempdir/topup_params.txt
        echo "-1 0 0 $totaltime" >> $tempdir/topup_params.txt
        ;;
        -x)
        dimsize=`fslinfo $cblip | grep "^dim1" | sed 's:dim1[ ]*::'`
        totaltime=`echo "$echospacing * ($dimsize-1)" | bc`
        echo "-1 0 0 $totaltime" > $tempdir/topup_params.txt
        echo "1 0 0 $totaltime" >> $tempdir/topup_params.txt
        ;;
        y)
        dimsize=`fslinfo $cblip | grep "^dim2" | sed 's:dim2[ ]*::'`
        totaltime=`echo "$echospacing * ($dimsize-1)" | bc`
        echo "0 1 0 $totaltime" > $tempdir/topup_params.txt
        echo "0 -1 0 $totaltime" >> $tempdir/topup_params.txt
        ;;
        -y)
        dimsize=`fslinfo $cblip | grep "^dim2" | sed 's:dim2[ ]*::'`
        totaltime=`echo "$echospacing * ($dimsize-1)" | bc`
        echo "0 -1 0 $totaltime" > $tempdir/topup_params.txt
        echo "0 1 0 $totaltime" >> $tempdir/topup_params.txt
        ;;
        z)
        dimsize=`fslinfo $cblip | grep "^dim3" | sed 's:dim3[ ]*::'`
        totaltime=`echo "$echospacing * ($dimsize-1)" | bc`
        echo "0 0 1 $totaltime" > $tempdir/topup_params.txt
        echo "0 0 -1 $totaltime" >> $tempdir/topup_params.txt
        ;;
        -z)
        dimsize=`fslinfo $cblip | grep "^dim3" | sed 's:dim3[ ]*::'`
        totaltime=`echo "$echospacing * ($dimsize-1)" | bc`
        echo "0 0 -1 $totaltime" > $tempdir/topup_params.txt
        echo "0 0 1 $totaltime" >> $tempdir/topup_params.txt
        ;;
    esac
    
    Log "For TOPUP using total read-out time of ($dimsize - 1) x $echospacing s = $totaltime s"

    # do topup
    fslmerge -t $tempdir/calib_blipped $tempdir/calib $tempdir/cblip 
    topup --imain=$tempdir/calib_blipped --datain=$tempdir/topup_params.txt --config=b02b0.cnf --out=$tempdir/topupresult --fout=$tempdir/topupresult_fmap
    topupresult=$tempdir/topupresult
    fi
fi

#Fieldmaps
if [ ! -z $topupresult ]; then
    # Moved the application of TOPUP correction below since otherwise it could be overwritten by user GDC warp
    :
elif [ ! -z $fmap ]; then
    # a fieldmap has been provided that needs registration - copy images over
    imcp $fmap $tempdir/fmap
    imcp $fmapmag $tempdir/fmapmag
    imcp $fmapmagbrain $tempdir/fmapmagbrain
    chmod u+w $tempdir/fmap.* $tempdir/fmapmag.* $tempdir/fmapmagbrain.*
fi

if [ -e $tempdir/fmap.* ]; then
    Log "Distortion Correction using asl_reg"
    
    # Do registration to T1 to get distortion correction warp
    # use whatever registration matrix we already have to initialise here 
    distbase=$tempdir/pwi # use the perfusion-weighted image (mean over all TIs) as the best basis we have for registration at this point
    if [ -z $finalreg ]; then
    distbase=$regfrom # use whatever image we have been using for (inital) registration
    fi

    Registration $distbase "-m $tempdir/mask --tissseg $tempdir/tissseg --imat $tempdir/asl2struct.mat --finalonly" distcorr
    
    # Generate the correction warp in ASL space
    convertwarp -r $tempdir/nativeref -o $tempdir/asldist_warp -w $tempdir/distcorr/asl2struct_warp.nii.gz --postmat=$tempdir/distcorr/struct2asl.mat --rel -j $tempdir/distcorr/jacobian_parts --constrainj

    # If we ALSO have a user supplied GDC warp then we need to merge that with the one we have found here
    # use convert warp to combine the warps. Note that we cannot do this in a single step becase we need
    # to use --postmat above to transform into ASL space  but this needs to happen before we apply
    # the user warp 
    if [ ! -z $gdcwarp ]; then
        Log "Adding user-supplied GDC warp to distortion correction warp"
        convertwarp -r $tempdir/nativeref -o $tempdir/asldist_warp -w $tempdir/asldist_warp --warp2=$tempdir/asldist_userwarp --rel -j $tempdir/distcorr/jacobian_parts --constrainj
    fi

    # Calculation of the jacobian for the warp - method suggested in:
    # https://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=FSL;d3fee1e5.0908
    fnirtfileutils -i $tempdir/asldist_warp -f spline -o $tempdir/asldist_warp_coef
    fnirtfileutils -i $tempdir/asldist_warp_coef -j $tempdir/distcorr/jacobian
    fslcpgeom $tempdir/asldist_warp $tempdir/distcorr/jacobian -d

elif [ ! -z $gdcwarp ]; then
    Log "Re-applying user-specified GDC warp"
    imcp $tempdir/asldist_userwarp $tempdir/asldist_warp
fi

if [ -e $tempdir/asldist_warp.* ]; then  
    # Re-apply combined warp/moco correction to the data
    imcp $tempdir/asldata $tempdir/asldata_origwarp
    ApplyWarps $tempdir/asldist_warp $tempdir/distcorr/jacobian $tempdir/asldata.cat $tempdir/calib2asl.mat 
    repeatsubtract=1;
fi

if [ ! -z $topupresult ]; then
    # Apply TOPUP correction here because otherwise it could be overwritten by user GDC warp. Note that we would like to combine these
    # in one transformation to avoid multiple interpolations, however this does not seem to be supported by applytopup

#    if [ -e $tempdir/struc.* ]; then
#    # currently DISABLED and applytopup is used with topup results
#    #we will do the distorition correction using epi_reg so that we can merge with motion correction matrices and also get the jacobian
#    # the fieldmap provided is from topup and will be in ASL space
#    # convert ot the correct units of rad/s from Hz
#    fslmaths ${topupresult}_fmap -mul 3.1459 -mul 2 $tempdir/fmap
#    # use the existing registration to get the fieldmap into T1 space
#    flirt -in $tempdir/fmap -out $tempdir/fmap -ref $tempdir/struc -applyxfm -init $tempdir/asl2struct.mat
#    # asl_reg/epi_reg will expect a fieldmap magnitude image (although doesn't really need it in this case) - just use the structural
#    imcp $tempdir/struc $tempdir/fmapmag
#    imcp $tempdir/struc_bet $tempdir/fmapmagbrain
#    nofmapreg=1
#    else
    Log "Distortion Correction using TOPUP"
    # we will use apply topup - this does not do the jacboian magntiude correction - therefore strictly only okay if using voxelwise calibration
    applytopup --imain=$tempdir/calib,$tempdir/cblip --inindex=1,2 --datain=$tempdir/topup_params.txt --topup=${topupresult} --out=$tempdir/calib --method=jac
    repeatsubtract=1;
    applytopup --imain=$tempdir/asldata --datain=$tempdir/topup_params.txt --inindex=1 --topup=${topupresult} --out=$tempdir/asldata --method=jac #ND using asldata as this has been motion corrected by this point (if requested)
    if [ ! $cmethod="voxel" ]; then
        Log "WARNING: Using apply_topup this does not correct for magntiude using the jocbian in distortion correction - this is not optimal when not using voxelwise calibration, to avoid this supply structural image(s)"
    fi
    if [ ! -z $cref ]; then
        applytopup --imain=$tempdir/cref --datain=$tempdir/topup_params.txt --inindex=1 --topup=${topupresult} --out=$tempdir/cref --method=jac
    fi
fi

# Repeat the label-control subtraction on the corrected data
if [ ! -z $repeatsubtract ]; then
    if [ $iaf = 'diff' ]; then
    # make sure the block format is correct for BASIL
    asl_file --data=$tempdir/asldata --ntis=$ntis --ibf=$ibf $afrptstr --iaf=$iaf --obf=tis --out=$tempdir/diffdata --mean=$tempdir/diffdata_mean
    else
    # create label-control difference data using asl_file - this gets it into the right block form for BASIL (blocks of TIs)
    asl_file --data=$tempdir/asldata --ntis=$ntis --ibf=$ibf $afrptstr --iaf=$iaf --obf=tis --diff --out=$tempdir/diffdata --mean=$tempdir/diffdata_mean
    fi
fi
### End of: Distortion Correction

# Mask the calibration images, saving the wholehead images (although currently unused)
if [ ! -z $calib ]; then
    imcp $tempdir/calib $tempdir/calib_wholehead
    fslmaths $tempdir/calib -mas $tempdir/mask_dil $tempdir/calib
fi
if [ ! -z $cref ]; then
    imcp $tempdir/cref $tempdir/cref_wholehead
    fslmaths $tempdir/cref -mas $tempdir/mask_dil $tempdir/cref
fi

# Copy or recalculate the sensitivity map (overwriting any created
# by FAST) if we have a cref image or a sensitivity image has been supplied directly
if [ ! -z $isen ]; then
    # User-supplied sensitivity image
    Log "WARNING: User supplied sens image"
    fslmaths $isen -mas $tempdir/mask_dil $outdir/native_space/sensitivity
    imcp $outdir/native_space/sensitivity $tempdir/sensitivity
elif [ ! -z $cref ]; then
    # Calculate sensitivty image using user-supplied cref image
    fslmaths $tempdir/calib -div $tempdir/cref -mas $tempdir/mask_dil $outdir/native_space/sensitivity
    imcp $outdir/native_space/sensitivity $tempdir/sensitivity
fi

# Senstivity correction cannot be done if this image hasn't been generated by this point
if [ -z $outdir/native_space/sensitivity ]; then
    senscorr=""
    Log "WARNING: sensitivity correction has been requested, but suitable map is not available, skipping that step"
fi

# Defaults for (some) parameters
# deal with T1
if [ -z $t1set ]; then
    t1set=1.3;
fi
Log "T1: $t1set"

if [ -z $t1bset ]; then
# follws the ASL white paper recommendation
    t1bset=1.65;
fi
Log "T1b: $t1bset"

#pre-processing for epochwise analysis
# separate data into epochs here
if [ ! -z $epoch ]; then
if [ -z $eol ]; then
    eol=0
fi
    asl_file --data=$tempdir/diffdata --ntis=$ntis --ibf=tis --iaf=diff --epoch=$tempdir/epoch --elen=$elen --eol=$eol --eunit=tis
    eadv=`expr $elen - $eol`
fi


# write options file for BASIL - these are the core options that are appropraite whether we are doing a single or epochwise analysis
Log "Setting up BASIL"

# T1 values
echo "--t1=$t1set --t1b=$t1bset" >> $tempdir/basil_options.txt

if [ ! -z $t1im ]; then
    basil_options=$basil_options"--t1im=$t1im "
    #echo "--t1im=$t1im" >> $tempdir/basil_options.txt
    Log "Using supplied T1 (tissue) image in BASIL: $t1im"
fi

if [ ! -z $batim ]; then
    basil_options=$basil_options"--batim=$batim "
    Log "Using supplied BAT (seconds) image in BASIL: $batim"
fi

if [ ! -z $artonly ]; then
    basil_options=$basil_options"--artonly "
    Log "BASIL fitting restricted to arterial component only"
fi

# data acquired using CASL?
if [ ! -z $casl ]; then
    echo "--casl" >> $tempdir/basil_options.txt;
    Log "cASL model"
else
    Log "pASL model"
fi

#echo "--tau=$boluslen" >> $tempdir/basil_options.txt
echo $tauslist >> $tempdir/basil_options.txt


# slice timing correction?
if [ ! -z $slicedt ]; then
    echo "--slicedt=$slicedt" >> $tempdir/basil_options.txt
    Log "Slice timing correction with delta: $slicedt"
fi

# Flip anlge for look-locker readout
if [ ! -z $fa ]; then
    echo "--FA=$fa" >> $tempdir/basil_options.txt
    Log "Flip angle (look-locker readout): $fa"
fi

# Multi-band setup (if not set then this is ignored)
if [ ! -z $sliceband ]; then
    echo "--sliceband=$sliceband" >> $tempdir/basil_options.txt
    Log "Multi-band setup with number of slices per band: $slicedband"
fi

# Infer arterial component?
if [ -z $artoff ]; then
    basil_options=$basil_options"--inferart "
    Log "Infer arterial component"
fi
# fix the bolus duration?
if [ -z $fixbolus ]; then
    basil_options=$basil_options"--infertau "
    Log "Varaiable bolus duration"
else
    Log "Fixed bolus duration"
fi

#deal with BAT
if [ -z $bat ]; then
    if [ -z $casl ]; then
    bat=0.7 #PASL default
    else
    bat=1.3 #CASL default
    fi
fi
echo "--bat=$bat" >> $tempdir/basil_options.txt  

if [ -z $fixbat ]; then
    Log "Variable arterial arrival time"
    Log "Setting prior/initial (tissue/gray matter) bolus arrival time to $bat"

    # Tissue BAT SD
    #defaults
    if [ -z $singleti ]; then
        if [ -z $batsd ]; then
            # multi TI/PLD data, default to a more liberal prior for tissue ATT since we should be able to determine from data
            # NB this leave the arterial BAT alone.
            batsd=1;
        fi
    fi

    # if required add batsd option for basil
    if [ ! -z $batsd ]; then
        Log "Setting std dev of the (tissue) BAT prior to $batsd"
        echo "--batsd=$batsd" >> $tempdir/basil_options.txt
    fi
else
    basil_options=$basil_options"--fixbat " 
    Log "Fixed arterial arrival time"
    Log "Setting arterial arrival time to $bat"
fi

# Noise specification
if [ -z $snr ]; then
    snr=10; #default SNR
fi
if [ $tpoints -le 5 ]; then
    # only small number of time points in data, will use informative noise prior
   noiseprior=1
   Log "Single volume: informative noise prior will be used"
fi

if [ ! -z $noiseprior ]; then
    # use an informative nosie prior
     if [ -z $noisesd ]; then
     Log "Using SNR of $snr to set noise std dev"
    # estimate signal magntiude
    fslmaths $tempdir/diffdata_mean -Tmax $tempdir/datamax
    brain_mag=`fslstats $tempdir/datamax -k $tempdir/mask -M`
    # this will correspond to whole brain CBF (roughly) - about 0.5 of GM
    noisesd=`echo "scale=2;sqrt( $brain_mag * 2 / $snr )" | bc`
    fi

    Log "Using a prior noise sd of: $noisesd"
    echo "--prior-noise-stddev=$noisesd" >> $tempdir/basil_options.txt
fi


# Exteneded options for BASIL
if [ ! -z $spatial ]; then
# if we are using spatial smoothing on CBF then we will also do the analysis in a single step
    Log "Instructing BASIL to use automated spatial smoothing"
    basil_options=$basil_options"--spatial "

fi

if [ ! -z $infert1 ]; then
    Log "Instructing BASIL to infer variable T1 values"
    basil_options=$basil_options"--infert1 "
fi


if [ ! -z $exch ]; then
    # use a specific exchange model in BASIL
    Log "Using exchange model: $exch"
    basil_options=$basil_options"--exch=$exch "
fi

if [ ! -z $disp ]; then
    # use a specific dispersion model in BASIL
    Log "Using dispersion model: $disp"
    basil_options=$basil_options"--disp=$disp "
fi

if [ ! -z $devel ]; then
    basil_options=$basil_options" --devel "
fi

# Save model fit by default
echo "--save-model-fit" >> $tempdir/basil_options.txt

if [ ! -z $model_options ]; then
    Log "Appending additional BASIL options to $tempdir/basil_options.txt"
    cat $model_options >> $tempdir/basil_options.txt
fi

Log "BASIL options ($tempdir/basil_options.txt):"
Log "----"
Log "`cat $tempdir/basil_options.txt`"
Log "----"
# -- end of main basil options setting

cp $tempdir/basil_options.txt $tempdir/basil_options_core.txt # keep a copy of the core options accumulated thus far (we might need these again for the epoch analysis)


##### Analyse data using BASIL
### First analysis on whole data, normal perfusion image
Log "Calling BASIL on data - conventional perusion image"
initbasil="" # Can be used to pass an intital MVN into the main run of BASIL from an intial run (below)
if [ $repeats -gt 1 ] || [ $repeats -eq 0 ]; then
    # do an initial analysis using the data averaged at each TI
    # NB repeats=0 is a special case of variable number of repeats at each TI
    Log "Initial run of BASIL on data where we have avareged all repeats at each TI"
    datafile=$tempdir/diffdata_mean
    mkdir $tempdir/init
    cat $tempdir/basil_options_core.txt > $tempdir/init/basil_options.txt
    echo "--repeats=1" >> $tempdir/init/basil_options.txt
    echo "$tislist" >> $tempdir/init/basil_options.txt
    Dobasil $datafile $tempdir/init
    initbasil=$tempdir/init/finalMVN
fi

# main analysis using full data
datafile=$tempdir/diffdata
if [ $repeats -gt 0 ]; then
    echo "--repeats=$repeats" >> $tempdir/basil_options.txt
else
    # variable number of repeats at each TI - tell basil
    echo "$rptslist" >> $tempdir/basil_options.txt
fi
echo "$tislist" >> $tempdir/basil_options.txt

Log "Main run of BASIL on ASL data"
Dobasil $datafile $tempdir $initbasil
### End of: First analysis on whole data

### Registration (2/2)
# Revisit the registration now that we have a pefusion image (with good GM/WM contrast) using BBR
# use existing registration for initial aligment
if [ $register -eq 1 ]; then
    if [ $finalreg -eq 1 ]; then
    Log "Performing final registration"
    cp $tempdir/asl2struct.mat $outdir/native_space/asl2struct_init.mat #preserve the intial registration for future reference
    Registration $tempdir/ftiss "-m $tempdir/mask_dil --tissseg $tempdir/tissseg --imat $tempdir/asl2struct.mat --finalonly"
    fi

    # Do a registration of the calibration image too as this may be easier to visualise
    if [ ! -z $calib ]; then
        Log "Generating calibration image in structural space"
        flirt -in $tempdir/calib -ref $tempdir/struc -applyxfm -init $tempdir/asl2struct.mat -out $tempdir/calib_struc $APPLYXFM_OPTS
    fi
fi

# Re-generate the mask from the structural image now that we have a better registration - as this might matter for PV correction
# NB we dont use the PVE here since we dont (necessarily) want to exclude the ventricles from the mask as this has implications for the spatial priors
# NB Also, only regenerate the mask if there was no user-specified mask given
if [ -f $tempdir/struc_bet.* -a -z "$mask" ]; then
    Log "Generating analysis mask from final structural registration"
    imcp $tempdir/mask $tempdir/mask_orig
    imcp $tempdir/mask_dil $tempdir/mask_dil_orig
    fslmaths $tempdir/struc_bet -bin $tempdir/struc_bet_mask
    flirt -in $tempdir/struc_bet_mask -ref $tempdir/nativeref -applyxfm -init $tempdir/struct2asl.mat -out $tempdir/mask $APPLYXFM_OPTS
    fslmaths $tempdir/mask -thr 0.25 -bin -fillh $tempdir/mask
    fslcpgeom $regfrom $tempdir/mask
    imcp $tempdir/mask $tempdir/mask_dil

    fslmaths $tempdir/mask -sub $tempdir/mask_dil_orig -thr 0 $tempdir/missed_voxels
    num_missed=`fslstats $tempdir/missed_voxels -V |awk '{print $1}'`
    if [ $num_missed -ne 0 ]; then
        Log "WARNING: Final analysis mask contains $num_missed voxels that were not in original model fitting mask"
    fi
fi


### End of: Registration (2/2)

### Partial Volume Estimates
# Note we do this here since we have the final registration now which we need to transform PV estimates into ASL space
if [ ! -z $fasthasrun ] && [ -z $pvgm ]; then
    # PVE in ASL space from strcutural segmentation results
    # invert the transformation matrix
    convert_xfm -omat $tempdir/struct2asl.mat -inverse $tempdir/asl2struct.mat
    
    # Gray matter - assume this will be PVE 1
    applywarp --ref=$tempdir/nativeref --in=$tempdir/pvgm_struct --out=$tempdir/pvgm_inasl --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
    # white matter  - assume this will be PVE 2
    applywarp --ref=$tempdir/nativeref --in=$tempdir/pvwm_struct --out=$tempdir/pvwm_inasl --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
    # threshold (upper and lower) the PVE to avoid artefacts of spline interpolation and also ignore very low PVE that could cause numerical issues.
    fslmaths $tempdir/pvgm_inasl -thr 0.1 -min 1 $tempdir/pvgm_inasl
    fslmaths $tempdir/pvwm_inasl -thr 0.1 -min 1 $tempdir/pvwm_inasl
    pvexist=1
fi

if [ ! -z $pvgm ]; then
    #using supplied PV images
    Log "Loading supplied PV images"
    if [ -z $pvwm ]; then
        Log "ERROR: no WM PV image has been supplied"
    fi
    Log "PV GM is: $pvgm"
    fslmaths $pvgm -thr 0.1 -min 1 $tempdir/pvgm_inasl
    Log "PV WM is: $pvwm"
    fslmaths $pvwm -thr 0.1 -min 1 $tempdir/pvwm_inasl
    pvexist=1
fi

# Thresholds for voxels to be considered 'pure' GM / WM
if [ -z $gm_thresh ]; then
    gm_thresh=0.8;
fi
if [ -z $wm_thresh ]; then
    wm_thresh=0.9;
fi

if [ ! -z $pvexist ]; then
    # make some masks 
    # these are currently used for masking after model fitting
    fslmaths $tempdir/pvgm_inasl -thr 0.1 -mas $tempdir/mask -bin $tempdir/gmmask
    fslmaths $tempdir/pvwm_inasl -thr 0.1 -mas $tempdir/mask -bin $tempdir/wmmask
    # these are for calculating mean perfusion within tissue types
    fslmaths $tempdir/pvgm_inasl -thr $gm_thresh -mas $tempdir/mask -bin $tempdir/gmmask_pure
    fslmaths $tempdir/pvwm_inasl -thr $wm_thresh -mas $tempdir/mask -bin $tempdir/wmmask_pure

    # If we have a structural -> std space transformation, generate cortical-only masks for GM
    # and WM using FSL atlases. This involves using the 'Cerebral Cortex' and 'Cerebral White Matter'
    # maps from the Harvard-Oxford subcortical mask and summing them to get a cortical mask that
    # excludes subcortical regions including the cerebellum
    if [ ! -z $warp ]; then
        struc2std_warp=$warp
        Log "Using user-supplied structural->standard space warp to generate cortical GM/WM masks"
    elif [ -f "$fslanat/T1_to_MNI_nonlin_coeff."* ]; then
        Log "Using FSL_ANAT structural->standard space warp to generate cortical GM/WM masks"
        struc2std_warp=$fslanat/T1_to_MNI_nonlin_coeff
    elif [ -e "$fslanat/T1_to_MNI_lin.mat" ]; then
        Log "Using FSL_ANAT structural->standard space linear transformation to generate cortical GM/WM masks"
        struc2std_trans=$fslanat/T1_to_MNI_lin.mat
    fi

    if [ ! -z $struc2std_warp ] || [ ! -z $struc2std_trans ]; then
        fslroi $FSLDIR/data/atlases/HarvardOxford/HarvardOxford-sub-prob-2mm.nii.gz $tempdir/lcort_std 0 2
        fslroi $FSLDIR/data/atlases/HarvardOxford/HarvardOxford-sub-prob-2mm.nii.gz $tempdir/rcort_std 11 2
        fslmaths $tempdir/lcort_std -add $tempdir/rcort_std -Tmean -mul 2 $tempdir/cort_std
        if [ ! -z $struc2std_warp ]; then
            invwarp --ref=$tempdir/struc --warp="$struc2std_warp" --out=$tempdir/std2struc_warp
            applywarp --in=$tempdir/cort_std --out=$tempdir/cort_asl --ref=$tempdir/nativeref --warp=$tempdir/std2struc_warp --postmat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
        else
            convert_xfm -omat "$struc2std_trans" -inverse $tempdir/std2struc_trans
            applywarp --in=$tempdir/cort_std --out=$tempdir/cort_asl --ref=$tempdir/nativeref --premat=$tempdir/std2struc_trans --postmat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
        fi

        fslmaths $tempdir/cort_asl -thr 50 -bin $tempdir/cort_mask_asl
        fslmaths $tempdir/gmmask_pure -mas $tempdir/cort_mask_asl $tempdir/gmmask_pure_cort
        fslmaths $tempdir/wmmask_pure -mas $tempdir/cort_mask_asl $tempdir/wmmask_pure_cereb
    else
        Log "WARNING: No structural->standard space transformation - cannot generate cortical-only GM/WM masks"
    fi
fi
### End of: Partial Volume Estimates

### Calibration
# Do calibration here becuase we do not need it before this point & if we are generating a CSF mask we have a better registration at this point
if [ -z $t1tset ]; then
    t1tset=1.3;
fi
Log "T1t (for calibration): $t1tset"

# TR (for calibration image)
if [ -z $tr ]; then
    tr=3.2
fi


# calibration image gain
if [ -z $cgain ]; then
    cgain=1;
fi

# Calibration if reqd
if [ -z $alpha ]; then
        # based on the ASL white paper
    if [ -z $casl ]; then
    alpha=0.98;
    else
    alpha=0.85;
    fi
fi

if [ -z $cmethod ]; then
# default calibration method is 'voxelwise' unless we have CSF PV estimates or CSF mask has been supplied
    if [ ! -z $fasthasrun ] || [ ! -z $csf ]; then
    cmethod=single
    else
    cmethod=voxel
    fi
fi

if [ ! -z $calib ]; then

    if [ $cmethod = 'single' ] || [ $cmethod = 'both' ]; then
        # Single M0 value for calibration
        Log "Doing calibration using a single M0 value with a CSF reference"
        if [ -z $csf ] && [ -z $fasthasrun ]; then
            Log "ERROR: Provide either a structural image or CSF mask for calibration when using --cmethod=single"
            exit 1
        fi
        # calcualte M0a from CSF
        Calibration
        Mo=`cat $outdir/calib/M0.txt`
    fi
    if [ $cmethod = 'voxel' ] || [ $cmethod = 'both' ]; then
        # Voxelwise M0 values for calibration
        Log "Doing calibration using voxelwise M0 map"
        mkdir $outdir/calib

        # Copy over the calibration image and apply the cgain setting - this increases the magntiude of M0 to match that of the ASL data (acquired with a higher gain - cgain>=1 normally)
        fslmaths $tempdir/calib -Tmean -mul $cgain $outdir/calib/M0
        Momap=$outdir/calib/M0
        if [ 1 -eq `echo "if($tr < 5){a=1};a" | bc`  ]; then
            # correct the M0 image for short TR using the equation from the white paper
            Log "Correcting the calibration (M0) image for short TR (using T1 of tissue $t1tset)"
            ccorr=`echo "1 / (1 - e(- $tr / $t1tset) )" | bc -l`
            fslmaths $Momap -mul $ccorr $Momap
        fi

        # Include partition co-effcient in M0 image to convert from M0 tissue to M0 arterial
        fslmaths $Momap -div 0.9 $Momap

        if [ ! -z $edgecorr ]; then
            # Correct for (partial volume) edge effects
            # Median smoothing and erosion
            fslmaths $Momap -fmedian -mas $tempdir/mask -ero $tempdir/calib_ero
            # Extrapolation to match mask
            asl_file --data=$tempdir/calib_ero --ntis=1 --mask=$tempdir/mask --extrapolate --neighbour=5 --out=$Momap
        fi
    elif [ $cmethod != 'single' ]; then
        Log "Error unrecognised calibration method: $cmethod, (use single, voxel or both)"
    fi 
elif [ ! -z $M0 ]; then
    # An M0 value has been supplied, use this
    cmethod=single # we are in 'single' mode as a single value has been supplied
    Mo=$M0
    Log "M0: $Mo"
    Log "Using supplied M0 value: $Mo"
fi

### End of: Calibration

### Output main BASIL results
if [ $cmethod = 'both' ]; then
    # Output with voxelwise and refregion calibration in separate subdirs
    Dooutput / voxelwise voxel
    Dooutput / refregion single
else
    # Only one calibration method so can use default output dir
    Dooutput / / $cmethod
fi

# save the mask used to the (native space) output directory. Note this is the final reg mask
imcp $tempdir/mask $outdir/native_space/mask

### End of: Output main BASIL results

### Partial Volume Correction BASIL
if [ ! -z $pvcorr ]; then
    # Structrual -> ASL registration may have changed slightly, so need to re-calculate sensitivity map, since it
    # either depends on the brain mask from the structural image, or the bias field extracted by FASL/FSL_ANAT
    # Save the original sensitivity map for comparison
    imcp $outdir/native_space/sensitivity $outdir/native_space/sensitivity_non_pvcorr
    if [ ! -z $isen ]; then
      # User-supplied sensitivity image
      Log "WARNING: User supplied sens image"
      fslmaths $isen -mas $tempdir/mask $outdir/native_space/sensitivity
    elif [ ! -z $cref ]; then
      # Calculate sensitivty image using user-supplied cref image
      fslmaths $tempdir/calib -div $tempdir/cref -mas $tempdir/mask $outdir/native_space/sensitivity
    elif [ ! -z $fasthasrun ]; then
      if [ ! -z $tempdir/biasfield_struct ]; then
        if [ ! -z $senscorr ]; then
          # Struct->ASL registration may have changed, so re-transform the bias field and invert to use for sensitivity correction
          Log "Re-registering bias field and calculating sensitivity map"
          applywarp --ref=$tempdir/nativeref --in=$tempdir/biasfield_struct --out=$outdir/native_space/biasfield --premat=$tempdir/struct2asl.mat $APPLYWARP_OPTS
          fslmaths $outdir/native_space/biasfield -recip $outdir/native_space/sensitivity
        fi
      fi
    fi
    
    # intructions for BASIL
    basil_options=$basil_options" --pgm $tempdir/pvgm_inasl --pwm $tempdir/pvwm_inasl "
    mkdir $tempdir/pvcorr
    cp $tempdir/basil_options.txt $tempdir/pvcorr/basil_options.txt #Dobasil expects the options file to be in the subdirectory
    # Run BASIL
    basil_pvcorr=$pvcorr
    Dobasil $datafile $tempdir/pvcorr

    imcp $tempdir/pvcorr/finalMVN $tempdir/finalMVN #just in case we are about to do a epochwise analysis

    # Output the PVC results
    if [ $cmethod = 'both' ]; then
        # Output with voxelwise and refregion calibration in separate subdirs
        Dooutput pvcorr voxelwise/pvcorr voxel
        Dooutput pvcorr refregion/pvcorr single
    else
        # Only one calibration method so can use default output dir
        Dooutput pvcorr pvcorr $cmethod
    fi
fi
### End of: Partial Volume Correction

### Epoch BASIL
if [ ! -z $epoch ]; then
    # epochwise analysis
    Log "Epochwise analysis"

    #genereate a list of epochs
    currdir=`pwd`
    cd $tempdir
    epochlist=`imglob epoch*`
    cd $currdir

    ecount=0
    for e in $epochlist; do
    Log "Processing epoch: $e"
    etislist=""
        # deal with the TIs
    for ((ei=0; ei<$elen; ei++)); do
        ethis=`expr $ecount \* $eadv + $ei`
        #echo $ethis
        eidx=`expr $ei + 1`
        #echo $ei
        #echo ${alltis[$ethis]}
        etislist=$etislist" --ti${eidx}=${alltis[$ethis]}"
    done
    Log "TIs for this epoch: "
    Log $etislist

    mkdir $tempdir/$e
    cp $tempdir/basil_options_core.txt $tempdir/$e/basil_options.txt # get the 'core' options and make a new basil_options file jsut for this TI
    echo "--repeats=1" >> $tempdir/$e/basil_options.txt #for epochs we specify all the TIs explicitly
    echo $etislist >>  $tempdir/$e/basil_options.txt #these are the basil options for this epoch

    fast=2 #we now switch BASIL to fast level '2' - this means it will only do analysis in a single step from here on in, but we will use our existing analysis for initialisation.
    
    Dobasil $tempdir/$e $tempdir/$e $tempdir/basil/finalMVN # init with results of first basil run

    # Output epochwise results
    Log "Saving results from epoch: $e"
    if [ $cmethod = 'both' ]; then
        # Output with voxelwise and refregion calibration in separate subdirs
        Dooutput / voxelwise/$e voxel
        Dooutput / refregion/$e single
    else
        # Only one calibration method so can use default output dir
        Dooutput / $e $cmethod
    fi

    ecount=`expr $ecount + 1`
    done
fi
### End of: Epoch BASIL




#OUTPUTS
# Setup option outputs - anything that would be common to all epochs
# note that we now do directory creation right at the start
#if [ ! -z $nativeout ]; then
#fi
if [ ! -z $structout ]; then
    #cp $tempdir/asl2struct.mat $outdir/struct_space/asl2struct.mat
    cp $tempdir/asl2struct.mat $outdir/native_space/asl2struct.mat #also provide the transformation matrix for reference
fi
#if [ ! -z $advout ]; then
#fi

#if [ -z $epoch ]; then
# normal single analysis of data
#Dooutput

##if [ ! -z $epoch ]; then
# epochwise analysis
#    for e in $epochlist; do
#    Log "Saving results from epoch: $e"
#    Dooutput $e
 #   done
#fi




if [ ! -z $pvcorr ]; then
# copy PVE in ASL space to output directory
imcp $tempdir/pvgm_inasl $outdir/native_space/pvgm_inasl
imcp $tempdir/pvwm_inasl $outdir/native_space/pvwm_inasl
fi

if [ ! -z $pvexist ]; then
    # copy PV masks to output directory
    imcp $tempdir/gmmask $outdir/native_space/gm_mask
    imcp $tempdir/wmmask $outdir/native_space/wm_mask
    imcp $tempdir/gmmask_pure $outdir/native_space/gm_roi
    imcp $tempdir/wmmask_pure $outdir/native_space/wm_roi
    imcp $tempdir/gmmask_pure_cort $outdir/native_space/cortical_gm_roi
    imcp $tempdir/wmmask_pure_cereb $outdir/native_space/cerebral_wm_roi
fi

# If using both calibration methods move all files into both subdirs and delete from generic output dir
if [ $cmethod = 'both' ]; then
    cp $outdir/native_space/* $outdir/native_space/voxelwise
    cp $outdir/native_space/* $outdir/native_space/refregion
    find $outdir/native_space -maxdepth 1 -type f -delete
fi

# Region analysis if requested
if [ ! -z $region_analysis ]; then

    PYTHON="fslpython"
    if ! [ -x "$(command -v fslpython)" ]; then
        Log "WARNING: fslpython not found - using generic python executable"
        PYTHON="python"
    fi

    roi_analysis_args="--add-arrival --gm-thresh $gm_thresh --wm-thresh $wm_thresh"
    roi_analysis_args="$roi_analysis_args --add-standard-atlases"
    if [ -f "$tempdir/gmmask_pure_cort."* ]; then
        # If we have cortical GM mask assume we have WM one too
        roi_analysis_args="$roi_analysis_args --roi-native $tempdir/gmmask_pure_cort $tempdir/wmmask_pure_cereb"
    fi

    for atlas in ${region_analysis_atlases[@]}; do
        Log "Adding custom region analysis atlas: $atlas"
        roi_analysis_add_atlases="$roi_analysis_add_atlases $atlas"
    done
    if [ ! -z "$roi_analysis_add_atlases" ]; then
        roi_analysis_args="$roi_analysis_args --add-mni-atlas $roi_analysis_add_atlases"
    fi

    for atlas_labels in ${region_analysis_atlases_labels[@]}; do
        Log "Using custom atlas labels in: $atlas_labels"
        roi_analysis_add_atlases_labels="$roi_analysis_add_atlases_labels $atlas_labels"
    done
    if [ ! -z "$roi_analysis_add_atlases_labels" ]; then
        roi_analysis_args="$roi_analysis_args --add-mni-atlas-labels $roi_analysis_add_atlases_labels"
    fi

    if [ ! -z $region_analysis_save_rois ]; then
        roi_analysis_args="$roi_analysis_args --save-mni-rois --save-struct-rois --save-native-rois"
    fi

    if [ ! -z $region_analysis_psf ]; then
        roi_analysis_args="$roi_analysis_args --psf=$region_analysis_psf"
    fi

    if [ -z $fslanat ]; then
        if [ -z $struc ]; then
            Log "WARNING: Region analysis requested but no structural data provided"
        elif [ -z $warp ]; then
            Log "WARNING: Region analysis requested but no structural->std space transformation provided"
        else
            if [ $cmethod = 'both' ]; then
                $PYTHON -u `which oxford_asl_roi_stats` --oxasl-output $outdir/native_space/voxelwise --output $outdir/native_space/voxelwise/region_analysis --gm-pve $tempdir/pvgm_struct --wm-pve $tempdir/pvwm_struct --struc $struc --struc2std $warp $roi_analysis_args
                $PYTHON -u `which oxford_asl_roi_stats` --oxasl-output $outdir/native_space/refregion --output $outdir/native_space/refregion/region_analysis --gm-pve $tempdir/pvgm_struct --wm-pve $tempdir/pvwm_struct --struc $struc --struc2std $warp $roi_analysis_args
            else
                $PYTHON -u `which oxford_asl_roi_stats` --oxasl-output $outdir/native_space --output $outdir/native_space/region_analysis --gm-pve $tempdir/pvgm_struct --wm-pve $tempdir/pvwm_struct --struc $struc --struc2std $warp $roi_analysis_args
            fi
        fi
    else
        if [ $cmethod = 'both' ]; then
            $PYTHON -u `which oxford_asl_roi_stats` --oxasl-output $outdir/native_space/voxelwise --output $outdir/native_space/voxelwise/region_analysis --fslanat $fslanat $roi_analysis_args
            $PYTHON -u `which oxford_asl_roi_stats` --oxasl-output $outdir/native_space/refregion --output $outdir/native_space/refregion/region_analysis --fslanat $fslanat $roi_analysis_args
        else
            $PYTHON -u `which oxford_asl_roi_stats` --oxasl-output $outdir/native_space --output $outdir/native_space/region_analysis --fslanat $fslanat $roi_analysis_args
        fi
    fi
fi

# Optional QC outputs - this is a selection of useful intermediate outputs that may be
# useful but not essential enough to be part of the default output, and avoids the
# need to save the entire temporary directory.
if [ ! -z $qc_output ]; then
    qcdir=$outdir/qc_output
    mkdir $qcdir
    Log "Generating QC output in $qcdir"
    
    # Corrected and uncorrected data
    imcp $tempdir/asldata $qcdir/asldata
    imcp $tempdir/asldata_orig $qcdir/asldata_orig
    imcp $tempdir/calib_wholehead $qcdir/calib_wholehead
    imcp $tempdir/calib $qcdir/calib_brain
    imcp $tempdir/calib_orig $qcdir/calib_orig
    imcp $tempdir/diffdata $qcdir/asldata_diff
    imcp $tempdir/diffdata_mean $qcdir/asldata_diff_mean
    
    # Motion parameters
    if [ ! -z $moco ]; then
        cp $tempdir/asldata.par $qcdir/moco_params.par
        fsl_tsplot -i $tempdir/asldata.par \
            -t "MCFLIRT estimated rotations (rads)" \
            -u 1 --start=1 --finish=3 -a x,y,z -w 640 -h 144 \
            -o $qcdir/moco_rotations.png
        fsl_tsplot -i $tempdir/asldata.par \
            -t "MCFLIRT estimated translations (mm)" \
            -u 1 --start=4 --finish=6 -a x,y,z -w 640 -h 144 \
            -o $qcdir/moco_translations.png
    fi

    # Structural registration including original mask
    imcp $tempdir/asl2struc_initial $qcdir/asl2struc_initial
    imcp $tempdir/calib2struc_initial $qcdir/calib2struc_initial
    imcp $tempdir/tmp_asl_reg/low2high_final $qcdir/asl2struc_final
    imcp $tempdir/calib_struc $qcdir/calib2struc_final
    imcp $regfrom $qcdir/regfrom
    imcp $tempdir/mask_orig $qcdir/mask_orig
    imcp $tempdir/mask_dil_orig $qcdir/mask_dil_orig

    # Model fit
    finalstep=`ls -d $tempdir/basil/step? | sed -n '$ p'`
    imcp $finalstep/modelfit $qcdir/modelfit
fi

# clearup
if [ ! -z $debug ]; then
    mv $tempdir $outdir
else
    rm -r $tempdir
fi


Log "Output is $outdir/"
Log "OXFORD_ASL - done."
