Initial public release; from gitlab d6d6497309e341d05c7afcb501e7d0fe6…

…05f1f64.
DCAN-Labs · Jan 7, 2020 · 8fe9f61 · 8fe9f61
1 parent c54f01d
commit 8fe9f61
Show file tree

Hide file tree

Showing 339 changed files with 101,166 additions and 0 deletions.
diff --git a/DiffusionPreprocessing/DiffPreprocPipeline.sh b/DiffusionPreprocessing/DiffPreprocPipeline.sh
@@ -0,0 +1,266 @@
+#!/bin/bash
+set -e
+
+# Preprocessing Pipeline for diffusion MRI. Generates the "data" directory that can be used as input to the fibre orientation estimation scripts.
+# Stamatios Sotiropoulos, Saad Jbabdi, Jesper Andersson, Analysis Group, FMRIB Centre, 2012.
+# Matt Glasser, Washington University, 2012.
+
+# Requirements for this script
+#  installed versions of: FSL5.0.1 or higher, FreeSurfer (version 5 or higher), gradunwarp (python code from MGH)
+#  environment: FSLDIR, FREESURFER_HOME, HCPPIPEDIR_dMRI, HCPPIPEDIR, HCPPIPEDIR_Global, HCPPIPEDIR_Bin, HCPPIPEDIR_Config, PATH (for gradient_unwarp.py)
+
+# make pipeline engine happy...
+if [ $# -eq 1 ] ; then
+    echo "Version unknown..."
+    exit 0
+fi
+
+
+########################################## Hard-Coded variables for the pipeline #################################
+
+b0dist=45     #Minimum distance in volumes between b0s considered for preprocessing
+b0maxbval=50  #Volumes with a bvalue smaller than that will be considered as b0s
+MissingFileFlag="EMPTY" #String used in the input arguments to indicate that a complete series is missing
+
+
+########################################## OUTPUT DIRECTORIES ####################################################
+
+## NB: NO assumption is made about the input paths with respect to the output directories - they can be totally different.  All input are taken directly from the input variables without additions or modifications.
+
+# Output path specifiers:
+#
+# ${StudyFolder} is an input parameter
+# ${Subject} is an input parameter
+
+# Main output directories
+# DiffFolder=${StudyFolder}/Diffusion
+# T1wDiffFolder=${StudyFolder}/T1w/Diffusion
+
+# All outputs are within the directory: ${StudyFolder}/
+# The full list of output directories are the following
+#    $DiffFolder/rawdata
+#    $DiffFolder/topup    
+#    $DiffFolder/eddy
+#    $DiffFolder/data
+#    $DiffFolder/reg
+#    $T1wDiffFolder
+
+# Also assumes that T1 preprocessing has been carried out with results in ${StudyFolder}/T1w
+
+########################################## SUPPORT FUNCTIONS #####################################################
+
+# function for parsing options
+getopt1() {
+    sopt="$1"
+    shift 1
+    for fn in $@ ; do
+	if [ `echo $fn | grep -- "^${sopt}=" | wc -w` -gt 0 ] ; then
+	    echo $fn | sed "s/^${sopt}=//"
+	    return 0
+	fi
+    done
+}
+
+defaultopt() {
+    echo $1
+}
+
+
+# function for finding the min between two numbers
+min(){
+  if [ $1 -le $2 ]; then
+     echo $1
+  else
+     echo $2
+  fi
+}
+
+
+Usage() {
+  echo "Usage: `basename $0` --posData=<dataRL1@dataRL2@...>"
+  echo "                     --negData=<dataLR1@dataLR2@...>"
+  echo "                     --path=<StudyFolder>"
+  echo "                     --subject=<SubjectID>"
+  echo "                     --echospacing=<Echo Spacing in msecs>"
+  echo "                     --PEdir=<Phase Encoding Direction (1 for LR/RL, 2 for AP/PA>"
+  echo "                     --gdcoeffs=<Coefficients for gradient nonlinearity distortion correction('NONE' to switch off)>"
+  echo "                     --printcom=<'' to run normally, 'echo' to just print and not run commands, or omit argument to run normally>"
+}
+
+################################################## OPTION PARSING ###################################################
+
+# Just give usage if no arguments specified
+if [ $# -eq 0 ] ; then Usage; exit 0; fi
+# check for correct options
+if [ $# -lt 7 ] ; then Usage; exit 1; fi
+
+# Input Variables
+PosInputImages=`getopt1 "--posData" $@`   # "$1" #dataRL1@[email protected]
+NegInputImages=`getopt1 "--negData" $@`   # "$2" #dataLR1@[email protected]
+StudyFolder=`getopt1 "--path" $@`         # "$3" #Path to subject's data folder
+Subject=`getopt1 "--subject" $@`          # "$4" #SubjectID
+echospacing=`getopt1 "--echospacing" $@`  # "$5" #Echo Spacing in msecs
+PEdir=`getopt1 "--PEdir" $@`              # "$6" #1 for LR/RL, 2 for AP/PA
+GdCoeffs=`getopt1 "--gdcoeffs" $@`        # "${7}" #Select correct coeffs for scanner gradient nonlinearities or "NONE" to turn off
+RUN=`getopt1 "--printcom" $@`             # use ="echo" for just printing everything and not running the commands (default is to run)
+RUN=${RUN:-''}
+
+# Path for scripts etc (uses variables defined in SetUpHCPPipeline.sh)
+scriptsdir=${HCPPIPEDIR_dMRI}
+globalscriptsdir=${HCPPIPEDIR_Global}
+
+# Build Paths 
+outdir=${StudyFolder}/Diffusion
+outdirT1w=${StudyFolder}/T1w/Diffusion
+if [ -d ${outdir} ]; then
+    ${RUN} rm -rf ${outdir}/rawdata
+    ${RUN} rm -rf ${outdir}/topup
+    ${RUN} rm -rf ${outdir}/eddy
+    ${RUN} rm -rf ${outdir}/data
+    ${RUN} rm -rf ${outdir}/reg
+fi
+${RUN} mkdir -p ${outdir}
+${RUN} mkdir -p ${outdirT1w}
+
+echo "OutputDir is ${outdir}"
+${RUN} mkdir ${outdir}/rawdata
+${RUN} mkdir ${outdir}/topup
+${RUN} mkdir ${outdir}/eddy
+${RUN} mkdir ${outdir}/data
+${RUN} mkdir ${outdir}/reg
+
+if [ ${PEdir} -eq 1 ]; then    #RL/LR phase encoding
+    basePos="RL"
+    baseNeg="LR"
+elif [ ${PEdir} -eq 2 ]; then  #AP/PA phase encoding
+    basePos="AP"
+    baseNeg="PA"
+fi
+
+
+########################################## DO WORK ###################################################################### 
+echo "Copying raw data"
+#Copy RL/AP images to workingdir
+PosInputImages=`echo ${PosInputImages} | sed 's/@/ /g'`
+Pos_count=1
+for Image in ${PosInputImages} ; do
+	if [[ ${Image} =~ ^.*EMPTY.*$  ]]  ;  
+	then
+		Image=EMPTY
+	fi
+
+    if [ ${Image} = ${MissingFileFlag} ];
+    then	
+        PosVols[${Pos_count}]=0
+    else
+	PosVols[${Pos_count}]=`${FSLDIR}/bin/fslval ${Image} dim4`
+	absname=`${FSLDIR}/bin/imglob ${Image}`
+	${RUN} ${FSLDIR}/bin/imcp ${absname} ${outdir}/rawdata/${basePos}_${Pos_count}
+	${RUN} cp ${absname}.bval ${outdir}/rawdata/${basePos}_${Pos_count}.bval
+	${RUN} cp ${absname}.bvec ${outdir}/rawdata/${basePos}_${Pos_count}.bvec
+    fi	
+    Pos_count=$((${Pos_count} + 1))
+done
+
+#Copy LR/PA images to workingdir
+NegInputImages=`echo ${NegInputImages} | sed 's/@/ /g'`
+Neg_count=1
+for Image in ${NegInputImages} ; do
+	if [[ ${Image} =~ ^.*EMPTY.*$  ]]  ;  
+	then
+		Image=EMPTY
+	fi
+
+    if [ ${Image} = ${MissingFileFlag} ];
+    then	
+	NegVols[${Neg_count}]=0
+    else
+	NegVols[${Neg_count}]=`${FSLDIR}/bin/fslval ${Image} dim4`
+	absname=`${FSLDIR}/bin/imglob ${Image}`
+	${RUN} ${FSLDIR}/bin/imcp ${absname} ${outdir}/rawdata/${baseNeg}_${Neg_count}
+	${RUN} cp ${absname}.bval ${outdir}/rawdata/${baseNeg}_${Neg_count}.bval
+	${RUN} cp ${absname}.bvec ${outdir}/rawdata/${baseNeg}_${Neg_count}.bvec
+    fi	
+    Neg_count=$((${Neg_count} + 1))
+done
+
+if [ ${Pos_count} -ne ${Neg_count} ]; then
+    echo "Wrong number of input datasets! Make sure that you provide pairs of input filenames."
+    echo "If the respective file does not exist, use EMPTY in the input arguments."
+    exit 1
+fi
+
+#Create two files for each phase encoding direction, that for each series contains the number of corresponding volumes and the number of actual volumes.
+#The file e.g. RL_SeriesCorrespVolNum.txt will contain as many rows as non-EMPTY series. The entry M in row J indicates that volumes 0-M from RLseries J
+#has corresponding LR pairs. This file is used in basic_preproc to generate topup/eddy indices and extract corresponding b0s for topup.
+#The file e.g. Pos_SeriesVolNum.txt will have as many rows as maximum series pairs (even unmatched pairs). The entry M N in row J indicates that the RLSeries J has its 0-M volumes corresponding to LRSeries J and RLJ has N volumes in total. This file is used in eddy_combine.
+Paired_flag=0
+for (( j=1; j<${Pos_count}; j++ )) ; do
+    CorrVols=`min ${NegVols[${j}]} ${PosVols[${j}]}`
+    ${RUN} echo ${CorrVols} ${PosVols[${j}]} >> ${outdir}/eddy/Pos_SeriesVolNum.txt
+    if [ ${PosVols[${j}]} -ne 0 ]; then
+	${RUN} echo ${CorrVols} >> ${outdir}/rawdata/${basePos}_SeriesCorrespVolNum.txt
+	if [ ${CorrVols} -ne 0 ]; then
+	    Paired_flag=1
+	fi
+    fi	
+done
+for (( j=1; j<${Neg_count}; j++ )) ; do
+    CorrVols=`min ${NegVols[${j}]} ${PosVols[${j}]}`
+    ${RUN} echo ${CorrVols} ${NegVols[${j}]} >> ${outdir}/eddy/Neg_SeriesVolNum.txt
+    if [ ${NegVols[${j}]} -ne 0 ]; then
+	${RUN} echo ${CorrVols} >> ${outdir}/rawdata/${baseNeg}_SeriesCorrespVolNum.txt
+    fi	
+done
+
+if [ ${Paired_flag} -eq 0 ]; then
+    echo "Wrong Input! No pairs of phase encoding directions have been found!"
+    echo "At least one pair is needed!"
+    exit 1
+fi
+
+echo "Running Basic Preprocessing"
+${RUN} ${scriptsdir}/basic_preproc.sh ${outdir} ${echospacing} ${PEdir} ${b0dist} ${b0maxbval}
+
+echo "Running Topup"
+${RUN} ${scriptsdir}/run_topup.sh ${outdir}/topup
+
+echo "Running Eddy"
+${RUN} ${scriptsdir}/run_eddy.sh ${outdir}/eddy
+
+GdFlag=0
+if [ ! ${GdCoeffs} = "NONE" ] ; then
+    echo "Gradient nonlinearity distortion correction coefficients found!"
+    GdFlag=1
+fi
+
+echo "Running Eddy PostProcessing"
+${RUN} ${scriptsdir}/eddy_postproc.sh ${outdir} ${GdCoeffs}
+
+#Naming Conventions
+T1wFolder="${StudyFolder}/T1w" #Location of T1w images
+T1wImage="${T1wFolder}/T1w_acpc_dc"
+T1wRestoreImage="${T1wFolder}/T1w_acpc_dc_restore"
+T1wRestoreImageBrain="${T1wFolder}/T1w_acpc_dc_restore_brain"
+BiasField="${T1wFolder}/BiasField_acpc_dc"
+FreeSurferBrainMask="${T1wFolder}/brainmask_fs"
+RegOutput="${outdir}"/reg/"Scout2T1w"
+QAImage="${outdir}"/reg/"T1wMulEPI"
+DiffRes=`${FSLDIR}/bin/fslval ${outdir}/data/data pixdim1`
+DiffRes=`printf "%0.2f" ${DiffRes}`
+
+echo "Running Diffusion to Structural Registration"
+${RUN} ${scriptsdir}/DiffusionToStructural.sh \
+  --t1folder="${T1wFolder}" \
+  --subject="${Subject}" \
+  --workingdir="${outdir}/reg" \
+  --datadiffdir="${outdir}/data" \
+  --t1="${T1wImage}" \
+  --t1restore="${T1wRestoreImage}" \
+  --t1restorebrain="${T1wRestoreImageBrain}" \
+  --biasfield="${BiasField}" \
+  --brainmask="${FreeSurferBrainMask}" \
+  --datadiffT1wdir="${outdirT1w}" \
+  --regoutput="${RegOutput}" \
+  --QAimage="${QAImage}" \
+  --gdflag=${GdFlag} --diffresol=${DiffRes}
diff --git a/DiffusionPreprocessing/scripts/DiffusionToStructural.sh b/DiffusionPreprocessing/scripts/DiffusionToStructural.sh
@@ -0,0 +1,109 @@
+#!/bin/bash
+
+set -e
+echo -e "\n START: DiffusionToStructural"
+
+
+########################################## SUPPORT FUNCTIONS ########################################## 
+
+# function for parsing options
+getopt1() {
+    sopt="$1"
+    shift 1
+    for fn in $@ ; do
+	if [ `echo $fn | grep -- "^${sopt}=" | wc -w` -gt 0 ] ; then
+	    echo $fn | sed "s/^${sopt}=//"
+	    return 0
+	fi
+    done
+}
+
+defaultopt() {
+    echo $1
+}
+
+################################################## OPTION PARSING #####################################################
+# Input Variables
+
+FreeSurferSubjectFolder=`getopt1 "--t1folder" $@`  # "$1" #${StudyFolder}/T1w
+FreeSurferSubjectID=`getopt1 "--subject" $@`       # "$2" #Subject ID
+WorkingDirectory=`getopt1 "--workingdir" $@`       # "$3" #Path to registration working dir, e.g. ${StudyFolder}/Diffusion/reg
+DataDirectory=`getopt1 "--datadiffdir" $@`         # "$4" #Path to diffusion space diffusion data, e.g. ${StudyFolder}/Diffusion/data
+T1wImage=`getopt1 "--t1" $@`                       # "$5" #T1w_acpc_dc image
+T1wRestoreImage=`getopt1 "--t1restore" $@`         # "$6" #T1w_acpc_dc_restore image
+T1wBrainImage=`getopt1 "--t1restorebrain" $@`      # "$7" #T1w_acpc_dc_restore_brain image
+BiasField=`getopt1 "--biasfield" $@`               # "$8" #Bias_Field_acpc_dc
+InputBrainMask=`getopt1 "--brainmask" $@`          # "$9" #Freesurfer Brain Mask, e.g. brainmask_fs
+GdcorrectionFlag=`getopt1 "--gdflag" $@`           # "$10"#Flag for gradient nonlinearity correction (0/1 for Off/On) 
+DiffRes=`getopt1 "--diffresol" $@`                 # "$11"#Diffusion resolution in mm (assume isotropic)
+
+# Output Variables
+T1wOutputDirectory=`getopt1 "--datadiffT1wdir" $@` # "$12" #Path to T1w space diffusion data (for producing output)
+RegOutput=`getopt1 "--regoutput" $@`               # "$13" #Temporary file for sanity checks 
+QAImage=`getopt1 "--QAimage" $@`                   # "$14" #Temporary file for sanity checks 
+
+echo $T1wOutputDirectory
+
+# Paths for scripts etc (uses variables defined in SetUpHCPPipeline.sh)
+GlobalScripts=${HCPPIPEDIR_Global}
+GlobalBinaries=${HCPPIPEDIR_Bin}
+
+T1wBrainImageFile=`basename $T1wBrainImage`
+regimg="nodif"
+
+${FSLDIR}/bin/imcp "$T1wBrainImage" "$WorkingDirectory"/"$T1wBrainImageFile"
+
+#b0 FLIRT BBR and bbregister to T1w
+${FSLDIR}/bin/epi_reg --epi="$DataDirectory"/"$regimg" --t1="$T1wImage" --t1brain="$WorkingDirectory"/"$T1wBrainImageFile" --out="$WorkingDirectory"/"$regimg"2T1w_initII
+
+${FSLDIR}/bin/applywarp --rel --interp=spline -i "$DataDirectory"/"$regimg" -r "$T1wImage" --premat="$WorkingDirectory"/"$regimg"2T1w_initII_init.mat -o "$WorkingDirectory"/"$regimg"2T1w_init.nii.gz
+${FSLDIR}/bin/applywarp --rel --interp=spline -i "$DataDirectory"/"$regimg" -r "$T1wImage" --premat="$WorkingDirectory"/"$regimg"2T1w_initII.mat -o "$WorkingDirectory"/"$regimg"2T1w_initII.nii.gz
+${FSLDIR}/bin/fslmaths "$WorkingDirectory"/"$regimg"2T1w_initII.nii.gz -div "$BiasField" "$WorkingDirectory"/"$regimg"2T1w_restore_initII.nii.gz
+
+SUBJECTS_DIR="$FreeSurferSubjectFolder"
+export SUBJECTS_DIR
+${FREESURFER_HOME}/bin/bbregister --s "$FreeSurferSubjectID" --mov "$WorkingDirectory"/"$regimg"2T1w_restore_initII.nii.gz --surf white.deformed --init-reg "$FreeSurferSubjectFolder"/"$FreeSurferSubjectID"/mri/transforms/eye.dat --bold --reg "$WorkingDirectory"/EPItoT1w.dat --o "$WorkingDirectory"/"$regimg"2T1w.nii.gz
+${FREESURFER_HOME}/bin/tkregister2 --noedit --reg "$WorkingDirectory"/EPItoT1w.dat --mov "$WorkingDirectory"/"$regimg"2T1w_restore_initII.nii.gz --targ "$T1wImage".nii.gz --fslregout "$WorkingDirectory"/diff2str_fs.mat
+
+${FSLDIR}/bin/convert_xfm -omat "$WorkingDirectory"/diff2str.mat -concat "$WorkingDirectory"/diff2str_fs.mat "$WorkingDirectory"/"$regimg"2T1w_initII.mat
+${FSLDIR}/bin/convert_xfm -omat "$WorkingDirectory"/str2diff.mat -inverse "$WorkingDirectory"/diff2str.mat
+
+${FSLDIR}/bin/applywarp --rel --interp=spline -i "$DataDirectory"/"$regimg" -r "$T1wImage".nii.gz --premat="$WorkingDirectory"/diff2str.mat -o "$WorkingDirectory"/"$regimg"2T1w
+${FSLDIR}/bin/fslmaths "$WorkingDirectory"/"$regimg"2T1w -div "$BiasField" "$WorkingDirectory"/"$regimg"2T1w_restore
+
+#Are the next two scripts needed?
+${FSLDIR}/bin/imcp "$WorkingDirectory"/"$regimg"2T1w_restore "$RegOutput"
+${FSLDIR}/bin/fslmaths "$T1wRestoreImage".nii.gz -mul "$WorkingDirectory"/"$regimg"2T1w_restore.nii.gz -sqrt "$QAImage"_"$regimg".nii.gz
+
+#Generate 1.25mm structural space for resampling the diffusion data into
+${FSLDIR}/bin/flirt -interp spline -in "$T1wRestoreImage" -ref "$T1wRestoreImage" -applyisoxfm ${DiffRes} -out "$T1wRestoreImage"_${DiffRes}
+${FSLDIR}/bin/applywarp --rel --interp=spline -i "$T1wRestoreImage" -r "$T1wRestoreImage"_${DiffRes} -o "$T1wRestoreImage"_${DiffRes}
+
+#Generate 1.25mm mask in structural space
+${FSLDIR}/bin/flirt -interp nearestneighbour -in "$InputBrainMask" -ref "$InputBrainMask" -applyisoxfm ${DiffRes} -out "$T1wOutputDirectory"/nodif_brain_mask
+${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/nodif_brain_mask -kernel 3D -dilM "$T1wOutputDirectory"/nodif_brain_mask
+
+#Rotate bvecs from diffusion to structural space
+${GlobalScripts}/Rotate_bvecs.sh "$DataDirectory"/bvecs "$WorkingDirectory"/diff2str.mat "$T1wOutputDirectory"/bvecs
+cp "$DataDirectory"/bvals "$T1wOutputDirectory"/bvals
+
+#Register diffusion data to T1w space. Account for gradient nonlinearities if requested
+if [ ${GdcorrectionFlag} -eq 1 ]; then
+    echo "Correcting Diffusion data for gradient nonlinearities and registering to structural space"
+    #In the future, we want this applywarp to be part of eddy and avoid second resampling step.
+    ${FSLDIR}/bin/convertwarp --rel --relout --warp1="$DataDirectory"/warped/fullWarp --postmat="$WorkingDirectory"/diff2str.mat --ref="$T1wRestoreImage"_${DiffRes} --out="$WorkingDirectory"/grad_unwarp_diff2str
+    ${FSLDIR}/bin/applywarp --rel -i "$DataDirectory"/warped/data_warped -r "$T1wRestoreImage"_${DiffRes} -w "$WorkingDirectory"/grad_unwarp_diff2str --interp=spline -o "$T1wOutputDirectory"/data
+
+    #Now register the grad_dev tensor 
+    ${GlobalBinaries}/vecreg -i "$DataDirectory"/grad_dev -o "$T1wOutputDirectory"/grad_dev -r "$T1wRestoreImage"_${DiffRes} -t "$WorkingDirectory"/diff2str.mat --interp=spline
+    ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/grad_dev -mas "$T1wOutputDirectory"/nodif_brain_mask "$T1wOutputDirectory"/grad_dev  #Mask-out values outside the brain 
+else
+#Register diffusion data to T1w space without considering gradient nonlinearities
+    ${FSLDIR}/bin/flirt -in "$DataDirectory"/data -ref "$T1wRestoreImage"_${DiffRes} -applyxfm -init "$WorkingDirectory"/diff2str.mat -interp spline -out "$T1wOutputDirectory"/data
+fi
+
+${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/data -mas "$T1wOutputDirectory"/nodif_brain_mask "$T1wOutputDirectory"/data  #Mask-out data outside the brain 
+${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/data -thr 0 "$T1wOutputDirectory"/data      #Remove negative intensity values (caused by spline interpolation) from final data
+
+
+echo " END: DiffusionToStructural"