asl_calib.in

#!/bin/sh

# ASL_CALIB: Calibration for ASL data
#
# Michael Chappell & Brad MacIntosh, FMRIB Image Analysis & Physics Groups
#
# Copyright (c) 2008-2013 Univerisity of Oxford
#
# SHCOPYRIGHT

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

#deal with options

Usage() {
    echo "ASL_CALIB"
#    echo "Version: 1.2b"
    echo "Calibration for ASL data"
    echo ""
    echo "Usage (optional parameters in {}):" 
    echo " -c         : specify calibration image (stacked form)"
    echo " -s         : specify structural image (already BETed)"
    echo " -t         : specify asl-->structural transformation matrix"
    echo " {--mode}   : Calibration mode: longtr or satrecov (see below) {default: longtr}"
    echo " {--tissref}: Tissue reference type: csf, wm, gm or none {default: csf}"
    echo " {--te}     : TE used in sequence (ms) - {default: 0 ms (i.e. negligible)}"
    echo " {-i}       : specify a CBF image for calibration (should be in ASL native space)"
    echo ""
    echo " Output options (set at least one of):"
    echo " > save all results to a given directory:"
    echo " {-o}       : specify output directory name."
    echo " > save specific results to indivdually anmed files:"
    echo " {--of}     : specify output filename for calibrated image - {default: <input_file>_calib}"
    echo "              requires -i option to have been set"
    echo " {--Mo}     : save the calculated M0 value to a specified file"  
    echo " {--om}     : save CSF mask to a specified file"
    echo ""
    echo " Extended options (all optional):"
    echo " -m         : Specify reference mask in calibration image space"
    echo "              - strucutral image & transformation matrix are not required"
    echo " --bmask    : Brain mask (in ASL data space) for sensitivity or tissue T1 estimation"
    echo " --t2star   : Correct with T2* rather than T2"
    echo "               (this alters the default values specified below to the T2* values)"
    echo " --t1r      : T1 of reference tissue (defaults: csf 4.3, gm 1.3, wm 1.0 s) "
    echo " --t2r      : T2(*) of reference tissue (defaults T2/T2*: csf 750/400, gm 100/60,  wm 50/50  ms)"
    echo " --t2b      : T2(*) of blood (default T2/T2*: 150/50 ms)"
    echo " --pc       : Partition co-efficient (defaults csf 1.15, gm 0.98,  wm 0.82)"
    echo " --alpha    : Specify inversion efficiency - only applied to final CBF image calculation"
    echo ""
    echo " CSF masking options (only for --tissref csf)"
    echo "  By default asl_calib extracts CSF from the structural image by segmentation and"
    echo "  this is then masked using the ventricles in MNI152 space."
    echo " --csfmaskingoff : turns off the ventricle masking, reference is based on segmentation only."
    echo "  Registration between structural image and MNI152 is done automatically unless:"
    echo "  --str2std  : Structural to MNI152 linear registration (.mat)"
    echo "  --warp     : Structural to MNI152 non-linear registration (warp)"
    echo ""
    echo "MODES:"
    echo "> longtr  Calibration image is a control image with a long TR."
    echo "  {--tr}     : TR used in calibration sequence - {default: 3.2s}"
    echo "  {--cgain}  : Relative gain between calibration and ASL data - {default: 1}"
    echo ""
    echo "> satrecov  Calibration image is a sequnce of control images at various TIs"
    echo "            M0 is to be determined from a saturation recovery"
    echo "            T1 of tissue (and FA correction) images are also calcualted"
    echo " --tis       : comma separated list of inversion times, e.g. --tis 0.2,0.4,0.6"
    echo " {--fa}      : Flip angle (in degrees) for Look-Locker readouts"
    echo "   >> Look-Locker flip angle correction - to perform this provide:"
    echo " {--lfa}     : Lower flip angle (in degrees) for dual FA calibration"
    echo " {--nphases} : Number of phases (repetitions) of higher FA"
    echo " {--fixa}    : Fix the saturation efficiency to 100% (useful if you have a low number of samples)"
    echo ""
    echo "Coil sensitivity correction:"
    echo " Calculate and apply a voxel-wise correction for coil sensitivity"
    echo " > using bias field from structural image (default)"
    echo " {--osen}    : save sensitivity image to specified file."
    echo " > using existing sensitivity image:"
    echo "  --isen     : input coil sensitivity image"
    echo " > using reference images (collected using same parameters):"
    echo "  --cref     : Reference image from coil with minimal variation e.g. body."
    echo "  {--cact}   : Image from coil used for actual ASL acquisition"
    echo "               {default: calibration image - only in longtr mode}"
    echo ""
}

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

if [ -z $1 ]; then
    Usage
    exit 1
fi

until [ -z $1 ]; do
    case $1 in
	--mode) mode=$2
	    shift;;
	-o) outdir=$2
	    shift;;
	--of) outflag=1  outfile=$2
	    shift;;
	-i) inflag=1 infile=$2
	    shift;;
	-c) calibflag=1 calib=$2
	    shift;;
	-s) strucflag=1 struc=$2
	    shift;;
	-t) transflag=1 trans=$2
	    shift;;
	--tr) trflag=1 tr=$2
	    shift;;
	--taq) taqflag=1 taq=$2
	    shift;;
	--te) teflag=1 te=$2
	    shift;;
	--tissref) tissref=$2
	    shift;;
	--t1r) T1rin=$2
	    shift;;
	--t2r) T2rin=$2
	    shift;;
	--t2b) T2bin=$2
	    shift;;
	--t2star) t2star=1
	    ;;
	--pc) pcin=$2
	      shift;;
	--alpha) alpha=$2
		 shift;;
	--tis) tis=$2
	    shift;;
	-m) maskflag=1 mask=$2
	    shift;;
	--bmask) bmask=$2
	    shift;;
	--refpve) refpve=$2 # this is the PVE image for the reference tissue (in strucural image space)
	    shift;;
	--of) offlag=1 outfact=$2
	    shift;;
	--Mo) Moflag=1 outMo=$2
	    shift;;
	--om) omflag=1 outmask=$2
	    shift;;
	--str2std) str2std=$2;
	    shift;;
	--warp) warp=$2
	    shift;;
	--csfmaskingoff) csfmaskingoff=1 # turn off the masking of CSF using MNI152 ventricles
	    ;;
	--cgain) cgain=$2
	shift;;
	--cref) crefflag=1; crefim=$2
	    shift;;
	--cact) cact=$2
	    shift;;
	--osen) osenflag=1; senout=$2
	    shift;;
	--isen) seninflag=1 senin=$2
	    shift;;
	--fa) fa=$2
	    shift;;
	--lfa) lfa=$2
	    shift;;
	--nphases) nphases=$2
	    shift;;
	--fixa) fixa=1
	    ;;
	--devel) devel=1
	    ;;
	--debug) debug=1
	    ;;
	--version) Version
	    ;;
	*)  Usage
	    echo "Error! Unrecognised option on command line: $1"
	    echo ""
	    exit 1;;
    esac
    shift
done

echo "ASL_CALIB"

# set the version of fabber to use
fabber=fabber_asl

#check for mandatory inputs
if [ -z $calib ]; then
    echo "ERROR: calibration image has not been specified"
    Usage
fi
if [ -z $mask ]; then
# if a reference mask has not bee supplied in $mask, then there should be a structural with transformation matrix from calibration space
    if [ -z $struc ]; then
	echo "ERROR: Structural image has not been supplied (alternatively supply a reference mask with -m)"
	Usage
    fi
    if [ -z $trans ]; then
	echo "ERROR: Transformation matrix has not been supplied (alternatively supply a reference mask with -m)"
	Usage
    fi
fi

# check if we have an input perfusion map
if [ ! -z $inflag ]; then
    echo "Input file is: $infile"
   #strip off extension from input file
   infile=`imglob $infile`

# set the output filename here if not specified
   if [ -z $outflag ]; then
       outfile=${infile}_calib;
   fi
fi

# create output directory if required
if [ ! -d $outdir ]; then
    mkdir $outdir
fi

# make a temporary directory to work in - delete at end
tmpbase=`tmpnam`
temp_calib=${tmpbase}_asl_calib
mkdir $temp_calib

# start the asl_calib logfile
log=$temp_calib/logfile
echo "ASL_CALIB logfile" > $log
echo $# >> $log

# set defaults
if [ -z $mode ]; then
    mode=longtr;
fi

if [ -z $trflag ]; then
    tr=3.2;
fi

if [ -z $teflag ]; then
    te=0;
fi
echo "TE: $te" >> $log

if [ -z $taq ]; then
    taq=0;
fi

# if we have a sensitvity map at input copy to right place
if [ ! -z $senin ]; then
    echo "Using sensitivity image: $senin" >> $log
    senson=1
    imcp $senin $temp_calib/sens
fi

# constants
T1csf=4.3
T2csf=750
T1gm=1.3
T2gm=100
T1wm=1.0
T2wm=50

T2b=150 # lu et a. 2012 MRM 67:42-49 have 154ms at 3T during normoxia

# partition coeffs
# based on Herscovitch and Raichle 1985
pccsf=1.15 # a blood water density of 0.87
pcwm=0.82
pcgm=0.98

if [ ! -z $t2star ]; then
# we need to correct for T2* not T2 so change the defaults
# NB these will still be overridden by specific values supplied
    T2csf=400
    T2gm=60 # from Foucher 2011 JMRI 34:785-790
    T2wm=50 # ditto

    T2b=50 #from Petersen 2006 MRM 55(2):219-232 see discussion
fi

# deal with tissue reference choice
if [ -z $tissref ]; then
    # set the default
	tissref=csf
	echo "Using default tissue reference type (CSF)" >> $log
fi

case $tissref in
	csf) T1r=$T1csf; T2r=$T2csf; fastpve=0; pc=$pccsf;;
	wm)  T1r=$T1wm;  T2r=$T2wm;  fastpve=2; pc=$pcwm;;
	gm)  T1r=$T1gm;  T2r=$T2gm;  fastpve=1; pc=$pcgm;;
	0)   fastpve=0;;
	1)   fastpve=1;;
	2)   fastpve=2;;
	none) ;;
	*)  echo "Error! Unrecognised tissue reference type: $tissref"
	    exit 1;;
    esac

echo "Tissue reference is: $tissref"
echo "Tissue reference : $tissref" >> $log

# command line override of default T1 and T2
if [ ! -z $T1rin ]; then
    if [ `${FSLDIR}/bin/imtest $T1rin` -gt 0 ]; then
	# we have been supplied with a T1 image - record that here to process later
	T1rim=$T1rin
	echo "Loading T1 image for reference region: T1im" >> $log
    else
	T1r=$T1rin
	echo "Setting T1r from command line: $T1r" >> $log
    fi
fi
if [ ! -z $T2rin ]; then
    if [ `${FSLDIR}/bin/imtest $T2rin` -gt 0 ]; then
	# we have been supplied with a T2 image - record that here to process later
	T2rim=$T2rin
	echo "Loading T2 image for reference region: T2im" >> $log
    else
	T2r=$T2rin
	echo "Setting T2r from command line: $T2r" >> $log
    fi
fi
if [ ! -z $T2bin ]; then
    T2b=$T2bin
    echo "Setting T2b from command line: $T2b" >> $log
fi
if [ ! -z $pcin ]; then
    pc=$pcin
    echo "Setting partition coefficient from command line: $pcin" >> $log
fi

# make sure the T1 has been set (either by user or above)
if [ -z $T1r ]; then
    echo "Error! T1 for reference tissue has not been set."
    exit 1
fi
# make sure the T2 has been set (either by user or above)
if [ -z $T2r ]; then
    # we only need T2 to have a meaningful value if we care about the TE
    if [ ! -z $teflag ]; then
	echo "Error! T2 for reference tissue has not been set."
	exit 1
    fi
    T2r=1.0 #a default value when TE=0
fi
if [ -z $pc ]; then
    echo "Error! Partition coefficient for reference tissue has not been set."
    exit 1
fi

echo "T1r $T1r; T2r $T2r; T2b $T2b; Part co-eff: $pc" >> $log

# sort out the M0 calib brain_mask
if [ -z $bmask ]; then
    echo "Creating brain mask from calibration image" >> $log
    #make a brain mask
    # take the mean
    fslmaths $calib -Tmean $temp_calib/calib_mean
    # bet
    bet $temp_calib/calib_mean $temp_calib/calib_mean -m #calib_mean_mask is the brain mask for the calib image
    bmask=$temp_calib/calib_mean_mask
fi

### Sensitivity image calculation (if reqd)
if [ ! -z $crefim ]; then
    echo "Calculate sensitivity image" >> $log
    senson=1
    # take the mean (and mask with the mask from the main calib image)
    fslmaths $crefim -Tmean -mas $bmask $temp_calib/crefim
    
    # take the ratio to give the sensitivity image
    if [ -z $cactim ]; then
        # if the cact image has not been supplied then use the mean of the calib image
	if [ ! $mode = longtr ]; then
	    echo "ERROR: You must supply an image from the actual coil used for ASL acquisition using --cact (unless you use longtr mode)"
	    exit 1
	fi
	fslmaths $calib -Tmean $temp_calib/cactim
    fi
    
    fslmaths $temp_calib/cactim -div $temp_calib/crefim -mas $bmask $temp_calib/sens
    
fi

if [ $tissref = "none" ]; then
# whole brain M0
# in this case use the brain mask
    imcp $bmask $temp_calib/refmask
    maskflag=1
    echo "Brain mask is being used as the reference tissue (beware!)" >> $log
fi

if [ -z $maskflag ]; then

# make brain mask from structural
    fslmaths $struc -bin $temp_calib/mask

    if [ -z $refpve ]; then
        # auto create tissue reference mask
	echo "FAST called to determine a reference tissue mask" >> $log
	
        # segment structural image
	fast -b -o $temp_calib/seg -p $struc
	fasthasrun=1;
	imcp  $temp_calib/seg_pve_$fastpve $temp_calib/refpve

    else
	# user supplied PV estimate for reference tissue
	echo "Using input reference PVE: $refpve" >> $log
	imcp $refpve $temp_calib/refpve
    fi

    if [ $tissref = "csf" ] && [ -z $csfmaskingoff ]; then
	echo "Ventricle selection" >> $log
	stdmaskfnirt=1  # by deafult now we do FNRIT transformation of ventricle mask

# cut down brain mask so that it only covers middle of brain
# sort out the roi
#	xsize=`fslinfo $struc | grep "^dim1" | sed 's:dim1[ ]*::'`
#	ysize=`fslinfo $struc | grep "^dim2" | sed 's:dim2[ ]*::'`
#	zsize=`fslinfo $struc | grep "^dim3" | sed 's:dim3[ ]*::'`
#	roisize="0.3";
#echo "$xsize $ysize $zsize $roisize"
#	delx=`echo "v = $roisize * $xsize; v /= 1; v" | bc`
#	xmin=`echo "v = 0.5 * $xsize - $delx / 2; v /= 1; v" | bc`
#	dely=`echo "v = $roisize * $ysize; v /= 1; v" | bc`
#	ymin=`echo "v = 0.5 * $ysize - $dely / 2; v /= 1; v" | bc`
#	delz=`echo "v = $roisize * $zsize; v /= 1; v" | bc`
#	zmin=`echo "v = 0.5 * $zsize - $delz / 2 + 0.1; v /= 1; v" | bc`
#	echo "$xmin $delx $ymin $dely $zmin $delz" >> $log
#	fslmaths $temp_calib/mask -roi $xmin $delx $ymin $dely $zmin $delz 0 1 $temp_calib/mask

        # select ventricles based on standard space atlas
	fslroi $FSLDIR/data/atlases/HarvardOxford/HarvardOxford-sub-prob-2mm $temp_calib/LVentricle 2 1
	fslroi $FSLDIR/data/atlases/HarvardOxford/HarvardOxford-sub-prob-2mm $temp_calib/RVentricle 13 1
	fslmaths $temp_calib/LVentricle -add $temp_calib/RVentricle -thr 0.1 -bin -ero $temp_calib/VentricleMask
	
        # register structural image to std space using FLIRT
	if [ -z $str2std ] && [ -z $warp ]; then
	    echo "Registering structural image to standard space using FLIRT" >> $log
	    flirt -in $struc -ref $FSLDIR/data/standard/MNI152_T1_2mm_brain.nii.gz -omat $temp_calib/struc2std.mat
	else
	    if [ -f $str2std ]; then
		cp $str2std $temp_calib/struc2std.mat
	    else
		echo "Error $str2std not found"
		exit 1
	    fi
	fi
	
	if [ -z $stdmaskfnirt ] && [ -z $warp ]; then
        # do std space masking with FLIRT registration
	    
	    convert_xfm -omat $temp_calib/std2struc.mat -inverse $temp_calib/struc2std.mat
	    flirt -in $temp_calib/VentricleMask -applyxfm -init $temp_calib/std2struc.mat -ref $temp_calib/refpve -out $temp_calib/mask
	else
            # do std space masking using FNIRT registration
            # Register the structural image to standard space using FNIRT
	    if [ -z $warp ]; then
		echo "Registering structural image to standard space using FNIRT" >> $log
		fnirt --in=$struc --aff=$temp_calib/struc2std.mat --config=T1_2_MNI152_2mm.cnf --cout=$temp_calib/coef_struc2MNI152
	    else
		if [ -f $warp ]; then
		    cp $warp $temp_calib/coef_struc2MNI152
		else
		    echo "Error: $warp not found"
		    exit 1
		fi
	    fi
            # Calculate the inverse warp using INVWARP and apply to standard space ventricle mask
	    invwarp --ref=$struc --warp=$temp_calib/coef_struc2MNI152 --out=$temp_calib/coef_MNI2struc
	    applywarp --ref=$temp_calib/refpve --in=$temp_calib/VentricleMask --warp=$temp_calib/coef_MNI2struc --out=$temp_calib/mask --interp=nn
	    # using refpve as reference for applywarp as this is the image we apply the mask to, nit the structural
	fi
	
	echo "Masking FAST output with standard space derrived ventricle mask" >> $log
	fslmaths $temp_calib/refpve -mas $temp_calib/mask $temp_calib/refpve
	if [ ! -z $debug ]; then
	    imcp $temp_calib/refmask $temp_calib/refmask_high_unthresh
	fi
    fi
    
    echo "Transforming tissue reference mask into perfusion space" >> $log
    #transform mask into perfusion space
    convert_xfm -omat $temp_calib/high2low.mat -inverse $trans
    #flirt -in $temp_calib/csf -applyxfm -init $temp_calib/high2low.mat -ref $calib -out $temp_calib/csf
    # new conversion using applywarp, supersmapling and integration
    applywarp --ref=$calib --in=$temp_calib/refpve --out=$temp_calib/refmask --premat=$temp_calib/high2low.mat --super --interp=spline --superlevel=4
    if [ ! -z $debug ]; then
	imcp $temp_calib/refmask $temp_calib/refmask_low_unthresh
    fi
    
    if [ ! -z $fasthasrun ] && [ -z $senson ]; then
        # also extract the bias field and convert to sensitivity image (as long as we have already been supplied by a sensivity iamge or reference)
	applywarp --ref=$calib --in=$temp_calib/seg_bias --out=$temp_calib/biasfield --premat=$temp_calib/high2low.mat --super --interp=spline --superlevel=4
	fslmaths $temp_calib/biasfield -recip $temp_calib/sens
	senson=1
	echo "Using bias field from structural image for sensitivity correction" >> $log
    fi

    
# threshold reference mask
    fslmaths $temp_calib/refmask -thr 0.9 -bin $temp_calib/refmask
    
# threshold reference mask and then keep top two clusters
#    cluster -i $temp_calib/refmask -t 0.9 -o $temp_calib/refmask_cluster --no_table
#    clusthr=`fslstats $temp_calib/refmask_cluster.nii.gz -R | awk '// {print $2}'`
#    clusthr=`echo "$clusthr - 1" | bc`
#    echo "Clustering threshold: $clusthr" >> $log
#    fslmaths $temp_calib/refmask_cluster -thr $clusthr $temp_calib/refmask
    
else
    #use supplied tissue reference mask
    imcp $mask $temp_calib/refmask
    echo "Using supplied reference tissue mask: $mask" >> $log
fi

#check there are some non-zero voxels
nzvox=`fslstats $temp_calib/refmask -V | awk '{print \$1}'`
echo "Number of voxels in tissue reference mask: $nzvox" >> $log
if [ "$nzvox" -lt 1 ]; then
    if [ -z $maskflag ]; then
	echo "ERROR: automatic masking has failed, check you have applied BET to the structural image, otherwise please provide a tissue reference mask."

	if [ ! -z $debug ]; then
	    cp -R $temp_calib ./temp_calib
	fi
	rm -r $temp_calib
	exit 1
    else
	echo "ERROR: no voxels left after transformation into perfusion space - check the supplied tissue reference mask."
	exit 1
    fi
    
fi
mask=$temp_calib/refmask

# ----

if [ $mode = longtr ]; then
    echo "MODE: longtr" >> $log
# Calibration data is a long TR acquisition - all we need to do here is take the mean 

# sort out cgain setting
    if [ -z $cgain ]; then
	cgain=1;        # default cgain is 1!
    fi
    echo "cgain is $cgain" >> $log

    #cut - throw away first volume
    #tsize=`fslinfo $calib | grep "^dim4" | sed 's:dim4[ ]*::'`
    #fslroi $calib $temp_calib/calib 1 $tsize
    # take the mean (again)
    fslmaths $calib -Tmean $temp_calib/calib

   if [ ! -z $senson ]; then
	echo "Apply sensitivity image" >> $log
        # apply sensitivity map to calibration image
	fslmaths $temp_calib/calib -div $temp_calib/sens $temp_calib/calib
    fi

    #mask M0 map with tissue reference
   fslmaths $temp_calib/calib -mas $mask $temp_calib/calib

   # calcualte T1 of reference region (if a T1 image has been supplied)
   if [ ! -z $T1rim ]; then
       fslmaths $T1rim -mas $mask $temp_calib/T1r
       T1r=`fslstats $temp_calib/T1r -M`
       echo "Calculated T1 of reference tissue: $T1r"
       echo "Calculated T1 of reference tissue: $T1r" >> $log
   fi

   # calcualte T2 of reference region (if a T2 image has been supplied)
   if [ ! -z $T2rim ]; then
       fslmaths $T2rim -mas $mask $temp_calib/T2r
       T2r=`fslstats $temp_calib/T2r -M`
       echo "Calculated T2 of reference tissue: $T2r"
       echo "Calculated T2 of reference tissue: $T2r" >> $log
   fi

    # calculate M0_ref value
    Moval=`fslstats $temp_calib/calib -M` # this is Mz of CSF
    echo "Signal in reference tissue: $Moval" >> $log
    corr=`echo "(1 - e(- ( $tr - $taq ) / $T1r) )" | bc -l`
    echo "Correction factor: $corr (tr=$tr, tsq=$taq, t1=$T1r)" >> $log
    Moval=`echo "$Moval / (1 - e(- ( $tr - $taq ) / $T1r) )" | bc -l` #this is now M0 of the reference
    echo "Mz of reference tissue: $Moval" >> $log

elif [ $mode = satrecov ]; then
    echo "MODE: satrecov" >> $log
    # Calibration image is control images and we want to do a saturation recovery fit
    # NB only do the fit in the CSF mask

    # sort out cgain setting
    if [ -z $cgain ]; then
	cgain=1;        # default cgain is 1!
    fi
    echo "cgain $cgain" >> $log
    
    #deal with TIs
    count=0
    tislist=""
    thetis=`echo $tis | sed 's:,: :g'`
    for ti in $thetis; do
	count=`expr ${count} + 1`
	tislist=`echo $tislist --ti${count}=$ti`
    done
    echo "TIs: tislist" >> $log

    # Extra options for Look Locker
    if [ ! -z $fa ]; then
	llopts="--FA=$fa"
	if [ ! -z $nphases ]; then
	    llopts=$llopts" --phases=$nphases"
	fi
	if [ ! -z $lfa ]; then
	    llopts=$llopts" --LFA=$lfa"
	fi
    fi
    echo "Look-Locker options: $llopts" >> $log

    # Extra sat recovery options
    if [ ! -z $fixa ]; then
	sropts="--fixa"
    fi
    echo "Extra saturation recovery options: $sropts" >> $log

    # do fabber within the tissue reference mask with a sensible T1 prior mean
    if [ ! -z $senson ]; then
	echo "Apply sensitivity image to data for reference tisse M0 estimation" >> $log
        # apply sensitivity map to calibration image - ONLY for the reference tissue calculations
	fslmaths $calib -div $temp_calib/sens $temp_calib/calib_senscorr
    else
	# no sensitivity correction required, but copy image over ready for next command
	imcp $calib $temp_calib/calib_senscorr
    fi
    echo "FABBER within reference tissue mask" >> $log
    $fabber --data=$temp_calib/calib_senscorr --mask=$mask --output=$temp_calib/satrecov --data-order=singlefile --model=satrecov --noise=white --method=vb $tislist $llopts $sropts --t1=$T1r

    # calculate M0 value
    Moval=`fslstats $temp_calib/satrecov/mean_M0t -M` # this is M0 of CSF at the TE of the sequence
    echo "M0 of reference tissue: $Moval" >> $log

    if [ ! -z $outdir ]; then
        # save useful results to specified output directory
	imcp $temp_calib/satrecov/mean_T1t $outdir/T1_ref
	imcp $temp_calib/satrecov/mean_M0t $outdir/M0_ref
    fi

    # do fabber again within whole brain to get estimated T1 of tissue and FA correction (if LL)
    # (note that we do not apply sensitivity correction to the data here - thius is 'built-into' the M0t map)
    echo "FABBER (again) within whole brain mask" >> $log
    if [ ! -z $outdir ]; then #NB we only bother with this if we have an output directory to put the results in
	$fabber --data=$calib --mask=$bmask --output=$temp_calib/satrecovT --data-order=singlefile --model=satrecov --noise=white --method=vb $tislist $llopts $sropts 


        # save useful results to specified output directory
	    imcp $temp_calib/satrecovT/mean_T1t $outdir/T1t
	    imcp $temp_calib/satrecovT/mean_M0t $outdir/M0t
	    if [ ! -z $lfa ]; then
		imcp $temp_calib/satrecovT/mean_g $outdir/facorr
	    fi
    fi


fi

# use equation to get the M0 value that is needed
#echo "$Moval $te $T2r $cgain $T2bl"
Moval=`echo "$Moval / e(- $te / $T2r )" | bc -l` # T2 correction for M0 reference
#echo $Moval
Moval=`echo "scale=2;$Moval*$cgain / $pc" | bc` #  this is M0 blood (at TE=0)
#echo $Moval
Moval=`echo "$Moval * e(- $te / $T2b )" | bc -l` # get M0 blood at TE used
echo "M0:$Moval"
echo "M0: $Moval" >> $log

# apply calibration to input image
if [ ! -z $inflag ]; then
    if [ ! -z $senson ]; then
	# apply sensitivity image
	fslmaths $infile -div $temp_calib/sens $temp_calib/infile
    else
	imcp $infile $temp_calib/infile
    fi

    if [ ! -z $alpha ]; then
	echo "Applying inversion efficiency of: $alpha" >> $log
	# apply the inversion efficiency supplied to M0 prior to final calculation
	Moval=`echo "scale=2; $Moval * $alpha" | bc`
    fi

    fslmaths $temp_calib/infile -mul 60 -mul 100 -div $Moval $outfile
fi

# save various things to the output directory (if specified)
if [ ! -z $outdir ]; then
    #save the M0 value in the output directory
    echo $Moval > $outdir/M0.txt

    # save the tissue reference mask
    imcp $temp_calib/refmask $outdir/refmask

    # copy the logfile across
    cp $temp_calib/logfile $outdir/logfile
fi

if [ ! -z $Moflag ]; then
#save the Mo value to the filenamed $outMo
#    echo "output M0: $Moval"
    echo $Moval > $outMo
fi

if [ ! -z $outmask ]; then
# save the tissue reference mask to a given file
    imcp $temp_calib/refmask $outmask
fi

if [ ! -z $offlag ]; then
#save the calibration factor to the filename $outfact
    factor=`echo "scale=2; 6000 / $Moval" | bc`
    echo $factor > $outfact
fi

if [ ! -z $senout ]; then
        # save sensitivity image
    if [ -f $temp_calib/sens.nii.gz ]; then
	imcp $temp_calib/sens $senout
	echo "Saving sensitivity image to: $senout" >> $log
    fi
fi

if [ ! -z $debug ]; then
    cp -R $temp_calib ./temp_calib
fi
rm -r $temp_calib

echo "ASL_calib - DONE."