#!/bin/bash

# TOAST: TurbO quasar asl bayesian AnalySis Tool
#
# Moss Zhao, Michael Chappell, QUBIC Group, IBME, University of Oxford Image
#
# Copyright (c) 2018-2018 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}

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

Usage() {
    # For future development, I have commented some potentially useful options
    echo "TOAST: TurbO quasar asl bayesian AnalySis Tool"
    echo "Beta version: it only works for AMC's file format:"
    echo "The input file should be one single file and the t dimension should be in the following order:"
    echo "Shift 1 Repeat 1 Tag, Shift 1 Repeat 1 Control, Shift 2 Repeat 1 Tag, Shift 2 Repeat 1 Control."
    echo ""
    echo "Usage (optional parameters in {}):"
    echo "Example (Acquired M0 method in the paper): toast -i <input file> -o <output directory> --struct <High resolution structural data> --calib <M0 file> --tr <TR of M0> --corrcal"
    echo "Example (Estimated M0 method in the paper): toast -i <input file> -o <output directory> --infert1  --corrcal"
    echo "Example (To estimate ABV, use the previous two methods with this additional option): --inferart"
    echo "Example (To correct bolus duration, use the previous two methods with this additional option): --tau <bolus duration from PC MRI data> --infertau"
    echo " -i         : specify data file"
    echo " {-o}       : specify output directory"
    echo " {-m}       : specify brain mask file"
    echo ""
    echo " Extended options:"
    echo " --calib    : include a calibration image (this will automatically tigger the --trans option)"
    echo " --transoff : transform the calibration image to the current data resolution"
    echo " --tr : TR of the calibration image. It needs to be corrected if less than 5.0 sec {default: 5.0 sec}"
    echo " --struct   : structural image (T1 or T2 weighted)"
    echo " --t1       : Set the value for T1 of tissue {default: 1.3 sec}"
    echo " --infert1  : Infer T1 of tissue from Turbo QUASAR ASL control data"
    echo " --t1b      : Set the value for T1 of arterial blood {default: 1.6 sec}"
    echo " --corrcal  : Correct partial volume effects on the edge of calibration image M0a"
    echo " --infertau : estimate bolus duration from data"
    echo " --taulowest: lowest possible bolus duration {default: 0.2 sec}"
    echo " --inferart : estimate arterial blood volume (ABV) from the arterial blood component"
    # Future development
    #echo " --disp     : include bolus dispersion in the model; choices: none(default), gamma"
    # Future development
    #echo " --mfree    : Do model-free (SVD deconvolution) analysis"
    # Future development of partial volume correction options
    #echo ""
    #echo " Partial volume effects correction options:"
    #echo " --pvcorr   : Set partial volume effect correction on. You should provide high resolution partial volume estimates and a structural image."
    #echo " --fslanat  : Name of the directory containing the output from fsl_anat"
    #echo " --t1wm     : T1 for WM {default: 1.1 s}"
    #echo ""
    echo " Sequence parameters:"
    echo " --tau      : Set bolus duration {default: 0.6 s}"
    echo " --slicedt  : Set the increase in TI with slice {default: 0.035 s}"
    echo " --alpha    : Inversion efficiency {default: 0.91}"
    echo " --lambda   : Blood-water coefficient factor {default: 0.9}"
    echo " --fa       : Flip angle for LL readout {default: 35 degrees}"
    echo " --lfa      : Lower flip angle for final phase of data {default: 11.7 degrees}"
    echo " --tis      : comma separated list of TI values"
    echo "            {default: 0.03,0.63,1.23,1.83,2.43,3.03,3.63,4.23,4.83,5.43,6.03}"
    echo " --shift    : slice shifting factor {default: 2}"
    echo " --break_1  : slice number of first acquisition point (start from 0) {default: 0}"
    echo " --break_2  : slice number of middle acquisition point (start from 0) {default: 7}"
    echo " --break_3  : slice number of last acquisition point (start from 0) {default: 14}"
    echo " --taupat   : comma separated list of bolus pattern. 1: label, 0: skip"
    echo "            {default: 1, 1, 1, 1, 1, 1}"
    echo ""
}

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

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

until [ -z $1 ]; do

    # look at this option and determine if has an argument specified by an =
    option=`echo $1 | sed s/=.*//`
    arg="" #specifies if an argument is to be read from next item on command line (=1 is is when = is used)

    if [ $option = $1 ]; then
    # no argument to this command has been found with it (i.e. after an =)
    # if there is an argument it will be the next option
        argument=$2
    else
        arg=1
        argument=`echo $1 | sed s/.*=//`
    fi
    takeargs=0;boolarg="";isbool="";

        case $option in
        -o) outflag=1 outdir=$argument
            shift;;
        -i) inflag=1 infile=$argument #input/data file
            shift;;
        -m) mask=$argument
            shift;;
        --calib) calib_data=$argument
                 trans_status=1 # Need to perform transformation of the calibration image
            shift;;
        --transoff) trans_status=0 # No need to perform transformation of the calibration image
            ;;
        --tr) TR_calib=$argument # TR of the calibration image (default is 5)
            shift;;
        --struct) data_struct=$argument
            shift;;
        --t1b) t1b=$argument
            shift;;
        --t1) t1=$argument
            shift;;
        --infert1) infert1=1
            ;;
        --t1wm) t1wm=$argument
            shift;;
        --slicedt) slicedt=$argument
            shift;;
        --fa) fa=$argument
            shift;;
        --lfa) lfa=$argument
            shift;;
        --shift) shift_factor=$argument
            shift;;
        --break_1) break_1=$argument
            shift;;
        --break_2) break_2=$argument
            shift;;
        --break_3) break_3=$argument
            shift;;
        --disp) disp=$argument
            shift;;
        #--mfree) mfree=1
        #    ;;
        #--edgecorr) isbool=1;
        #   boolarg=edgecorr;
        #   ;;
        --tis) tis=$argument
            shift;;
        #--iform) iform=$argument # to choose the input form of the data
        #    shift;;
        --tau) tau=$argument
            shift;;
        --corrcal) corrcal=1
            ;;
        --alpha) alpha=$argument
            shift;;
        --lambda) lambda=$argument
            shift;;
        --infertau) infertau=1
            ;;
        --taulowest) tau_lowest=$argument
            shift;;
        --inferart) infer_arterial=1
            ;;
        --pvcorr) pvcorr=1
            ;;
        #--pvgm) pvgm_highres=$argument
            #shift;;
        #--pvwm) pvwm_highres=$argument
            #shift;;
        #--s) strim=$argument # structure images for registration and partial volume estimation
            #shift;;
        #--fslanat) fslanat=$argument # Directory containing the output from fsl_anat
        #    shift;;

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

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

        shift
done

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

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

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

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

# save the starting directory
stdir=`pwd`

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

# make a temp directory to work in
#tmpbase=`tmpnam`
#tmpbase="haha"
tmpbase="turquat_temp_"
timestamp=$(date +"%H%M%S")
#tempdir=${tmpbase}_turquat
tempdir=$tmpbase$timestamp
mkdir $tempdir

# deal with the TIs
if [ -z $tis ]; then
    # default QUASAR list of TIs
    tis="0.03,0.63,1.23,1.83,2.43,3.03,3.63,4.23,4.83,5.43,6.03"
fi

# deal with bolus patterns
if [ -z $taupat ]; then
    # default QUASAR list of TIs
    taupat="1,1,1,1,1,1"
fi

# Create TI list
count=0
tislist=""
tislist_calib="" # only the first four TIs
thetis=`echo $tis | sed 's:,: :g'`
for ti in $thetis; do
    count=`expr ${count} + 1`
    tislist=`echo $tislist --ti${count}=$ti`
    if [ "$count" -le "4" ]; then
        tislist_calib=`echo $tislist_calib --ti${count}=$ti`
    fi
done
# echo "TIs list: $tislist" >> $log
ntis=$count;
# TIs for MT and non MT effects

# Hard coded. Can Martin fix this?
tislist_calib_MT="--ti1=0.03 --ti2=0.63 --ti3=1.23 --ti4=1.83 --ti5=2.43 --ti6=3.03"
tislist_calib_non_MT="--ti1=3.63 --ti2=4.23 --ti3=4.83 --ti4=5.34 --ti5=6.03"
TI_start_of_MT=0 # TI = 0. Essentially the last signal with MT effects. Also the starting point of non MT saturation recovery
TI_start_of_non_MT=3.03 # TI = 6. Essentially the last signal with MT effects. Also the starting point of non MT saturation recovery


# Create bolus pattern list
count=0
taulist=""
thetaus=`echo $taupat | sed 's:,: :g'`
for current_taupat in $thetaus; do
    count=`expr ${count} + 1`
    bolus_pattern_list=`echo $bolus_pattern_list --bolus_${count}=$current_taupat`
done
# echo "bolus pattern list: $bolus_pattern_list" >> $log
ntis=$count;


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

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

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

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

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

# calibration parameters
if [ -z $TR_calib ]; then
    TR_calib=5;
fi
# sequence parameters
# slicedt
if [ -z $slicedt ]; then
    slicedt=0.036;
fi
# Flip angles
if [ -z $fa ]; then
    fa=35;
fi
if [ -z $lfa ]; then
    lfa=11.7;
fi
if [ -z $alpha ]; then
    alpha=0.91;
fi
if [ -z $lambda ]; then
    lambda=0.9;
fi

if [ -z $disp ]; then
    disp=none;
fi

# Slice shifting factor
if [ -z $shift_factor ]; then
    shift_factor=1;
fi
if [ -z $break_1 ]; then
    break_1=0;
    break_1_result=$break_1;
fi
if [ -z $break_2 ]; then
    break_2=7;
    break_2_result=`expr $break_2 + 1`;
fi
if [ -z $break_3 ]; then
    break_3=14;
    break_3_result=$break_3;
fi

# Compute segment length
segment_1_length=`expr $break_2 - $break_1`
segment_2_length=`expr $break_3 - $break_2 + 1`
segment_1_length_result=`expr $break_2_result - $break_1_result`
segment_2_length_result=`expr $break_3_result - $break_2_result + 1`


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

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

cd $tempdir
# Re-arrange the file
# Split the input file into different repeats and shifts
shift_1_repeat_1_tag="shift_1_repeat_1_tag"
shift_1_repeat_1_control="shift_1_repeat_1_control"
shift_1_repeat_1_diff="shift_1_repeat_1_diff"
shift_1_repeat_1_tissue="shift_1_repeat_1_tissue"
shift_1_repeat_1_blood="shift_1_repeat_1_blood"

shift_2_repeat_1_tag="shift_2_repeat_1_tag"
shift_2_repeat_1_control="shift_2_repeat_1_control"
shift_2_repeat_1_diff="shift_2_repeat_1_diff"
shift_2_repeat_1_control_effective="shift_2_repeat_1_control_effective"
shift_2_repeat_1_tag_effective="shift_2_repeat_1_tag_effective"
shift_2_repeat_1_diff_effective="shift_2_repeat_1_diff_effective"
mask_shift_2_effective="mask_shift_2_effective"
shift_2_repeat_1_tissue="shift_2_repeat_1_tissue"
shift_2_repeat_1_blood="shift_2_repeat_1_blood"


fslroi data $shift_1_repeat_1_tag 0 77
fslroi data $shift_1_repeat_1_control 77 77
fslroi data $shift_2_repeat_1_tag 154 77
fslroi data $shift_2_repeat_1_control 231 77

# Create a mask for the whole data. We use shift 1 repeat 1 control data
if [ -z $mask ]; then
# auto generate mask
    echo "Creating mask from ASL control data"
    fslmaths $shift_1_repeat_1_control -Tmean aslmean
    bet aslmean mask
    fslmaths mask -bin mask
else
    cd $stdir
    imcp $mask $tempdir/mask
    cd $tempdir
fi
mask=mask

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

# Rearrange the slices for shift 2
# Do control image first
# Shift 2 We need to split and rearrange the z direction
# Control
fslroi $shift_2_repeat_1_control data_segment_1 0 64 0 64 $break_1 $segment_1_length
fslroi $shift_2_repeat_1_control data_segment_2 0 64 0 64 $break_2 $segment_2_length
fslmerge -z $shift_2_repeat_1_control_effective data_segment_2 data_segment_1 # Segment 2 was scanned first   
# Tag
fslroi $shift_2_repeat_1_tag data_segment_1 0 64 0 64 $break_1 $segment_1_length
fslroi $shift_2_repeat_1_tag data_segment_2 0 64 0 64 $break_2 $segment_2_length
fslmerge -z $shift_2_repeat_1_tag_effective data_segment_2 data_segment_1 # Segment 2 was scanned first 
# Mask
fslroi $mask data_segment_1 0 64 0 64 $break_1 $segment_1_length
fslroi $mask data_segment_2 0 64 0 64 $break_2 $segment_2_length
fslmerge -z $mask_shift_2_effective data_segment_2 data_segment_1 # Segment 2 was scanned first


# Take tag control differences
fslmaths $shift_1_repeat_1_tag -sub $shift_1_repeat_1_control $shift_1_repeat_1_diff
fslmaths $shift_2_repeat_1_tag_effective -sub $shift_2_repeat_1_control_effective $shift_2_repeat_1_diff_effective


# Extract Tissue and Arterial blood components
n_dynamics=7
# There are seven dynamics (last one is low flip angle, discard)
# Shift 1
$asl_file --data=$shift_1_repeat_1_diff --ntis=$n_dynamics --ibf=tis --iaf=diff --split=data_diff_ph_
# Tissue component
fslmaths data_diff_ph_000 -add data_diff_ph_001 -add data_diff_ph_003 -add data_diff_ph_004 -div 4 $shift_1_repeat_1_tissue
# Arterial component
fslmaths data_diff_ph_002 -add data_diff_ph_005 -div 2 -sub $shift_1_repeat_1_tissue $shift_1_repeat_1_blood
# Remove the intermediate images
imrm data_diff_ph_000 data_diff_ph_001 data_diff_ph_002 data_diff_ph_003 data_diff_ph_004 data_diff_ph_005 data_diff_ph_006

# Shift 2
$asl_file --data=$shift_2_repeat_1_diff_effective --ntis=$n_dynamics --ibf=tis --iaf=diff --split=data_diff_ph_
# Tissue component
fslmaths data_diff_ph_000 -add data_diff_ph_001 -add data_diff_ph_003 -add data_diff_ph_004 -div 4 $shift_2_repeat_1_tissue
# Arterial component
fslmaths data_diff_ph_002 -add data_diff_ph_005 -div 2 -sub $shift_2_repeat_1_tissue $shift_2_repeat_1_blood
# Remove the intermediate images
imrm data_diff_ph_000 data_diff_ph_001 data_diff_ph_002 data_diff_ph_003 data_diff_ph_004 data_diff_ph_005 data_diff_ph_006



# Here we need to estimate T1 and flip angle correction (g) values using the control data
# We need to split the control data in to MT (first 6 TI) and non MT (last 5 TI) parts
# Now estimate T1 and M0 from Shift 1
$asl_file --data=$shift_1_repeat_1_control --ntis=$n_dynamics --ibf=tis --iaf=diff --split=data_control_

# Split the TIs into MT and non-MT parts. First six TIs have MT, last five TIs no MT
last_TI_1=0
last_TI_2=6
length_1=6
length_2=5
# First six dynamics contains high flip angle signal
fslroi data_control_000 data_control_0_m_1 $last_TI_1 $length_1
fslroi data_control_000 data_control_0_m_2 $last_TI_2 $length_2
fslroi data_control_001 data_control_1_m_1 $last_TI_1 $length_1
fslroi data_control_001 data_control_1_m_2 $last_TI_2 $length_2
fslroi data_control_002 data_control_2_m_1 $last_TI_1 $length_1
fslroi data_control_002 data_control_2_m_2 $last_TI_2 $length_2
fslroi data_control_003 data_control_3_m_1 $last_TI_1 $length_1
fslroi data_control_003 data_control_3_m_2 $last_TI_2 $length_2
fslroi data_control_004 data_control_4_m_1 $last_TI_1 $length_1
fslroi data_control_004 data_control_4_m_2 $last_TI_2 $length_2
fslroi data_control_005 data_control_5_m_1 $last_TI_1 $length_1
fslroi data_control_005 data_control_5_m_2 $last_TI_2 $length_2
fslroi data_control_006 data_control_shift_1_MT_low_fa $last_TI_1 $length_1 # Dynamic 7 contains low flip angle signal
fslroi data_control_006 data_control_shift_1_non_MT_low_fa $last_TI_2 $length_2

echo "Split data into MT and non MT parts..."
# Now take the mean of MT and non-MT part
# You have to use two dynamics (repeats) of the high flip angle data to minimize the MT effects (See Turbo QUASAR paper)
fslmaths data_control_0_m_1 -add data_control_1_m_1 -add data_control_2_m_1 -add data_control_3_m_1 -add data_control_4_m_1 -add data_control_5_m_1 -div 6 data_control_shift_1_MT_high_fa
fslmaths data_control_0_m_2 -add data_control_1_m_2 -add data_control_2_m_2 -add data_control_3_m_2 -add data_control_4_m_2 -add data_control_5_m_2 -div 6 data_control_shift_1_non_MT_high_fa
# Now merge with the low flip angle data
fslmerge -t data_control_shift_1_MT data_control_shift_1_MT_high_fa data_control_shift_1_MT_low_fa
fslmerge -t data_control_shift_1_non_MT data_control_shift_1_non_MT_high_fa data_control_shift_1_non_MT_low_fa

# You also need two initial points for the two part of the data, ie. where the saturation signal recovers from
# For the data with MT effects, the initial signal is zero
# For the data without MT effects, the initial signal is TI = 6
initial_point=5 # TI = 6, starts from zero
initial_length=1 # We only need one point
fslroi data_control_000 data_initial_000 $initial_point $initial_length
fslroi data_control_001 data_initial_001 $initial_point $initial_length
fslroi data_control_002 data_initial_002 $initial_point $initial_length
fslroi data_control_003 data_initial_003 $initial_point $initial_length
fslroi data_control_004 data_initial_004 $initial_point $initial_length
fslroi data_control_005 data_initial_005 $initial_point $initial_length
fslroi data_control_006 data_initial_shift_1_non_MT_low_fa $initial_point $initial_length # The initial signal of the low FA curve without MT effects
# Now compute the mean of the first six dynamics (high FA)
fslmaths data_initial_000 -add data_initial_001 -add data_initial_002 -add data_initial_003 -add data_initial_004 -add data_initial_005 -div 6 data_initial_shift_1_non_MT_high_fa
# Now do the initial signal with MT effects. It should starts from zero but we make it a very low value for easier fitting
fslmaths data_initial_shift_1_non_MT_high_fa -mul 0 -add 0.0000001 -mas mask data_initial_shift_1_MT_high_fa
fslmaths data_initial_shift_1_non_MT_low_fa -mul 0 -add 0.0000001 -mas mask data_initial_shift_1_MT_low_fa

echo "Finished spliting data into MT and non MT parts"

# Merge the first seven TIs
#fslmerge -t shift_1_repeat_1_control_TI_4 data_control_0_m_0 data_control_1_m_0 data_control_2_m_0 data_control_3_m_0 data_control_4_m_0 data_control_5_m_0 data_control_6_m_0

# Remove the intermediate files
#imrm data_control_000 data_control_001 data_control_002 data_control_003 data_control_004 data_control_005 data_control_006
#imrm data_control_0_m_0 data_control_1_m_0 data_control_2_m_0 data_control_3_m_0 data_control_4_m_0 data_control_5_m_0 data_control_6_m_0

# TIs for MT and non MT effects
tislist_calib_MT="--ti1=0.03 --ti2=0.63 --ti3=1.23 --ti4=1.83 --ti5=2.43 --ti6=3.03"
tislist_calib_non_MT="--ti1=3.63 --ti2=4.23 --ti3=4.83 --ti4=5.34 --ti5=6.03"
# Done Moss Edit

# Prepare for a fabber option file
# We estimate these paramters:
# T1, M0t, g(FA correction factor), A(saturation efficiency)
current_options_file="options_calib_shift_1.txt"
echo "Begin T1 and M0 estimation for Shift 1 of the data (non MT part, last 5 TIs)"
echo "# Turbo QUASAR analysis calibration options" >> $current_options_file
echo "--mask=mask" >> $current_options_file
echo "--method=spatialvb" >> $current_options_file
echo "--noise=white" >> $current_options_file
echo "--model=satrecovdualfa" >> $current_options_file
echo "--repeats=1" >> $current_options_file
echo "--phases=1" >> $current_options_file
echo $tislist_calib_non_MT >> $current_options_file
echo "--t1=$t1 --FA=$fa --LFA=$lfa " >> $current_options_file
echo "--slicedt=$slicedt" >> $current_options_file
echo "--t_initial_high_FA=$TI_start_of_non_MT" >> $current_options_file # hard coded value
echo "--t_initial_low_FA=$TI_start_of_non_MT" >> $current_options_file # hard coded value
echo "--print-free-energy" >> $current_options_file
echo "--save-residuals" >> $current_options_file
echo "--save-model-fit" >> $current_options_file
echo "--PSP_byname1=M0t" >> $current_options_file
echo "--PSP_byname1_type=M" >> $current_options_file
echo "--PSP_byname2=T1t" >> $current_options_file
echo "--PSP_byname2_type=N" >> $current_options_file
echo "--PSP_byname3=A" >> $current_options_file
echo "--PSP_byname3_type=N" >> $current_options_file
echo "--PSP_byname4=g" >> $current_options_file
echo "--PSP_byname4_type=N" >> $current_options_file
echo "--PSP_byname5=M0_initial_high_fa" >> $current_options_file
echo "--PSP_byname5_type=I" >> $current_options_file
echo "--PSP_byname5_image=data_initial_shift_1_non_MT_high_fa" >> $current_options_file
echo "--PSP_byname6=M0_initial_low_fa" >> $current_options_file
echo "--PSP_byname6_type=I" >> $current_options_file
echo "--PSP_byname6_image=data_initial_shift_1_non_MT_low_fa" >> $current_options_file

# Now perform model fitting to estimate these parameters
$fabber --data=data_control_shift_1_non_MT --data-order=singlefile --output=output_calib_shift_1 -@ $current_options_file

# Here we need to estimate T1 and flip angle correction (g) values using the control data
# Create a mask for the whole data. We use shift 1 repeat 1 control data
if [ ! -z $infert1 ]; then
    echo "Estimating T1 of tissue from ASL control data..."
    # Now do the non MT part

    # T1, M0t, g(FA correction factor), A(saturation efficiency)
    current_options_file="options_calib_shift_1_MT.txt"
    echo "Begin T1 and M0 estimation for Shift 1 of the data (MT part, first 6 TIs)"
    echo "# Turbo QUASAR analysis calibration options" >> $current_options_file
    echo "--mask=mask" >> $current_options_file
    echo "--method=spatialvb" >> $current_options_file
    echo "--noise=white" >> $current_options_file
    echo "--model=satrecovdualfa" >> $current_options_file
    echo "--repeats=1" >> $current_options_file
    echo "--phases=1" >> $current_options_file
    echo $tislist_calib_MT >> $current_options_file
    echo "--t1=$t1 --FA=$fa --LFA=$lfa " >> $current_options_file
    echo "--slicedt=$slicedt" >> $current_options_file
    echo "--t_initial_high_FA=$TI_start_of_MT" >> $current_options_file # hard coded value
    echo "--t_initial_low_FA=$TI_start_of_MT" >> $current_options_file # hard coded value
    echo "--print-free-energy" >> $current_options_file
    echo "--save-residuals" >> $current_options_file
    echo "--save-model-fit" >> $current_options_file
    echo "--PSP_byname1=M0t" >> $current_options_file
    echo "--PSP_byname1_type=M" >> $current_options_file
    echo "--PSP_byname2=T1t" >> $current_options_file
    echo "--PSP_byname2_type=N" >> $current_options_file
    echo "--PSP_byname3=A" >> $current_options_file
    echo "--PSP_byname3_type=N" >> $current_options_file
    echo "--PSP_byname4=g" >> $current_options_file
    echo "--PSP_byname4_type=N" >> $current_options_file
    echo "--PSP_byname5=M0_initial_high_fa" >> $current_options_file
    echo "--PSP_byname5_type=I" >> $current_options_file
    echo "--PSP_byname5_image=data_initial_shift_1_MT_high_fa" >> $current_options_file
    echo "--PSP_byname6=M0_initial_low_fa" >> $current_options_file
    echo "--PSP_byname6_type=I" >> $current_options_file
    echo "--PSP_byname6_image=data_initial_shift_1_MT_low_fa" >> $current_options_file

    # Now perform model fitting to estimate these parameters
    $fabber --data=data_control_shift_1_MT --data-order=singlefile --output=output_calib_shift_1_MT -@ $current_options_file

else
    echo "Use user provided T1 of tissue in CBF quantification..."
    fslmaths output_calib_shift_1/mean_T1t -mul 0 -add $t1 -mas mask output_calib_shift_1/mean_T1t
fi




# Now estimate T1 and M0 from shift 2
$asl_file --data=$shift_2_repeat_1_control_effective --ntis=$n_dynamics --ibf=tis --iaf=diff --split=data_control_

# Split the TIs into MT and non-MT parts. First six TIs have MT, last five TIs no MT
last_TI_1=0
last_TI_2=6
length_1=6
length_2=5
fslroi data_control_000 data_control_0_m_1 $last_TI_1 $length_1
fslroi data_control_000 data_control_0_m_2 $last_TI_2 $length_2
fslroi data_control_001 data_control_1_m_1 $last_TI_1 $length_1
fslroi data_control_001 data_control_1_m_2 $last_TI_2 $length_2
fslroi data_control_002 data_control_2_m_1 $last_TI_1 $length_1
fslroi data_control_002 data_control_2_m_2 $last_TI_2 $length_2
fslroi data_control_003 data_control_3_m_1 $last_TI_1 $length_1
fslroi data_control_003 data_control_3_m_2 $last_TI_2 $length_2
fslroi data_control_004 data_control_4_m_1 $last_TI_1 $length_1
fslroi data_control_004 data_control_4_m_2 $last_TI_2 $length_2
fslroi data_control_005 data_control_5_m_1 $last_TI_1 $length_1
fslroi data_control_005 data_control_5_m_2 $last_TI_2 $length_2
fslroi data_control_006 data_control_shift_2_MT_low_fa $last_TI_1 $length_1
fslroi data_control_006 data_control_shift_2_non_MT_low_fa $last_TI_2 $length_2

echo "Split data into MT and non MT parts..."
# Now take the mean of MT and non-MT part
# You have to use two dynamics (repeats) of the high flip angle data to minimize the MT effects (See Turbo QUASAR paper)
fslmaths data_control_0_m_1 -add data_control_1_m_1 -add data_control_2_m_1 -add data_control_3_m_1 -add data_control_4_m_1 -add data_control_5_m_1 -div 6 data_control_shift_2_MT_high_fa
fslmaths data_control_0_m_2 -add data_control_1_m_2 -add data_control_2_m_2 -add data_control_3_m_2 -add data_control_4_m_2 -add data_control_5_m_2 -div 6 data_control_shift_2_non_MT_high_fa
# Now merge with the low flip angle data
fslmerge -t data_control_shift_2_MT data_control_shift_2_MT_high_fa data_control_shift_2_MT_low_fa
fslmerge -t data_control_shift_2_non_MT data_control_shift_2_non_MT_high_fa data_control_shift_2_non_MT_low_fa

# You also need two initial points for the two part of the data, ie. where the saturation signal recovers from
# For the data with MT effects, the initial signal is zero
# For the data without MT effects, the initial signal is TI = 6
initial_point=5 # TI = 6, starts from zero
initial_length=1 # We only need one point
fslroi data_control_000 data_initial_000 $initial_point $initial_length
fslroi data_control_001 data_initial_001 $initial_point $initial_length
fslroi data_control_002 data_initial_002 $initial_point $initial_length
fslroi data_control_003 data_initial_003 $initial_point $initial_length
fslroi data_control_004 data_initial_004 $initial_point $initial_length
fslroi data_control_005 data_initial_005 $initial_point $initial_length
fslroi data_control_006 data_initial_shift_2_non_MT_low_fa $initial_point $initial_length # The initial signal of the low FA curve without MT effects
# Now compute the mean of the first six dynamics (high FA)
fslmaths data_initial_000 -add data_initial_001 -add data_initial_002 -add data_initial_003 -add data_initial_004 -add data_initial_005 -div 6 data_initial_shift_2_non_MT_high_fa
# Now do the initial signal with MT effects. It should starts from zero but we make it a very low value for easier fitting
fslmaths data_initial_shift_2_non_MT_high_fa -mul 0 -add 0.0000001 -mas mask_shift_2_effective data_initial_shift_2_MT_high_fa
fslmaths data_initial_shift_2_non_MT_low_fa -mul 0 -add 0.0000001 -mas mask_shift_2_effective data_initial_shift_2_MT_low_fa

echo "Finished spliting data into MT and non MT parts"

# Merge the first seven TIs
#fslmerge -t shift_1_repeat_1_control_TI_4 data_control_0_m_0 data_control_1_m_0 data_control_2_m_0 data_control_3_m_0 data_control_4_m_0 data_control_5_m_0 data_control_6_m_0

# Remove the intermediate files
#imrm data_control_000 data_control_001 data_control_002 data_control_003 data_control_004 data_control_005 data_control_006
#imrm data_control_0_m_0 data_control_1_m_0 data_control_2_m_0 data_control_3_m_0 data_control_4_m_0 data_control_5_m_0 data_control_6_m_0


# Done Moss Edit

# Prepare for a fabber option file
# We estimate these paramters:
# T1, M0t, g(FA correction factor), A(saturation efficiency)
current_options_file="options_calib_shift_2.txt"
echo "Begin T1 and M0 estimation for Shift 2 of the data (non MT part, last 5 TIs)"
echo "# Turbo QUASAR analysis calibration options" >> $current_options_file
echo "--mask=mask_shift_2_effective" >> $current_options_file
echo "--method=spatialvb" >> $current_options_file
echo "--noise=white" >> $current_options_file
echo "--model=satrecovdualfa" >> $current_options_file
echo "--repeats=1" >> $current_options_file
echo "--phases=1" >> $current_options_file
echo $tislist_calib_non_MT >> $current_options_file
echo "--t1=$t1 --FA=$fa --LFA=$lfa " >> $current_options_file
echo "--slicedt=$slicedt" >> $current_options_file
echo "--t_initial_high_FA=$TI_start_of_non_MT" >> $current_options_file # hard coded value
echo "--t_initial_low_FA=$TI_start_of_non_MT" >> $current_options_file # hard coded value
echo "--print-free-energy" >> $current_options_file
echo "--save-residuals" >> $current_options_file
echo "--save-model-fit" >> $current_options_file
echo "--PSP_byname1=M0t" >> $current_options_file
echo "--PSP_byname1_type=M" >> $current_options_file
echo "--PSP_byname2=T1t" >> $current_options_file
echo "--PSP_byname2_type=N" >> $current_options_file
echo "--PSP_byname3=A" >> $current_options_file
echo "--PSP_byname3_type=N" >> $current_options_file
echo "--PSP_byname4=g" >> $current_options_file
echo "--PSP_byname4_type=N" >> $current_options_file
echo "--PSP_byname5=M0_initial_high_fa" >> $current_options_file
echo "--PSP_byname5_type=I" >> $current_options_file
echo "--PSP_byname5_image=data_initial_shift_2_non_MT_high_fa" >> $current_options_file
echo "--PSP_byname6=M0_initial_low_fa" >> $current_options_file
echo "--PSP_byname6_type=I" >> $current_options_file
echo "--PSP_byname6_image=data_initial_shift_2_non_MT_low_fa" >> $current_options_file

# Now perform model fitting to estimate these parameters
$fabber --data=data_control_shift_2_non_MT --data-order=singlefile --output=output_calib_shift_2 -@ $current_options_file

# Here we need to estimate T1 and flip angle correction (g) values using the control data
# Create a mask for the whole data. We use shift 1 repeat 1 control data
if [ ! -z $infert1 ]; then
    echo "Estimating T1 of tissue from ASL control data..."
    # Now do the non MT part

    # T1, M0t, g(FA correction factor), A(saturation efficiency)
    current_options_file="options_calib_shift_2_MT.txt"
    echo "Begin T1 and M0 estimation for Shift 2 of the data (MT part, first 5 TIs)"
    echo "# Turbo QUASAR analysis calibration options" >> $current_options_file
    echo "--mask=mask_shift_2_effective" >> $current_options_file
    echo "--method=spatialvb" >> $current_options_file
    echo "--noise=white" >> $current_options_file
    echo "--model=satrecovdualfa" >> $current_options_file
    echo "--repeats=1" >> $current_options_file
    echo "--phases=1" >> $current_options_file
    echo $tislist_calib_MT >> $current_options_file
    echo "--t1=$t1 --FA=$fa --LFA=$lfa " >> $current_options_file
    echo "--slicedt=$slicedt" >> $current_options_file
    echo "--t_initial_high_FA=$TI_start_of_MT" >> $current_options_file # hard coded value
    echo "--t_initial_low_FA=$TI_start_of_MT" >> $current_options_file # hard coded value
    echo "--print-free-energy" >> $current_options_file
    echo "--save-residuals" >> $current_options_file
    echo "--save-model-fit" >> $current_options_file
    echo "--PSP_byname1=M0t" >> $current_options_file
    echo "--PSP_byname1_type=M" >> $current_options_file
    echo "--PSP_byname2=T1t" >> $current_options_file
    echo "--PSP_byname2_type=N" >> $current_options_file
    echo "--PSP_byname3=A" >> $current_options_file
    echo "--PSP_byname3_type=N" >> $current_options_file
    echo "--PSP_byname4=g" >> $current_options_file
    echo "--PSP_byname4_type=N" >> $current_options_file
    echo "--PSP_byname5=M0_initial_high_fa" >> $current_options_file
    echo "--PSP_byname5_type=I" >> $current_options_file
    echo "--PSP_byname5_image=data_initial_shift_2_MT_high_fa" >> $current_options_file
    echo "--PSP_byname6=M0_initial_low_fa" >> $current_options_file
    echo "--PSP_byname6_type=I" >> $current_options_file
    echo "--PSP_byname6_image=data_initial_shift_2_MT_low_fa" >> $current_options_file

    # Now perform model fitting to estimate these parameters
    $fabber --data=data_control_shift_2_MT --data-order=singlefile --output=output_calib_shift_2_MT -@ $current_options_file

else
    echo "Use user provided T1 of tissue in CBF quantification..."
    fslmaths output_calib_shift_2/mean_T1t -mul 0 -add $t1 -mas mask output_calib_shift_2/mean_T1t
fi
# Done Moss edit

# Done estimating T1 and M0


# Model based analysis
# Prepare for a fabber option file
if [ ! -z $infert1 ]; then

    # First we copy the original images
    imcp $shift_1_repeat_1_tissue shift_1_repeat_1_tissue_original
    imcp $shift_2_repeat_1_tissue shift_2_repeat_1_tissue_original
    imcp $shift_1_repeat_1_blood shift_1_repeat_1_blood_original
    imcp $shift_2_repeat_1_blood shift_2_repeat_1_blood_original
    # Get the tissue with MT signals
    fslroi $shift_1_repeat_1_tissue shift_1_repeat_1_tissue_MT $last_TI_1 $length_1
    fslroi $shift_2_repeat_1_tissue shift_2_repeat_1_tissue_MT $last_TI_1 $length_1
    fslroi $shift_1_repeat_1_blood shift_1_repeat_1_blood_MT $last_TI_1 $length_1
    fslroi $shift_2_repeat_1_blood shift_2_repeat_1_blood_MT $last_TI_1 $length_1
    shift_1_repeat_1_tissue_MT=shift_1_repeat_1_tissue_MT
    shift_2_repeat_1_tissue_MT=shift_2_repeat_1_tissue_MT
    shift_1_repeat_1_blood_MT=shift_1_repeat_1_blood_MT
    shift_2_repeat_1_blood_MT=shift_2_repeat_1_blood_MT
    # Get the tissue with non MT signals
    fslroi $shift_1_repeat_1_tissue $shift_1_repeat_1_tissue $last_TI_2 $length_2
    fslroi $shift_2_repeat_1_tissue $shift_2_repeat_1_tissue $last_TI_2 $length_2
    fslroi $shift_1_repeat_1_blood $shift_1_repeat_1_blood $last_TI_2 $length_2
    fslroi $shift_2_repeat_1_blood $shift_2_repeat_1_blood $last_TI_2 $length_2

    # Use only TIs with MT effects
    tislist_MT=$tislist_calib_MT
    tislist=$tislist_calib_non_MT
    #tislist_MT=$tislist_calib_MT
    # Priors
    calib_shift_1_dir_MT=output_calib_shift_1_MT
    calib_shift_2_dir_MT=output_calib_shift_2_MT
    calib_shift_1_dir=output_calib_shift_1
    calib_shift_2_dir=output_calib_shift_2
    # We also need a T1 of blood under MT effects
    # We assume this relationship T1_tissue/T1_tissue_MT = T1_blood/T1_blood_MT
    # So T1_blood_MT = T1_tissue_MT * T1_blood / T1_tissue
    #fslmaths output_calib_shift_1_MT/mean_T1t -mul $t1b -div $t1 -mas mask output_calib_shift_1_MT/mean_T1b
    #fslmaths output_calib_shift_2_MT/mean_T1t -mul $t1b -div $t1 -mas mask_shift_2_effective output_calib_shift_2_MT/mean_T1b
    
    #fslmaths output_calib_shift_1_MT/mean_T1t -mul $t1b -div output_calib_shift_1/mean_T1t -mas mask output_calib_shift_1_MT/mean_T1b
    #fslmaths output_calib_shift_2_MT/mean_T1t -mul $t1b -div output_calib_shift_2/mean_T1t -mas mask_shift_2_effective output_calib_shift_2_MT/mean_T1b

# Use full data simple quantification
else
    # Priors
    calib_shift_1_dir=output_calib_shift_1
    calib_shift_2_dir=output_calib_shift_2
fi

echo "Begin model-based analysis"

    # Fixed bolus duration
    # Condition not to infer bolus duration
    # Shift 1
    current_options_file="options_tissue_shift_1.txt"
    echo "# Turbo QUASAR ASL tissue component analysis options for Shift 1 of the data" > $current_options_file
    echo "--mask=mask" >> $current_options_file
    echo "--method=spatialvb" >> $current_options_file
    echo "--noise=white" >> $current_options_file
    echo "--model=turboquasar" >> $current_options_file
    echo "--disp=$disp" >> $current_options_file
    #echo "--inferart" >> $current_options_file
    echo "--repeats=1" >> $current_options_file
    echo $tislist >> $current_options_file
    echo "--t1=$t1 --t1b=$t1b --t1wm=$t1wm --tau=$tau --fa=$fa " >> $current_options_file
    #if [ ! -z $infert1 ]; then
    #    #echo "--t1_MT=$t1 --t1b_MT=$t1b " >> $current_options_file
    #fi
    echo "--slicedt=$slicedt" >> $current_options_file
    echo $bolus_pattern_list >> $current_options_file
    echo "--slice_shift=$shift_factor" >> $current_options_file
    echo "--onephase" >> $current_options_file
    if [ ! -z $infert1 ]; then
        echo "--usecalib " >> $current_options_file # We incorporate previously estimated T1 and g value as priors
        echo "--infert1" >> $current_options_file # If we infer T1 we must infer T1 of tissue and blood together (must be two of them)
    fi
    #echo "--artdir" >> $current_options_file
    # Save model fitting results and residue
    echo "--save-model-fit" >> $current_options_file
    echo "--print-free-energy" >> $current_options_file
    echo "--save-residuals" >> $current_options_file
    # Here we create a shortcut to the latest results directory
    #echo "--link-to-latest" >> $current_options_file

    # Here are the settings
    echo "--PSP_byname1=ftiss" >> $current_options_file # Estimate CBF
    echo "--PSP_byname1_type=M" >> $current_options_file
    echo "--PSP_byname2=delttiss" >> $current_options_file # Estimate ATT
    echo "--PSP_byname2_type=N" >> $current_options_file
    echo "--PSP_byname3=sp_log" >> $current_options_file # Estimate sp log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname3_type=N" >> $current_options_file
    echo "--PSP_byname4=s_log" >> $current_options_file # Estimate s log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname4_type=N" >> $current_options_file
    if [ ! -z $infert1 ]; then
        echo "--PSP_byname5=g" >> $current_options_file # Estimate g (FA correction parameter)
        echo "--PSP_byname5_type=I" >> $current_options_file
        echo "--PSP_byname5_image=output_calib_shift_1/mean_g" >> $current_options_file
        echo "--PSP_byname6=T_1" >> $current_options_file # Estimate T1 of tissue
        echo "--PSP_byname6_type=I" >> $current_options_file
        echo "--PSP_byname6_image=output_calib_shift_1/mean_T1t" >> $current_options_file
        #echo "--PSP_byname7=T_1b" >> $current_options_file # Estimate T1 of blood
        #echo "--PSP_byname7_type=I" >> $current_options_file
        #echo "--PSP_byname7_image=output_calib_shift_1/mean_T1b" >> $current_options_file

        #echo "--PSP_byname7=T_1_MT" >> $current_options_file # Estimate T1 of tissue
        #echo "--PSP_byname7_type=I" >> $current_options_file
        #echo "--PSP_byname7_image=output_calib_shift_1_MT/mean_T1t" >> $current_options_file
        #echo "--PSP_byname8=T_1b_MT" >> $current_options_file # Estimate T1 of blood
        #echo "--PSP_byname8_type=I" >> $current_options_file
        #echo "--PSP_byname8_image=output_calib_shift_1_MT/mean_T1b" >> $current_options_file
    fi

    # Now perform model fitting to estimate CBF, ATT, etc
    $fabber --data=$shift_1_repeat_1_tissue --data-order=singlefile --output=output_tissue_shift_1 -@ $current_options_file
    


    # Shift 2
    # Fixed bolus duration
    # Condition not to infer bolus duration
    current_options_file="options_tissue_shift_2.txt"
    echo "# Turbo QUASAR ASL tissue component analysis options for Shift 2 of the data" > $current_options_file
    echo "--mask=mask_shift_2_effective" >> $current_options_file
    echo "--method=spatialvb" >> $current_options_file
    echo "--noise=white" >> $current_options_file
    echo "--model=turboquasar" >> $current_options_file
    echo "--disp=$disp" >> $current_options_file
    #echo "--inferart" >> $current_options_file
    echo "--repeats=1" >> $current_options_file
    echo $tislist >> $current_options_file
    echo "--t1=$t1 --t1b=$t1b --t1wm=$t1wm --tau=$tau --fa=$fa " >> $current_options_file
    #if [ ! -z $infert1 ]; then
    #    #echo "--t1_MT=$t1 --t1b_MT=$t1b " >> $current_options_file
    #fi
    echo "--slicedt=$slicedt" >> $current_options_file
    echo $bolus_pattern_list >> $current_options_file
    echo "--slice_shift=$shift_factor" >> $current_options_file
    echo "--onephase" >> $current_options_file
    if [ ! -z $infert1 ]; then
        echo "--usecalib " >> $current_options_file # We incorporate previously estimated T1 and g value as priors
        echo "--infert1" >> $current_options_file # If we infer T1 we must infer T1 of tissue and blood together (must be two of them)
    fi
    #echo "--artdir" >> $current_options_file
    # Save model fitting results and residue
    echo "--save-model-fit" >> $current_options_file
    echo "--print-free-energy" >> $current_options_file
    echo "--save-residuals" >> $current_options_file
    # Here we create a shortcut to the latest results directory
    #echo "--link-to-latest" >> $current_options_file

    # Here are the settings
    echo "--PSP_byname1=ftiss" >> $current_options_file # Estimate CBF
    echo "--PSP_byname1_type=M" >> $current_options_file
    echo "--PSP_byname2=delttiss" >> $current_options_file # Estimate ATT
    echo "--PSP_byname2_type=N" >> $current_options_file
    echo "--PSP_byname3=sp_log" >> $current_options_file # Estimate sp log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname3_type=N" >> $current_options_file
    echo "--PSP_byname4=s_log" >> $current_options_file # Estimate s log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname4_type=N" >> $current_options_file
    if [ ! -z $infert1 ]; then
        echo "--PSP_byname5=g" >> $current_options_file # Estimate g (FA correction parameter)
        echo "--PSP_byname5_type=I" >> $current_options_file
        echo "--PSP_byname5_image=output_calib_shift_2/mean_g" >> $current_options_file
        echo "--PSP_byname6=T_1" >> $current_options_file # Estimate T1 of tissue
        echo "--PSP_byname6_type=I" >> $current_options_file
        echo "--PSP_byname6_image=output_calib_shift_2/mean_T1t" >> $current_options_file
        #echo "--PSP_byname7=T_1b" >> $current_options_file # Estimate T1 of blood
        #echo "--PSP_byname7_type=I" >> $current_options_file
        #echo "--PSP_byname7_image=output_calib_shift_2/mean_T1b" >> $current_options_file

        #echo "--PSP_byname7=T_1_MT" >> $current_options_file # Estimate T1 of tissue
        #echo "--PSP_byname7_type=I" >> $current_options_file
        #echo "--PSP_byname7_image=output_calib_shift_2_MT/mean_T1t" >> $current_options_file
        #echo "--PSP_byname8=T_1b_MT" >> $current_options_file # Estimate T1 of blood
        #echo "--PSP_byname8_type=I" >> $current_options_file
        #echo "--PSP_byname8_image=output_calib_shift_2_MT/mean_T1b" >> $current_options_file
    fi


    # Now perform model fitting to estimate CBF, ATT, etc
    $fabber --data=$shift_2_repeat_1_tissue --data-order=singlefile --output=output_tissue_shift_2 -@ $current_options_file

    # Condition to estimate bolus duration
    if [ ! -z $infertau ]; then

        if [ -z $tau_lowest ]; then
            tau_lowest=0.2;
        fi
        # We need two steps to estimate it
        # Step 1: Estimate CBF and ATT (as the previous step)
        # Step 2: Incorporate results from Step 1 to estimate bolus duration and refine CBF and ATT 
        echo "Inferring bolus duration..."

        # Shift 1
        echo "Estimating bolus duration from Shift 1 of the data..."
        # Rename the results of the previous step
        mv output_tissue_shift_1 output_tissue_shift_1_step_1
        # Now we need to create a new MVN file to include the new bolus duration parameter
        mvntool --input=output_tissue_shift_1_step_1/finalMVN --output=output_tissue_shift_1_step_1/finalMVN2 --mask=mask --param=3 --new --val=1 --var=1

        # Now create the new option files
        current_options_file="options_tissue_shift_1.txt"
        echo "--infertau" >> $current_options_file # Estimate bolus duration
        echo "--tau_lowest=$tau_lowest" >> $current_options_file # Lowest limit of bolus duration
        if [ ! -z $infert1 ]; then
            echo "--PSP_byname7=tautiss" >> $current_options_file # Estimate bolus duration
            echo "--PSP_byname7_type=N" >> $current_options_file
        else
            echo "--PSP_byname5=tautiss" >> $current_options_file # Estimate bolus duration
            echo "--PSP_byname5_type=N" >> $current_options_file
        fi


        $fabber --data=$shift_1_repeat_1_tissue --data-order=singlefile --output=output_tissue_shift_1 -@ $current_options_file --continue-from-mvn=output_tissue_shift_1_step_1/finalMVN2
    
        # Shift 2
        echo "Estimating bolus duration from Shift 1 of the data..."
        # Rename the results of the previous step
        mv output_tissue_shift_2 output_tissue_shift_2_step_1
        # Now we need to create a new MVN file to include the new bolus duration parameter
        mvntool --input=output_tissue_shift_2_step_1/finalMVN --output=output_tissue_shift_2_step_1/finalMVN2 --mask=mask_shift_2_effective --param=3 --new --val=1 --var=1

        # Now create the new option files
        current_options_file="options_tissue_shift_2.txt"
        echo "--infertau" >> $current_options_file # Estimate bolus duration
        echo "--tau_lowest=$tau_lowest" >> $current_options_file # Lowest limit of bolus duration
        if [ ! -z $infert1 ]; then
            echo "--PSP_byname7=tautiss" >> $current_options_file # Estimate bolus duration
            echo "--PSP_byname7_type=N" >> $current_options_file
        else
            echo "--PSP_byname5=tautiss" >> $current_options_file # Estimate bolus duration
            echo "--PSP_byname5_type=N" >> $current_options_file
        fi

        $fabber --data=$shift_2_repeat_1_tissue --data-order=singlefile --output=output_tissue_shift_2 -@ $current_options_file --continue-from-mvn=output_tissue_shift_2_step_1/finalMVN2


    fi


# Here we correct MT effects for the tissue component
if [ ! -z $infert1 ]; then

    echo "Correcting for MT effects in tissue data..."

    # Fixed bolus duration
    # Condition not to infer bolus duration
    # Shift 1
    current_options_file="options_tissue_shift_1_MT.txt"
    echo "# Turbo QUASAR ASL tissue component analysis options for Shift 1 of the data with MT effects" > $current_options_file
    echo "--mask=mask" >> $current_options_file
    echo "--method=spatialvb" >> $current_options_file
    echo "--noise=white" >> $current_options_file
    echo "--model=turboquasar" >> $current_options_file
    echo "--disp=$disp" >> $current_options_file
    #echo "--inferart" >> $current_options_file
    echo "--repeats=1" >> $current_options_file
    echo $tislist_calib_MT >> $current_options_file
    echo "--t1=$t1 --t1b=$t1b --t1wm=$t1wm --tau=$tau --fa=$fa " >> $current_options_file
    #if [ ! -z $infert1 ]; then
    #    #echo "--t1_MT=$t1 --t1b_MT=$t1b " >> $current_options_file
    #fi
    echo "--slicedt=$slicedt" >> $current_options_file
    echo $bolus_pattern_list >> $current_options_file
    echo "--slice_shift=$shift_factor" >> $current_options_file
    echo "--onephase" >> $current_options_file
    if [ ! -z $infert1 ]; then
        echo "--usecalib " >> $current_options_file # We incorporate previously estimated T1 and g value as priors
        echo "--infert1" >> $current_options_file # If we infer T1 we must infer T1 of tissue and blood together (must be two of them)
    fi
    #echo "--artdir" >> $current_options_file
    # Save model fitting results and residue
    echo "--save-model-fit" >> $current_options_file
    echo "--print-free-energy" >> $current_options_file
    echo "--save-residuals" >> $current_options_file
    # Here we create a shortcut to the latest results directory
    #echo "--link-to-latest" >> $current_options_file

    # Here are the settings
    echo "--PSP_byname1=ftiss" >> $current_options_file # Estimate CBF
    echo "--PSP_byname1_type=M" >> $current_options_file
    echo "--PSP_byname2=delttiss" >> $current_options_file # Estimate ATT
    echo "--PSP_byname2_type=N" >> $current_options_file
    echo "--PSP_byname3=sp_log" >> $current_options_file # Estimate sp log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname3_type=N" >> $current_options_file
    echo "--PSP_byname4=s_log" >> $current_options_file # Estimate s log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname4_type=N" >> $current_options_file
    if [ ! -z $infert1 ]; then
        echo "--PSP_byname5=g" >> $current_options_file # Estimate g (FA correction parameter)
        echo "--PSP_byname5_type=I" >> $current_options_file
        echo "--PSP_byname5_image=output_calib_shift_1_MT/mean_g" >> $current_options_file
        echo "--PSP_byname6=T_1" >> $current_options_file # Estimate T1 of tissue
        echo "--PSP_byname6_type=I" >> $current_options_file
        echo "--PSP_byname6_image=output_calib_shift_1_MT/mean_T1t" >> $current_options_file
        #echo "--PSP_byname7=T_1b" >> $current_options_file # Estimate T1 of blood
        #echo "--PSP_byname7_type=I" >> $current_options_file
        #echo "--PSP_byname7_image=output_calib_shift_1/mean_T1b" >> $current_options_file

        #echo "--PSP_byname7=T_1_MT" >> $current_options_file # Estimate T1 of tissue
        #echo "--PSP_byname7_type=I" >> $current_options_file
        #echo "--PSP_byname7_image=output_calib_shift_1_MT/mean_T1t" >> $current_options_file
        #echo "--PSP_byname8=T_1b_MT" >> $current_options_file # Estimate T1 of blood
        #echo "--PSP_byname8_type=I" >> $current_options_file
        #echo "--PSP_byname8_image=output_calib_shift_1_MT/mean_T1b" >> $current_options_file
    fi

    # Now perform model fitting to estimate CBF, ATT, etc
    echo $shift_1_repeat_1_tissue_MT
    echo $current_options_file
    $fabber --data=$shift_1_repeat_1_tissue_MT --data-order=singlefile --output=output_tissue_shift_1_MT -@ $current_options_file
    


    # Shift 2
    # Fixed bolus duration
    # Condition not to infer bolus duration
    current_options_file="options_tissue_shift_2_MT.txt"
    echo "# Turbo QUASAR ASL tissue component analysis options for Shift 2 of the data with MT effects" > $current_options_file
    echo "--mask=mask_shift_2_effective" >> $current_options_file
    echo "--method=spatialvb" >> $current_options_file
    echo "--noise=white" >> $current_options_file
    echo "--model=turboquasar" >> $current_options_file
    echo "--disp=$disp" >> $current_options_file
    #echo "--inferart" >> $current_options_file
    echo "--repeats=1" >> $current_options_file
    echo $tislist_calib_MT >> $current_options_file
    echo "--t1=$t1 --t1b=$t1b --t1wm=$t1wm --tau=$tau --fa=$fa " >> $current_options_file
    #if [ ! -z $infert1 ]; then
    #    #echo "--t1_MT=$t1 --t1b_MT=$t1b " >> $current_options_file
    #fi
    echo "--slicedt=$slicedt" >> $current_options_file
    echo $bolus_pattern_list >> $current_options_file
    echo "--slice_shift=$shift_factor" >> $current_options_file
    echo "--onephase" >> $current_options_file
    if [ ! -z $infert1 ]; then
        echo "--usecalib " >> $current_options_file # We incorporate previously estimated T1 and g value as priors
        echo "--infert1" >> $current_options_file # If we infer T1 we must infer T1 of tissue and blood together (must be two of them)
    fi
    #echo "--artdir" >> $current_options_file
    # Save model fitting results and residue
    echo "--save-model-fit" >> $current_options_file
    echo "--print-free-energy" >> $current_options_file
    echo "--save-residuals" >> $current_options_file
    # Here we create a shortcut to the latest results directory
    #echo "--link-to-latest" >> $current_options_file

    # Here are the settings
    echo "--PSP_byname1=ftiss" >> $current_options_file # Estimate CBF
    echo "--PSP_byname1_type=M" >> $current_options_file
    echo "--PSP_byname2=delttiss" >> $current_options_file # Estimate ATT
    echo "--PSP_byname2_type=N" >> $current_options_file
    echo "--PSP_byname3=sp_log" >> $current_options_file # Estimate sp log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname3_type=N" >> $current_options_file
    echo "--PSP_byname4=s_log" >> $current_options_file # Estimate s log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname4_type=N" >> $current_options_file
    if [ ! -z $infert1 ]; then
        echo "--PSP_byname5=g" >> $current_options_file # Estimate g (FA correction parameter)
        echo "--PSP_byname5_type=I" >> $current_options_file
        echo "--PSP_byname5_image=output_calib_shift_2_MT/mean_g" >> $current_options_file
        echo "--PSP_byname6=T_1" >> $current_options_file # Estimate T1 of tissue
        echo "--PSP_byname6_type=I" >> $current_options_file
        echo "--PSP_byname6_image=output_calib_shift_2_MT/mean_T1t" >> $current_options_file
        #echo "--PSP_byname7=T_1b" >> $current_options_file # Estimate T1 of blood
        #echo "--PSP_byname7_type=I" >> $current_options_file
        #echo "--PSP_byname7_image=output_calib_shift_2/mean_T1b" >> $current_options_file

        #echo "--PSP_byname7=T_1_MT" >> $current_options_file # Estimate T1 of tissue
        #echo "--PSP_byname7_type=I" >> $current_options_file
        #echo "--PSP_byname7_image=output_calib_shift_2_MT/mean_T1t" >> $current_options_file
        #echo "--PSP_byname8=T_1b_MT" >> $current_options_file # Estimate T1 of blood
        #echo "--PSP_byname8_type=I" >> $current_options_file
        #echo "--PSP_byname8_image=output_calib_shift_2_MT/mean_T1b" >> $current_options_file
    fi


    # Now perform model fitting to estimate CBF, ATT, etc
    $fabber --data=$shift_2_repeat_1_tissue_MT --data-order=singlefile --output=output_tissue_shift_2_MT -@ $current_options_file

    # Condition to estimate bolus duration
    if [ ! -z $infertau ]; then

        if [ -z $tau_lowest ]; then
            tau_lowest=0.2;
        fi
        # We need two steps to estimate it
        # Step 1: Estimate CBF and ATT (as the previous step)
        # Step 2: Incorporate results from Step 1 to estimate bolus duration and refine CBF and ATT 
        echo "Inferring bolus duration..."

        # Shift 1
        echo "Estimating bolus duration from Shift 1 of the data with MT effects..."
        # Rename the results of the previous step
        mv output_tissue_shift_1_MT output_tissue_shift_1_MT_step_1
        # Now we need to create a new MVN file to include the new bolus duration parameter
        mvntool --input=output_tissue_shift_1_MT_step_1/finalMVN --output=output_tissue_shift_1_MT_step_1/finalMVN2 --mask=mask --param=3 --new --val=1 --var=1

        # Now create the new option files
        current_options_file="options_tissue_shift_1_MT.txt"
        echo "--infertau" >> $current_options_file # Estimate bolus duration
        echo "--tau_lowest=$tau_lowest" >> $current_options_file # Lowest limit of bolus duration
        if [ ! -z $infert1 ]; then
            echo "--PSP_byname7=tautiss" >> $current_options_file # Estimate bolus duration
            echo "--PSP_byname7_type=N" >> $current_options_file
        else
            echo "--PSP_byname5=tautiss" >> $current_options_file # Estimate bolus duration
            echo "--PSP_byname5_type=N" >> $current_options_file
        fi


        $fabber --data=$shift_1_repeat_1_tissue_MT --data-order=singlefile --output=output_tissue_shift_1_MT -@ $current_options_file --continue-from-mvn=output_tissue_shift_1_step_1/finalMVN2
    
        # Shift 2
        echo "Estimating bolus duration from Shift 2 of the data with MT effects..."
        # Rename the results of the previous step
        mv output_tissue_shift_2_MT output_tissue_shift_2_MT_step_1
        # Now we need to create a new MVN file to include the new bolus duration parameter
        mvntool --input=output_tissue_shift_2_MT_step_1/finalMVN --output=output_tissue_shift_2_MT_step_1/finalMVN2 --mask=mask_shift_2_effective --param=3 --new --val=1 --var=1

        # Now create the new option files
        current_options_file="options_tissue_shift_2_MT.txt"
        echo "--infertau" >> $current_options_file # Estimate bolus duration
        echo "--tau_lowest=$tau_lowest" >> $current_options_file # Lowest limit of bolus duration
        if [ ! -z $infert1 ]; then
            echo "--PSP_byname7=tautiss" >> $current_options_file # Estimate bolus duration
            echo "--PSP_byname7_type=N" >> $current_options_file
        else
            echo "--PSP_byname5=tautiss" >> $current_options_file # Estimate bolus duration
            echo "--PSP_byname5_type=N" >> $current_options_file
        fi

        $fabber --data=$shift_2_repeat_1_tissue_MT --data-order=singlefile --output=output_tissue_shift_2_MT -@ $current_options_file --continue-from-mvn=output_tissue_shift_2_step_1/finalMVN2


    fi

fi
# Done model based analysis


# Model free analysis (to be implemented)

# Done model free analysis


# Estimate ABV and the arterial blood components
if [ ! -z $infer_arterial ]; then
    echo "Estimateing ABV and arrival time of arterial blood (macrovasculature)"
    # Prepare a fabber option file
    # Shift 1
    current_options_file="options_blood_shift_1.txt"
    echo "# Turbo QUASAR ASL arterial blood component analysis options for Shift 1 of the data" > $current_options_file
    echo "--mask=mask" >> $current_options_file
    echo "--method=spatialvb" >> $current_options_file
    echo "--noise=white" >> $current_options_file
    echo "--model=turboquasar" >> $current_options_file
    echo "--disp=$disp" >> $current_options_file
    echo "--inferart" >> $current_options_file
    echo "--artdir" >> $current_options_file
    echo "--tissoff" >> $current_options_file
    echo "--repeats=1" >> $current_options_file
    echo $tislist >> $current_options_file
    echo "--t1=$t1 --t1b=$t1b --tau=$tau --fa=$fa " >> $current_options_file
    #if [ ! -z $infert1 ]; then
    #    #echo "--t1_MT=$t1 --t1b_MT=$t1b " >> $current_options_file
    #fi
    echo "--slicedt=$slicedt" >> $current_options_file
    echo $bolus_pattern_list >> $current_options_file
    echo "--slice_shift=$shift_factor" >> $current_options_file
    echo "--onephase" >> $current_options_file
    if [ ! -z $infert1 ]; then
        echo "--usecalib " >> $current_options_file # We incorporate previously estimated T1 and g value as priors
    fi
    #echo "--infert1" >> $current_options_file # If we infer T1 we must infer T1 of tissue and blood together (must be two of them)
    # Save model fitting results and residue
    echo "--save-model-fit" >> $current_options_file
    echo "--print-free-energy" >> $current_options_file
    echo "--save-residuals" >> $current_options_file
    # Here we create a shortcut to the latest results directory
    #echo "--link-to-latest" >> $current_options_file

    # Here are the settings
    echo "--PSP_byname1=fblood" >> $current_options_file # Estimate CBF
    echo "--PSP_byname1_type=M" >> $current_options_file
    echo "--PSP_byname2=deltblood" >> $current_options_file # Estimate ATT
    echo "--PSP_byname2_type=N" >> $current_options_file
    echo "--PSP_byname3=sp_log" >> $current_options_file # Estimate sp log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname3_type=N" >> $current_options_file
    echo "--PSP_byname4=s_log" >> $current_options_file # Estimate s log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname4_type=N" >> $current_options_file
    echo "--PSP_byname5=thblood" >> $current_options_file # Estimate angle theta of crushing gradient
    echo "--PSP_byname5_type=N" >> $current_options_file
    echo "--PSP_byname6=phiblood" >> $current_options_file # Estimate angle phi of crushing gradient
    echo "--PSP_byname6_type=N" >> $current_options_file
    echo "--PSP_byname7=bvblood" >> $current_options_file # Estimate volume of blood
    echo "--PSP_byname7_type=N" >> $current_options_file
    if [ ! -z $infert1 ]; then
        echo "--PSP_byname8=g" >> $current_options_file # Incorporate FA correction factor g
        echo "--PSP_byname8_type=I" >> $current_options_file
        echo "--PSP_byname8_image=output_calib_shift_1/mean_g" >> $current_options_file
        #echo "--PSP_byname9=T_1b_MT" >> $current_options_file # T1 blood under MT effects
        #echo "--PSP_byname9_type=I" >> $current_options_file
        #echo "--PSP_byname9_image=output_calib_shift_1_MT/mean_T1b" >> $current_options_file
    fi

    # Now perform model fitting to estimate ABV and arrival time to macrovasculature
    $fabber --data=$shift_1_repeat_1_blood --data-order=singlefile --output=output_blood_shift_1 -@ $current_options_file
    
    # Shift 2 data
    current_options_file="options_blood_shift_2.txt"
    echo "# Turbo QUASAR ASL arterial blood component analysis options for Shift 2 of the data" > $current_options_file
    echo "--mask=mask_shift_2_effective" >> $current_options_file
    echo "--method=spatialvb" >> $current_options_file
    echo "--noise=white" >> $current_options_file
    echo "--model=turboquasar" >> $current_options_file
    echo "--disp=$disp" >> $current_options_file
    echo "--inferart" >> $current_options_file
    echo "--artdir" >> $current_options_file
    echo "--tissoff" >> $current_options_file
    echo "--repeats=1" >> $current_options_file
    echo $tislist >> $current_options_file
    echo "--t1=$t1 --t1b=$t1b --tau=$tau --fa=$fa " >> $current_options_file
    #if [ ! -z $infert1 ]; then
    #    #echo "--t1_MT=$t1 --t1b_MT=$t1b " >> $current_options_file
    #fi
    echo "--slicedt=$slicedt" >> $current_options_file
    echo $bolus_pattern_list >> $current_options_file
    echo "--slice_shift=$shift_factor" >> $current_options_file
    echo "--onephase" >> $current_options_file
    if [ ! -z $infert1 ]; then    
        echo "--usecalib " >> $current_options_file # We incorporate previously estimated T1 and g value as priors
    fi
    #echo "--infert1" >> $current_options_file # If we infer T1 we must infer T1 of tissue and blood together (must be two of them)
    #echo "--artdir" >> $current_options_file
    # Save model fitting results and residue
    echo "--save-model-fit" >> $current_options_file
    echo "--print-free-energy" >> $current_options_file
    echo "--save-residuals" >> $current_options_file
    # Here we create a shortcut to the latest results directory
    #echo "--link-to-latest" >> $current_options_file

    # Here are the settings
    echo "--PSP_byname1=fblood" >> $current_options_file # Estimate CBF
    echo "--PSP_byname1_type=M" >> $current_options_file
    echo "--PSP_byname2=deltblood" >> $current_options_file # Estimate ATT
    echo "--PSP_byname2_type=N" >> $current_options_file
    echo "--PSP_byname3=sp_log" >> $current_options_file # Estimate sp log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname3_type=N" >> $current_options_file
    echo "--PSP_byname4=s_log" >> $current_options_file # Estimate s log (dispersion parameter, even if dispersion is none)
    echo "--PSP_byname4_type=N" >> $current_options_file
    echo "--PSP_byname5=thblood" >> $current_options_file # Estimate angle theta of crushing gradient
    echo "--PSP_byname5_type=N" >> $current_options_file
    echo "--PSP_byname6=phiblood" >> $current_options_file # Estimate angle phi of crushing gradient
    echo "--PSP_byname6_type=N" >> $current_options_file
    echo "--PSP_byname7=bvblood" >> $current_options_file # Estimate volume of blood
    echo "--PSP_byname7_type=N" >> $current_options_file
    if [ ! -z $infert1 ]; then
        echo "--PSP_byname8=g" >> $current_options_file # Incorporate FA correction factor g
        echo "--PSP_byname8_type=I" >> $current_options_file
        echo "--PSP_byname8_image=output_calib_shift_2/mean_g" >> $current_options_file
        #echo "--PSP_byname9=T_1b_MT" >> $current_options_file # T1 blood under MT effects
        #echo "--PSP_byname9_type=I" >> $current_options_file
        #echo "--PSP_byname9_image=output_calib_shift_2_MT/mean_T1b" >> $current_options_file
    fi

    # Now perform model fitting to estimate ABV and arrival time to macrovasculature
    $fabber --data=$shift_2_repeat_1_blood --data-order=singlefile --output=output_blood_shift_2 -@ $current_options_file

fi



# Combine the two shifts (take the average)
# Create an output directory
output_dir_temp="output_dir_temp"
mkdir $output_dir_temp

# Rearrage shift 2
# M0t
fslroi output_calib_shift_2/mean_M0t data_segment_1 0 64 0 64 $break_1_result $segment_1_length_result
fslroi output_calib_shift_2/mean_M0t data_segment_2 0 64 0 64 $break_2_result $segment_2_length_result
fslmerge -z output_calib_shift_2/mean_M0t_result data_segment_2 data_segment_1 # Segment 2 was scanned first
# T1t
fslroi output_calib_shift_2/mean_T1t data_segment_1 0 64 0 64 $break_1_result $segment_1_length_result
fslroi output_calib_shift_2/mean_T1t data_segment_2 0 64 0 64 $break_2_result $segment_2_length_result
fslmerge -z output_calib_shift_2/mean_T1t_result data_segment_2 data_segment_1 # Segment 2 was scanned first
# CBF_relative
fslroi output_tissue_shift_2/mean_ftiss data_segment_1 0 64 0 64 $break_1_result $segment_1_length_result
fslroi output_tissue_shift_2/mean_ftiss data_segment_2 0 64 0 64 $break_2_result $segment_2_length_result
fslmerge -z output_tissue_shift_2/mean_ftiss_result data_segment_2 data_segment_1 # Segment 2 was scanned first

# ATT
fslroi output_tissue_shift_2/mean_delttiss data_segment_1 0 64 0 64 $break_1_result $segment_1_length_result
fslroi output_tissue_shift_2/mean_delttiss data_segment_2 0 64 0 64 $break_2_result $segment_2_length_result
fslmerge -z output_tissue_shift_2/mean_delttiss_result data_segment_2 data_segment_1 # Segment 2 was scanned first
# Bolus duration
if [ ! -z $infertau ]; then
    fslroi output_tissue_shift_2/mean_tautiss data_segment_1 0 64 0 64 $break_1_result $segment_1_length_result
    fslroi output_tissue_shift_2/mean_tautiss data_segment_2 0 64 0 64 $break_2_result $segment_2_length_result
    fslmerge -z output_tissue_shift_2/mean_tautiss_result data_segment_2 data_segment_1 # Segment 2 was scanned first
fi
# T1 and M0t of MT effects
if [ ! -z $infert1 ]; then
    # T1
    fslroi output_calib_shift_2_MT/mean_T1t data_segment_1 0 64 0 64 $break_1_result $segment_1_length_result
    fslroi output_calib_shift_2_MT/mean_T1t data_segment_2 0 64 0 64 $break_2_result $segment_2_length_result
    fslmerge -z output_calib_shift_2_MT/mean_T1t_result data_segment_2 data_segment_1 # Segment 2 was scanned first
    # M0t
    fslroi output_calib_shift_2_MT/mean_M0t data_segment_1 0 64 0 64 $break_1_result $segment_1_length_result
    fslroi output_calib_shift_2_MT/mean_M0t data_segment_2 0 64 0 64 $break_2_result $segment_2_length_result
    fslmerge -z output_calib_shift_2_MT/mean_M0t_result data_segment_2 data_segment_1 # Segment 2 was scanned first
    # CBF relative MT
    fslroi output_tissue_shift_2_MT/mean_ftiss data_segment_1 0 64 0 64 $break_1_result $segment_1_length_result
    fslroi output_tissue_shift_2_MT/mean_ftiss data_segment_2 0 64 0 64 $break_2_result $segment_2_length_result
    fslmerge -z output_tissue_shift_2_MT/mean_ftiss_result data_segment_2 data_segment_1 # Segment 2 was scanned first  
    # T1b
    #fslroi output_calib_shift_2_MT/mean_T1b data_segment_1 0 64 0 64 $break_1_result $segment_1_length_result
    #fslroi output_calib_shift_2_MT/mean_T1b data_segment_2 0 64 0 64 $break_2_result $segment_2_length_result
    #fslmerge -z output_calib_shift_2_MT/mean_T1b_result data_segment_2 data_segment_1 # Segment 2 was scanned first   
fi
# Arterial blood colume and arrival time to macrovasculature
if [ ! -z $infer_arterial ]; then
    # ABV
    fslroi output_blood_shift_2/mean_fblood data_segment_1 0 64 0 64 $break_1_result $segment_1_length_result
    fslroi output_blood_shift_2/mean_fblood data_segment_2 0 64 0 64 $break_2_result $segment_2_length_result
    fslmerge -z output_blood_shift_2/mean_fblood_result data_segment_2 data_segment_1 # Segment 2 was scanned first
    # Arrival time to macrovasculature
    fslroi output_blood_shift_2/mean_deltblood data_segment_1 0 64 0 64 $break_1_result $segment_1_length_result
    fslroi output_blood_shift_2/mean_deltblood data_segment_2 0 64 0 64 $break_2_result $segment_2_length_result
    fslmerge -z output_blood_shift_2/mean_deltblood_result data_segment_2 data_segment_1 # Segment 2 was scanned first    
fi

# Now compute average CBF, ATT, M0t, T1t, bolus duration
fslmaths output_tissue_shift_1/mean_ftiss -add output_tissue_shift_2/mean_ftiss_result -div 2 -thr 0 $output_dir_temp/CBF_relative
# If corrected MT effects, we need to compute a relative CBF under MT effect
if [ ! -z $infert1 ]; then
    fslmaths output_tissue_shift_1_MT/mean_ftiss -add output_tissue_shift_2_MT/mean_ftiss_result -div 2 -thr 0 $output_dir_temp/CBF_relative_MT
fi
fslmaths output_tissue_shift_1/mean_delttiss -add output_tissue_shift_2/mean_delttiss_result -div 2 -thr 0 $output_dir_temp/ATT
fslmaths output_calib_shift_1/mean_T1t -add output_calib_shift_2/mean_T1t_result -div 2 -thr 0 $output_dir_temp/T1_tissue
# Here we need to handle MT effects M0
if [ ! -z $infert1 ]; then
    fslmaths output_calib_shift_1_MT/mean_M0t -add output_calib_shift_2_MT/mean_M0t_result -div 2 -thr 0 $output_dir_temp/M0t_from_control_data
# If we dont correct MT effects, then just use normal M0
else
    fslmaths output_calib_shift_1/mean_M0t -add output_calib_shift_2/mean_M0t_result -div 2 -thr 0 $output_dir_temp/M0t_from_control_data
fi
if [ ! -z $infertau ]; then
    fslmaths output_tissue_shift_1/mean_tautiss -add output_tissue_shift_2/mean_tautiss_result -div 2 -thr 0 $output_dir_temp/Bolus_duration_raw
    # In the analysis the estimated (observed) bolus duration is actually the results of an operation of the actual bolus duration
    # ie. actual_tau_to_be_estimated = (dti - tau_lowest) / 2 * tanh(mean_tautiss) + (dti + tau_lowest) / 2
    # So we need to apply function of tanh to obtain the actual bolus duration to be estimated
    # tanh(x) = (exp(2x)-1)/(exp(2x)+1)
    # Formula: (dti - tau_lowest) / 2 * tanh(x) + (dti + tau_lowest) / 2
    # slope = (dti - tau_lowest) / 2
    # tanh(x) = (exp(2x)-1)/(exp(2x)+1)
    # intercept = (dti + tau_lowest) / 2
    # This needs to be updated
    # dti=0.6 (not always the case)
    dti=0.6
    slope_sub=`echo print $dti-$tau_lowest | perl`
    slope=`echo print $slope_sub/2 | perl`
    intercept_sub=`echo print $dti+$tau_lowest | perl`
    intercept=`echo print $intercept_sub/2 | perl`
    fslmaths output_dir_temp/Bolus_duration_raw -mul 2 -exp -sub 1 output_dir_temp/numerator
    fslmaths output_dir_temp/Bolus_duration_raw -mul 2 -exp -add 1 output_dir_temp/denominator
    fslmaths output_dir_temp/numerator -div output_dir_temp/denominator output_dir_temp/tanh
    fslmaths output_dir_temp/tanh -mul $slope -add $intercept -mas mask output_dir_temp/Bolus_duration

    # We now use the estimated bolus duration as the final result
    fslmaths output_tissue_shift_1/mean_tautiss -add output_tissue_shift_2/mean_tautiss_result -div 2 -thr 0 $output_dir_temp/Bolus_duration
    #imrm Bolus_duration_raw numerator denominator tanh
fi
# T1 and M0t of MT effects
if [ ! -z $infert1 ]; then
    # T1
    fslmaths output_calib_shift_1_MT/mean_T1t -add output_calib_shift_2_MT/mean_T1t_result -div 2 -thr 0 $output_dir_temp/T1_tissue_MT
    # M0t
    fslmaths output_calib_shift_1_MT/mean_M0t -add output_calib_shift_2_MT/mean_M0t_result -div 2 -thr 0 $output_dir_temp/M0t_MT
    # T1b under MT effects
    #fslmaths output_calib_shift_1_MT/mean_T1b -add output_calib_shift_2_MT/mean_T1b_result -div 2 -thr 0 $output_dir_temp/T1_blood_MT
    # T1b normal
    #fslmaths $output_dir_temp/T1_blood_MT -mul 0 -add $t1b -mas $output_dir_temp/T1_blood_MT $output_dir_temp/T1_blood
fi
# Arterial blood colume and arrival time to macrovasculature infert1
if [ ! -z $infer_arterial ]; then
    # ABV
    fslmaths output_blood_shift_1/mean_fblood -add output_blood_shift_2/mean_fblood_result -div 2 -thr 0 $output_dir_temp/ABV_relative
    # arrival time to macrovasculature
    fslmaths output_blood_shift_1/mean_deltblood -add output_blood_shift_2/mean_deltblood_result -div 2 -thr 0 $output_dir_temp/Arrival_time_blood
fi



# Additional calibration data is provided
if [ ! -z $calib_data ]; then
    echo "Using user provided calibration data..."
    # Move the calib data into the working directory
    imcp $stdir"/"$calib_data ./

    # We need to perform transformation on the calibration data
    if [ ! -z $trans_status ]; then
        echo "Transforming calibration data to Turbo QUASAR ASL resolution..."
        # Move the structural data into the working directory
        imcp $stdir"/"$data_struct ./
        # Extract brain from the structure and calibration data
        struct_brain="struct_brain"
        bet $data_struct $struct_brain
        calib_brain="calib_brain"
        bet $calib_data $calib_brain
        #echo b
        # Register ASL to Structure
        echo "Running Turbo QUASAR ASL to structure registration..." 
        flirt -ref $struct_brain -in $output_dir_temp/CBF_relative -out CBF_high -omat turbo_to_struct.mat
        #flirt -ref $struct_brain -in $output_dir_temp/T1_tissue -out CBF_high -omat turbo_to_struct.mat
        echo "Running calibration image to structure registration..."
        flirt -ref $struct_brain -in $calib_brain -out calib_high -out calib_high -omat calib_to_struct.mat
        echo "Matrix inversion and concatenation..."
        convert_xfm -omat struct_to_turbo.mat -inverse turbo_to_struct.mat
        convert_xfm -omat struct_to_calib.mat -inverse calib_to_struct.mat
        convert_xfm -omat calib_to_turbo.mat -concat struct_to_turbo.mat calib_to_struct.mat
        convert_xfm -omat turbo_to_calib.mat -concat struct_to_calib.mat turbo_to_struct.mat
        echo "Transforming calibration image to Turbo QUASAR ASL resolution..."
        applywarp --ref=$output_dir_temp/CBF_relative --in=$calib_brain --out=$output_dir_temp/calib_M0t --premat=calib_to_turbo.mat --super --interp=spline --superlevel=4
        #applywarp --ref=$output_dir_temp/T1_tissue --in=$calib_brain --out=$output_dir_temp/calib_M0t --premat=calib_to_turbo.mat --super --interp=spline --superlevel=4



        # Copy the trasnformation matrices to result dirctory
        cp turbo_to_struct.mat calib_to_struct.mat struct_to_turbo.mat struct_to_calib.mat calib_to_turbo.mat turbo_to_calib.mat $output_dir_temp

    # We do not need to perform transformation
    # Copy the calibration image to temp output folder
    else
        imcp $stdir"/"$calib_data $output_dir_temp"/"calib_M0t
    fi

    # Correct TR if the TR of the calibration image is less than 5s
    # We only do this if the calibration image is provided by user
    TR_calib_in_ms=`echo print $TR_calib*1000 | perl`
    # If else in bash can only handle integer operations
    if [ $TR_calib_in_ms -lt 5000 ]; then
        echo "Correcting TR of calibration image..."
        #statements
        # 10.1002/mrm.25197
        # multiple_factor=(1/(1-exp(-TR/T1_of_tissue)))
        correction_factor=`echo "(1/(1-e((-1)*$TR_calib/$t1)))" | bc -l`
        fslmaths $output_dir_temp/calib_M0t -mul $correction_factor $output_dir_temp/calib_M0t
    fi

# Use the estimated M0t image from the previous step (saturation recovery model)
else
    echo "Using calibration data estimated from fitting ASL control data to the saturation recovery signal..."
    imcp $output_dir_temp"/"M0t_from_control_data $output_dir_temp/calib_M0t


fi

# Correct the calibration image by median filter, erosion, and extrapolation
# Ref: http://www.sciencedirect.com/science/article/pii/S1053811917307103
fslmaths $output_dir_temp/calib_M0t -mas mask $output_dir_temp/calib_M0t
if [ ! -z $corrcal ]; then
    echo "Correcting partial volume effects on the edge of the calibration data..."
    fslmaths $output_dir_temp/calib_M0t -fmedian M0t_median
    fslmaths M0t_median -ero M0t_ero
    asl_file --data=M0t_ero --ntis=1 --mask=mask --extrapolate --neighbour=5 --out=M0t_extra
    fslmaths M0t_extra -div $lambda $output_dir_temp/M0a_for_absolute_CBF
    #imrm M0t_median M0t_ero M0t_extra

    # If we corrected MT effects
    if [ ! -z $infert1 ]; then
        fslmaths $output_dir_temp/M0t_MT -fmedian M0t_median
        fslmaths M0t_median -ero M0t_ero
        asl_file --data=M0t_ero --ntis=1 --mask=mask --extrapolate --neighbour=5 --out=M0t_extra
        fslmaths M0t_extra -div $lambda $output_dir_temp/M0a_for_absolute_CBF_MT
    fi

# Otherwise, we dont correct the edge problems
else
    fslmaths $output_dir_temp/calib_M0t -div $lambda $output_dir_temp/M0a_for_absolute_CBF
fi



# Now compuate CBF in absolute unit
fslmaths $output_dir_temp/CBF_relative -div $output_dir_temp/M0a_for_absolute_CBF -div $alpha -mul 6000 $output_dir_temp/CBF_absolute
# If we corrected MT effects
if [ ! -z $infert1 ]; then
    fslmaths $output_dir_temp/CBF_relative_MT -div $output_dir_temp/M0a_for_absolute_CBF_MT -div $alpha -mul 6000 $output_dir_temp/CBF_absolute_MT
    # Now we take the average of the two CBFs
    fslmaths $output_dir_temp/CBF_absolute -add $output_dir_temp/CBF_absolute_MT -div 2 $output_dir_temp/CBF_absolute
fi
# Arterial blood colume in absolute unit
if [ ! -z $infer_arterial ]; then
    echo "Output ABV is in decimal format"
    fslmaths $output_dir_temp/ABV_relative -div $output_dir_temp/M0a_for_absolute_CBF -div $alpha $output_dir_temp/ABV_absolute
fi

# Copy results to output directory
#cp -R $output_dir_temp $outdir
cp -a $output_dir_temp/. $outdir
#imcp $tempdir/mask $outdir/mask

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

echo "Turbo QUSAR ASL analysis done!"

