FreeSurferPipeline

#!/bin/bash 
set -e

scriptDir=$( dirname "$0" )
source "${scriptDir}/preproc_functions/helper_functions"
logFile="freesurfer.log"

# Requirements for this script
#  installed versions of: FSL (version 5.0.6), FreeSurfer (version 5.3.0-HCP)
#  environment: FSLDIR , FREESURFER_HOME , HCPPIPEDIR , CARET7DIR 


show_usage() {
cat <<EndOfHelp
----------------------------------------------
FreeSurferPipeline is an adaptation of the HCP FreeSurfer pipeline described in Glasser et al. 2013

Command line options:
  -log: name of logfile. Default: freesurfer.log
  -ncpus: number of CPUs to use in recon-all main processing. Default: 1
  -T1: bias-corrected, gradient-undistorted T1-weighted image [Required]
  -T1brain: skull-stripped version of -T1 [Required]
  -subject: Subject ID to be used as identifier to FreeSurfer [Required]
  -subjectDir: Master directory for FreeSurfer outputs (will place within a folder named by Subject ID) [Required]
  -T2: T2-weighted image to combine with T1 (not vetted at this point)
  -seed: random seed to be used for recon-all (ensures completely consistent results across re-runs)

Example call: FreeSurferPipeline -T1 mprage_biascorr.nii.gz -T1brain mprage_bet.nii.gz -subject 1005 -subjectDir ~/FS_Subjects
----------------------------------------------
EndOfHelp

}

################################################## OPTION PARSING #####################################################

datefmt='+%F+%I:%M'
thiscommandinfo="$0 $@\nRun started $(date $datefmt)"

# write to an incomplete file (later to be changed to complete if process finishes)
ncpus=1

#if no parameters are passed in, then print help and exit.
if [ $# -eq 0 ]; then
    show_usage
    exit 0
fi


while [ _$1 != _ ] ; do
    if [[ $1 = -log ]]; then
	logFile="${2}"
	shift 2
    elif [[ $1 = -ncpus ]]; then
	ncpus="${1}"
	shift 2
    elif [[ $1 = -t1 || $1 = -T1 ]] ; then
	T1wImage="${2}"
	shift 2
    elif [[ $1 = -t1brain || $1 = -T1brain ]]; then
	T1wImageBrain="${2}"
	shift 2
    elif [[ $1 = -t2 || $1 = -T2 ]]; then
	T2wImage="${2}"
	shift 2
    elif [[ $1 = -subject ]]; then
	SubjectID="${2}"
	shift 2
    elif [[ $1 = -subjectDir ]]; then
	SubjectDIR="${2}"
	shift 2
    elif [[ $1 = -seed ]]; then
	recon_all_seed="${2}"
	shift 2
    else
	echo -e "----------------\n\n"
	echo "Unrecognized command line parameter: ${1}"
	exit 1
    fi
done

if [ -z "$T1wImage" ]; then
    echo "Missing required argument: -t1 (bias-corrected structural scan)"
    exit 1
fi

if [ ! -r "$T1wImage" ]; then
    echo "Unable to read file: $T1wImage"
    exit 1
fi

if [ -z "$T1wImageBrain" ]; then
    echo "Missing required argument: -t1brain (skull-stripped structural scan)"
    exit 1
fi

if [ ! -r "$T1wImageBrain" ]; then
    echo "Unable to read file: $T1wImageBrain"
    exit 1
fi

rel "$thiscommandinfo" c

# Input Variables
#SubjectID=`opts_GetOpt1 "--subject" $@` #FreeSurfer Subject ID Name
#SubjectDIR=`opts_GetOpt1 "--subjectDIR" $@` #Location to Put FreeSurfer Subject's Folder
#T1wImage=`opts_GetOpt1 "--t1" $@` #T1w FreeSurfer Input (Full Resolution)
#T1wImageBrain=`opts_GetOpt1 "--t1brain" $@` 
#T2wImage=`opts_GetOpt1 "--t2" $@` #T2w FreeSurfer Input (Full Resolution)
#recon_all_seed=`opts_GetOpt1 "--seed" $@`

# ------------------------------------------------------------------------------
#  Show Command Line Options
# ------------------------------------------------------------------------------

rel "Finished Parsing Command Line Options" c
rel "SubjectID: ${SubjectID}" c
rel "SubjectDIR: ${SubjectDIR}" c
rel "T1wImage: ${T1wImage}" c
rel "T1wImageBrain: ${T1wImageBrain}" c
rel "T2wImage: ${T2wImage}" c
rel "recon_all_seed: ${recon_all_seed}" c

# figure out whether to include a random seed generator seed in all the recon-all command lines
seed_cmd_appendix=""
if [ -z "${recon_all_seed}" ] ; then
	seed_cmd_appendix=""
else
	seed_cmd_appendix="-norandomness -rng-seed ${recon_all_seed}"
fi
rel "seed_cmd_appendix: ${seed_cmd_appendix}" c


# ------------------------------------------------------------------------------
#  Identify Tools
# ------------------------------------------------------------------------------

which_flirt=$( which flirt )
flirt_version=$( flirt -version )
rel "which flirt: ${which_flirt}" c
rel "flirt -version: ${flirt_version}" c

which_applywarp=$( which applywarp )
rel "which applywarp: ${which_applywarp}" c

which_fslstats=$( which fslstats )
rel "which fslstats: ${which_fslstats}" c

which_fslmaths=$( which fslmaths )
rel "which fslmaths: ${which_fslmaths}" c

which_recon_all=$( which recon-all )
recon_all_version=$( recon-all --version )
rel "which recon-all: ${which_recon_all}" c
rel "recon-all --version: ${recon_all_version}" c

which_mri_convert=$( which mri_convert )
rel "which mri_convert: ${which_mri_convert}" c

which_mri_em_register=$( which mri_em_register )
mri_em_register_version=$( mri_em_register --version )
rel "which mri_em_register: ${which_mri_em_register}" c
rel "mri_em_register --version: ${mri_em_register_version}" c

which_mri_watershed=$( which mri_watershed )
mri_watershed_version=$( mri_watershed --version )
rel "which mri_watershed: ${which_mri_watershed}" c
rel "mri_watershed --version: ${mri_watershed_version}" c

#figure out if we are using freesurfer 6, which requires different steps
[ $( cat $FREESURFER_HOME/build-stamp.txt )  == "freesurfer-Linux-centos6_x86_64-stable-pub-v6.0.0-2beb96c" ] && fs6=1 || fs6=0

# Start work

T1wImageFile=`remove_ext $T1wImage`;
T1wImageBrainFile=`remove_ext $T1wImageBrain`;

#PipelineScripts=${HCPPIPEDIR_FS}

if [ -e "$SubjectDIR"/"$SubjectID"/scripts/IsRunning.lh+rh ] ; then
  rm "$SubjectDIR"/"$SubjectID"/scripts/IsRunning.lh+rh
fi

##Make Spline Interpolated Downsample to 1mm
##rel "Make Spline Interpolated Downsample to 1mm" c

Mean=$( fslstats $T1wImageBrain -M ) #normalize image based on mean intensity

#HCP downsamples to 1mm to be compatible with FreeSurfer. Because our images are already at 1mm, no need to do this
#flirt -interp spline -in "$T1wImage" -ref "$T1wImage" -applyisoxfm 1 -out "$T1wImageFile"_1mm.nii.gz
#applywarp --rel --interp=spline -i "$T1wImage" -r "$T1wImageFile"_1mm.nii.gz --premat=$FSLDIR/etc/flirtsch/ident.mat -o "$T1wImageFile"_1mm.nii.gz
#applywarp --rel --interp=nn -i "$T1wImageBrain" -r "$T1wImageFile"_1mm.nii.gz --premat=$FSLDIR/etc/flirtsch/ident.mat -o "$T1wImageBrainFile"_1mm.nii.gz

#mean normalize
fslmaths "$T1wImageFile" -div $Mean -mul 150 -abs "${T1wImageFile}_normalized"

#Initial Recon-all Steps
rel "Initial Recon-all Steps" c

# Call recon-all with flags that are part of "-autorecon1", with the exception of -skullstrip.
# -skullstrip of FreeSurfer not reliable for Phase II data because of poor FreeSurfer mri_em_register registrations with Skull on, 
# so run registration with PreFreeSurfer masked data and then generate brain mask as usual.
rel "recon-all -i \"${T1wImageFile}_normalized.nii.gz\" -subjid $SubjectID -sd $SubjectDIR -motioncor -talairach -nuintensitycor -normalization ${seed_cmd_appendix}"

# Generate brain mask
rel "mri_convert \"${T1wImageBrainFile}.nii.gz\" \"$SubjectDIR\"/\"$SubjectID\"/mri/brainmask.mgz --conform"
if [[ $fs6 -eq 1 ]]; then
    rel "mri_em_register -mask \"$SubjectDIR\"/\"$SubjectID\"/mri/brainmask.mgz \"$SubjectDIR\"/\"$SubjectID\"/mri/nu.mgz $FREESURFER_HOME/average/RB_all_2016-05-10.vc700.gca \"$SubjectDIR\"/\"$SubjectID\"/mri/transforms/talairach_with_skull.lta"
    rel "mri_watershed -T1 -brain_atlas $FREESURFER_HOME/average/RB_all_withskull_2016-05-10.vc700.gca \"$SubjectDIR\"/\"$SubjectID\"/mri/transforms/talairach_with_skull.lta \"$SubjectDIR\"/\"$SubjectID\"/mri/T1.mgz \"$SubjectDIR\"/\"$SubjectID\"/mri/brainmask.auto.mgz"  
else
    #assume 5.x
    rel "mri_em_register -mask \"$SubjectDIR\"/\"$SubjectID\"/mri/brainmask.mgz \"$SubjectDIR\"/\"$SubjectID\"/mri/nu.mgz $FREESURFER_HOME/average/RB_all_2008-03-26.gca \"$SubjectDIR\"/\"$SubjectID\"/mri/transforms/talairach_with_skull.lta"
    rel "mri_watershed -T1 -brain_atlas $FREESURFER_HOME/average/RB_all_withskull_2008-03-26.gca \"$SubjectDIR\"/\"$SubjectID\"/mri/transforms/talairach_with_skull.lta \"$SubjectDIR\"/\"$SubjectID\"/mri/T1.mgz \"$SubjectDIR\"/\"$SubjectID\"/mri/brainmask.auto.mgz"
fi

rel "cp \"$SubjectDIR\"/\"$SubjectID\"/mri/brainmask.auto.mgz \"$SubjectDIR\"/\"$SubjectID\"/mri/brainmask.mgz"

if [ $ncpus -gt 1 ]; then
    openmp_appendix="-openmp ${ncpus}"
fi

if [ $fs6 -eq 1 ]; then
    #depart from HCP because it leads to crashes at this stage under FS6
    #instead, just run autorecon2 and autorecon3
    rel "recon-all -subjid $SubjectID -sd $SubjectDIR -autorecon2 ${openmp_appendix} ${seed_cmd_appendix}"
    rel "recon-all -subjid $SubjectID -sd $SubjectDIR -autorecon3 ${openmp_appendix} ${seed_cmd_appendix}"
else
    rel "Call recon-all to run most of the -autorecon2 stages, but turning off smooth2, inflate2, curvstats, and segstats stages" c
    rel "recon-all -subjid $SubjectID -sd $SubjectDIR -autorecon2 -nosmooth2 -noinflate2 -nocurvstats -nosegstats ${openmp_appendix} ${seed_cmd_appendix}"

    #Highres white stuff and Fine Tune T2w to T1w Reg
    if [ -n "$T2wImage" ]; then
	rel "High resolution white matter and fine tune T2w to T1w registration" c
	echo "UNTESTED... will now blow up!"
	"$PipelineScripts"/FreeSurferHiresWhite.sh "$SubjectID" "$SubjectDIR" "$T1wImage" "$T2wImage"
    fi

    #Intermediate Recon-all Steps
    rel "Intermediate Recon-all Steps" c
    rel "recon-all -subjid $SubjectID -sd $SubjectDIR -smooth2 -inflate2 -curvstats -sphere -surfreg -jacobian_white -avgcurv -cortparc ${openmp_appendix} ${seed_cmd_appendix}"

    #Highres pial stuff (this module adjusts the pial surface based on the the T2w image)
    if [ -n "$T2wImage" ]; then
	rel "High Resolution pial surface" c
	echo "UNTESTED... will now blow up!"
	"$PipelineScripts"/FreeSurferHiresPial.sh "$SubjectID" "$SubjectDIR" "$T1wImage" "$T2wImage"
    fi

    #Final Recon-all Steps
    rel "Final Recon-all Steps" c
    rel "recon-all -subjid $SubjectID -sd $SubjectDIR -surfvolume -parcstats -cortparc2 -parcstats2 -cortparc3 -parcstats3 -cortribbon -segstats -aparc2aseg -wmparc -balabels -label-exvivo-ec ${openmp_appendix} ${seed_cmd_appendix}"
fi

rel "Completed FreeSurferPipeline for subject: $SubjectID" c