#!/bin/sh # ASL_CALIB: Calibration for GRASE-ASL data # # Michael Chappell & Brad MacIntosh, FMRIB Image Analysis & Physics Groups # # Copyright (c) 2008-2011 Univerisity 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 5.0 (c) 2012, The University of # Oxford (the "Software") # # The Software remains the property of the University of Oxford ("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 Isis # Innovation Limited ("Isis"), the technology transfer company of the # University, to negotiate a licence. Contact details are: # innovation@isis.ox.ac.uk quoting reference DE/9564. #deal with options Usage() { echo "ASL_CALIB" echo "Version: 1.2b" echo "Calibration for GRASE-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: _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 "" 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 "" echo "Coil sensitivity correction (optional):" echo " Calculate and apply a voxel-wise correction for coil sensitivity" 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 " {--osen} : save sensitivity image to specified file." echo "" } 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 ;; --tis) tis=$2 shift;; -m) maskflag=1 mask=$2 shift;; --bmask) bmask=$2 shift;; --of) offlag=1 outfact=$2 shift;; --Mo) Moflag=1 outMo=$2 shift;; --om) omflag=1 outmask=$2 shift;; --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;; --stdmaskfnirt) stdmaskfnirt=1 ;; --devel) devel=1 ;; --debug) debug=1 ;; *) 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 #fabber=$FSLDIR/bin/fabber if [ ! -z $devel ]; then fabber=~chappell/cproject/fabber/fabber if [ $FSLMACHTYPE = linux_64-gcc4.1 ]; then #assumption that the architecuture tell me that we are on a jalapeno node echo "We are on jalapeno" fabber=~chappell/cproject/jalapeno/fabber/fabber fi fi #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 $calib ]; 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=`$FSLDIR/bin/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" 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 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;; wm) T1r=$T1wm; T2r=$T2wm; fastpve=2;; gm) T1r=$T1gm; T2r=$T2gm; fastpve=1;; 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 T1r=$T1rin echo "Setting T1r from command line: $T1r" >> $log fi if [ ! -z $T2rin ]; then T2r=$T2rin echo "Setting T2r from command line: $T2r" >> $log fi if [ ! -z $T2bin ]; then T2b=$T2bin echo "Setting T2b from command line: $T2b" >> $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 echo "T1r $T1r; T2r $T2r; T2b $T2b" >> $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 if [ ! -z $senout ]; then # save sensitivity image imcp $temp_calib/sens $senout fi fi if [ $tissref = "none" ]; then # whole brain M0 # in this case use the brain mask ${FSLDIR}/bin/imcp $bmask $temp_calib/refmask maskflag=1 echo "Brain mask is being used as the reference tissue (beware!)" >> $log fi if [ -z $maskflag ]; then # auto create tissue reference mask echo "FAST called to determine a reference tissue mask" >> $log # segment structural image fast -o $temp_calib/seg -p $struc # make brain mask from structural fslmaths $struc -bin $temp_calib/mask if [ $tissref = "csf" ]; then echo "Ventricle selection" >> $log # 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 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 if [ -z $stdmaskfnirt ]; 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 $struc -out $temp_calib/mask else # do std space masking using FNIRT registration # Register the structural image to standard space using FNIRT 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 # 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=$struc --in=$temp_calib/VentricleMask --warp=$temp_calib/coef_MNI2struc --out=$temp_calib/mask --interp=nn fi fi #tissue reference mask echo "Masking FAST output with standard space derrived ventricle mask" >> $log fslmaths $temp_calib/seg_pve_$fastpve -mas $temp_calib/mask $temp_calib/refmask if [ ! -z $debug ]; then imcp $temp_calib/refmask $temp_calib/refmask_high_unthresh fi if [ ! -z $transflag ]; then echo "Transforming tissue reference mask into perfusion space" #transform CSF 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/refmask --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 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 ${FSLDIR}/bin/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 # calculate M0_ref value Moval=`fslstats $temp_calib/calib -M` # this is Mz of CSF 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 # do fabber within the tissue reference mask with a sensible T1 prior mean echo "FABBER within reference tissue mask" >> $log $fabber --data=$calib --mask=$mask --output=$temp_calib/satrecov --data-order=singlefile --model=satrecov --noise=white --method=vb $tislist $llopts --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) 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 # 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*0.87" | bc` # 0.87 based on Hct, this is M0 blood (at TE=0) #echo $Moval Moval=`echo "$Moval * e(- $te / $T2b )" | bc -l` # get M0 blood at TE used echo "M0b:$Moval" echo "M0b: $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 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 $debug ]; then cp -R $temp_calib ./temp_calib fi rm -r $temp_calib echo "ASL_calib - DONE."