preprocessMprage

#!/bin/bash
function printHelp() {

set -e #stop script on error

cat <<EndOfHelp
----------------------------------------------
preprocessMprage is a shell script that computes the nonlinear warp of a
  participant's structural scan to a standard stereotaxic template (e.g., MNI).

It is intended to be run within an mprage directory containing raw data (or a
  nifti of the raw data if you use the -n option).

The basic pipeline is:
  1) Convert dicom files to NIFTI
  2) Optional bias field correction to reduce spatial variation in intensity
       due to coil sensitivity. Particularly important for 32-channel T1 images.
  3) Brain extract the structural image
  4) Warp structural to standard stereotactic space using affine (linear)
      transformation: flirt. (Used only for warp coefficients)
  5) Warp structural image to stereotactic space using nonlinear transformation
       with fnirt. Affine coefficients from the linear warping are used as
       starting values (helps optimization). In the case of MNI, uses FSL's
       config file with optimal parameters for fnirt. In the case of Talairach,
       use a custom settings file adapted from FSL's MNI configuration.

Command line options:
  -b -bet_opts: bet options. Passes value to fsl bet (brain extraction).
        Must be quoted (e.g., -b "-f 0.7"). Defaults to "-R -f 0.5 -v"
  -bright_skull: apply a correction to reduce the intensity of bright voxels in the skull prior to FNIRT.
        Skull voxels that are bright (esp. near top of brain) can lead to artifactual stretching in FNIRT.
        This correction is recommended if the skull is noticeably brighter (2-3x) than brain voxels.
  -check_dependencies: what tools are avaible on this system. Exit w/ status 0 if everything needed is found
  -cleanup: cleanup intermediate files not needed for further processing.
  -cleanup_only: run cleanup of intermediate files, then exit. (after checking files)
  -custom_brainmask: a NIfTI image containing a custom brain mask to be used instead of running a skull strip internally.
        Must be the same dimensions as the structural scan. Technically, this can even be a skull-stripped volume such that
        non-zero values are used to mask the structural image.
  -cut_zindex Z:  keep index Z to end. use Z1-Z2 to specify range. (e.g. '0-100' to keep up until 100 instead of 100 and up)
        useful when brain is only part of image.
  -d -delete_dicom: delete or archive DICOM files. Options are -d n ("no": leaves DICOM
        untouched), -d y ("yes": deletes DICOM files), -d a ("archive":
        compresses files into archive file: mprage_dicom.tar.gz. If not
        passed, user will be prompted for action. Can also use "yes", "no", and "archive."
  -deface: use pydeface to remove face and maybe also hands and arms near the face (esp. babies)
  -deneck: remove from the neck down with local remove_neck.py (implemented for baby data)
  -fnirt_mask: controls how dilated the --refmask is for fnirt. Changing this may help with undesirable stretching at edge. Options:
  	vtight: no dilation of template brain mask
	tight: 1x dilation (-dilF) of mask
	normal: 2x dilation (-dilF -dilF) of mask (default)
	loose: 3x dilation (-dilF -dilF -dilF) of mask
  -grad_unwarp: a file containing gradient nonlinearity coefficients for a SIEMENS scanner (e.g., trio.coeff.grad).
        Structural image will be corrected for gradient distortion prior to transformation to template.
  -h -help: print command help
  -log: Name for log file that documents each command that is run. Default: preprocessMprage.log
  -n -nifti: skip DICOM conversion and start with specified NIFTI file. 
        Example: -n mprage.nii.gz (will complete bet, flirt, and
        fnirt, but skip DICOM -> NIFTI conversion).
  -no_bias: skip bias field correction.
  -no_robustfov: skip robustfov neck cleanup
  -o -output: output file name. The filename for the final warped mprage image.
        Defaults to mprage_nonlinear_warp_<referenceName>.nii.gz
  -p -dicom: file pattern for dicom MR files. Defaults to "MR*". Please enclose the pattern in quotes or it may not function properly
  -post_bet_skullmask: provide a mask to remove extra skull even after skull stripping. likely created by hand using afni's draw plugin
        See https://labneurocogdevel.github.io/tuts/afni_draw
  -r -template_brain: reference brain. Currently supports "MNI_2mm", "MNI_FSL_2mm", "SPM_2mm", "Tal_2mm", 1YO_2mm, neo_2mm.
        Default is "MNI_2mm".
  -startover: ignore .preprocessmprage_complete. useful if changing/adding arguments like "-ss_method bet -b '-f .6 -R'"
  -ss_method: skull-stripping method used to create final brain-extracted data. Default is fnirt-inv. Options:
        fnirt-inv: uses the template brain mask warped onto the subject's brain to define brain mask
        bet: uses the bet program (FSL) to skull strip (also respects -bet_opts).
        3dSkullStrip: uses the 3dSkullStrip program (can pass additional parameters using -bet_opts).
        noisy: use @NoisySkullStrip
        ROBEX: uses the ROBEX program (which doesn't accept any parameters).
        freesurfer=/path/to/subj/mri/brainmask.mgz
  -strong_bias: used for images with very strong bias fields
  -unifize:  run afni's 3dUnifize on input mprage
  -use_old_mni: support legacy preprocessing
  -w -warpres: resolution of warp basis (nonlinear neighborhood size) in mm. Default is 8 mm.
  -weak_bias: used for images with little and/or smooth bias fields (cf. fsl_anat). This is the default.

Example call: preprocessMprage -r MNI_2mm -b "-R -f 0.5 -g 0.2" -d a -o mprage_final.nii.gz -w 6
----------------------------------------------
EndOfHelp

} #end of printHelp

#Author: Michael Hallquist
#Written: 5/2/2010
#Last defaults updated: 7/21/2017
#
#Changelog:
# 11/05/2021(WF)
#  - add ss_method option freesurfer=$SUBJECT_DIR/mri/brainmask.mgz (b/c 7T GRAPPA T1 are difficult)
# 09/07/2021(WF)
#  - add DIL_FNIRT_INPUT_BY env variable. testing nonlinear warp with baby brains
# 17/06/2021(WF)
#  - add '-deneck' and '-deface' for OHSU baby dataset
#    added preproc_functions/mprage_utils. moved reorient to there
#    '.cur_step' file tracks "pre-preprocess" steps, things before skullstripping
# 12/05/2021(WF)
#  - update version reporting to include git+dirty/clean & user@host:pwd
# 10/05/2021(WF)
#  - '-post_bet_skullmask' added for applying hand drawn skull mask 
#22/04/2021 (WF) 
# - '-unifize' switch runs 3dUnifize (remove shading artificat)
# - save mprage nifti notes to 'mprage_history.txt' (e.g. lncdtools/mknii)
#11/0l/2018 (WF) 
# - can use octave, add old_template_check (-use_old_mni)
#10/29/2018 (WF) 
# - add -check_dependencies to check install
#7/21/2017 (MH)
# - Added a correction (-bright_skull) for reducing bright skull voxels prior to FNIRT. This was leading to distortion (stretching) in some cases
# - Added an option (-fnirt_mask) to choose the tightness of the refmask for fnirt to improve warp in some subjects
# - Revert to 2x dilated mask for --refmask in fnirt. This fits FEAT setup and tends to give better results (less edge stretching) for most subjects
#7/11/2017 (MH)
# - Add more QA images and create qa_images folder
# - Default to weak bias field correction. Strong was harming fnirt in some datasets and fsl_anat has also switched to weak by default
# - Support ROBEX for brain extraction (outperforms bet and 3dSkullStrip in general)
#11/14/2016
# - Added robustfov into pipeline to help with neck cleanup on larger structural scans
# - Allow bet or 3dSkullStrip to carry forward as brain-extracted image for subsequent processing (-ss_method)
# - Allow for a custom brain mask to be used to help with manual tweaks of skull strip (-custom_brainmask)
# - Allow for a .preproc_override file to supplant parameters passed on the command line. Useful for 
#     manual custom processing to still be picked up by preprocessMprage being run in large batch (autopreproc).
# - Implement .preprocessmprage_incomplete, consistent with preprocessFunctional approach
# - Archive preprocessMprage.log on re-run, rather than deleting it (consistent with preprocessFunctional)
# - Support gradient nonlinearity distortion correction: -grad_unwarp
#10/30/2014
# - added check for initial file extension for NIfTI to make sure that it matches FSLOUTPUTYPE. Otherwise, pipeline may fail to use reoriented file.
#4/5/2014
# - added -log parameter to keep a log of processing steps.
# - generate brain mask from inverse MNI warp to reduce sensitivity to BET (follows fsl_anat)
# - generate struct_to_template.png image to display quality of struct -> template warp.
# - provide option for bias field correction (to reduce spatial variability in intensity due to coil sensitivity).
#     Bias field is also incorporated into fnirt to improve warp to template. (adapted from fsl_anat).
# - use fslreorient2std to ensure that mprage matches template (MNI) orientation (LPI/RPI)
# - provide -w option for setting warp resolution. Default to 8mm.
# - align related arguments between preprocessMprage and preprocessFunctional (e.g., -dicom).
# - switch to remove_ext for T1
# - switch to more general arg parsing to allow for longer names
#6/4/2012
# - provide option to use SPM8 canonical MNI template
#10/27/2011
# - put back in 3dresample to orient to LPI. With some of Kirsten's data, warp was failing due to some orientation problem.
#10/18/2011
# - removed 3dresample command to orient original functional to LPI. Bug in AFNI that was flipping storage but not header.
# - Added -f to gzip to overwrite old file.
#08/31/2011
# - Changed default MNI template to new nonlinear template from Vladimir Fonov
#07/28/2011
# - Cleaned up DICOM to NIFTI conversion sections to adopt Dimon
# - Removed sleep commands (just used for monitoring output)
#06/28/2011:
# - Changed Dimon parameters used for DICOM > NIFTI conversion
#   it was failing on a couple datasets
#5/4/2011:
# - Switched from using the oriented output file from dcm2nii;
#   its usage is counterindicated (had messed up qform/sform values)
# - Changed the code that determines which nii.gz file to use as anatomical.
#   now uses regex to determine correct file 
#3/10/2011:
# - Switched default template brain to MNI.
# - Exit script on unrecognized template brain parameter.

set -e #exit if any error occurs (stop processing)

scriptDir=$( dirname "$0" )

source "${scriptDir}/preproc_functions/helper_functions" #contains rel for running/logging; also sets up stddir env variable, preproc_git_ver
source "${scriptDir}/preproc_functions/parse_args" #contains function for parsing command line args
source "${scriptDir}/preproc_functions/template_funcs" # old_template_check
source "${scriptDir}/preproc_functions/mprage_utils" # reorient, unifize, deface, deneck, current_step_file


##load necessary modules if running on cluster (specific to ACI at the moment)
if command -v module >/dev/null && uname -a | grep -q aci.ics.psu.edu ; then
    #setup DEPENd lab environment and programs
    source /gpfs/group/mnh5174/default/lab_resources/ni_path.bash
elif command -v module >/dev/null && uname -a | grep -q its.unc.edu ; then
    source /proj/mnhallqlab/lab_resources/ni_path.bash
fi

## NB. these defaults might be overwritten when parse_mprage_args is called
#      (even if the user doesn't call any flag)
#set defaults for dicomPattern, reference, betOpts, and outputFile
deface=0
deneck=0
bright_skull=0
cleanup=0
cleanup_only=0
CUT_ZIDX=""
biasCorrect=1
strongBias=0 #0 for weak bias correction, 1 for strong bias
dicomPattern="MR*"
logFile="preprocessMprage.log"
reference="MNI_2mm"
wr=8 #warp resolution in mm
fnirt_mask="normal"
robustfov=1
ssmethod=fnirt-inv
gcoeffs= #no gradient nonlinearity correction by default

dotfile=.preprocessmprage
unifize=0
post_bet_skullmask="" # 20210510. hand drawn mask of troublesome skull

#exit script if processing complete
if [[ -f "${dotfile}_complete" && "$@" != *-startover* && "$@" != *-cleanup_only* && "$@" != *-help* ]]; then
    echo -e "\n--- preprocessMprage already complete.\n"
    exit 0
fi

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

datefmt='+%F+%I:%M'
#grab command, arguments, script version (date), and start time.
#needs to come before while loop because the positional parameters get popped off.
thiscommandinfo="$0 $@ \n# preprocessMprage v$(perl -ne 'print $1 if /^#Last defaults updated: (.*)/' $0) git $(preproc_git_ver) \n# Run started $(date $datefmt)\n# $(whoami)@$(hostname):$(pwd)"

# write to an incomplete file (later to be changed to complete if process finishes)
echo -e "$thiscommandinfo " > ${dotfile}_incomplete

#now in preproc_functions/parse_args
parse_mprage_args "$@" #parse command line inputs

# make sure we are using the correct template
old_template_check

#if a .preproc_override file exists, parse args in that file, overriding the command line args in preceding step
if [ -f .preproc_override ]; then
    #there is a challenge of resplitting quoted arguments after read. This solution (eval) worked for this case
    #http://stackoverflow.com/questions/17338863/split-a-string-stored-in-a-variable-into-multiple-words-using-spaces-but-not-t
    read override_string < .preproc_override || echo "Did you forget a trailing newline in .preproc_override?"
    echo "Overriding the following command line arguments from .preproc_override: $override_string"
    eval override=($override_string)
    parse_mprage_args "${override[@]}"
fi

#dumpenv #output all variables to screen for debugging

case "$fnirt_mask" in
    vtight) refMaskSuffix=;;
    tight) refMaskSuffix=_dil;;
    normal) refMaskSuffix=_dil2;;
    loose) refMaskSuffix=_dil3;;
    *) echo -e "-fnirt_mask argument $fnirt_mask not recognized. Options are tight, normal, loose, vloose."; exit 1 ;;
esac
	   
case "$reference" in
    Tal_2mm)  bettedRefBrain=${stddir}/talairach_fsl_mni152/TalFSL_MNI152_T1_2mm_brain;
	unbettedRefBrain=${stddir}/talairach_fsl_mni152/TalFSL_MNI152_T1_2mm;
	refMask=${stddir}/talairach_fsl_mni152/TalFSL_MNI152_T1_2mm_brain_mask${refMaskSuffix}; #TODO: support other masks
	fnirtConfig=${stddir}/fnirtTalairachSettings.cnf;;
    MNI_1mm) bettedRefBrain=${stddir}/mni_icbm152_nlin_asym_09c/mni_icbm152_t1_tal_nlin_asym_09c_brain;
	unbettedRefBrain=${stddir}/mni_icbm152_nlin_asym_09c/mni_icbm152_t1_tal_nlin_asym_09c;
	refMask=${stddir}/mni_icbm152_nlin_asym_09c/mni_icbm152_t1_tal_nlin_asym_09c_mask_dil; #Jul2017: 3dmask_tool dil of brain mask
	fnirtConfig=${stddir}/fnirtMNI1mm.cnf;; # DONT TRUST THIS FILE # TODO: check/fix
    MNI_2mm) bettedRefBrain=${stddir}/mni_icbm152_nlin_asym_09c/mni_icbm152_t1_tal_nlin_asym_09c_brain_2mm;
	unbettedRefBrain=${stddir}/mni_icbm152_nlin_asym_09c/mni_icbm152_t1_tal_nlin_asym_09c_2mm;
	refMask=${stddir}/mni_icbm152_nlin_asym_09c/mni_icbm152_t1_tal_nlin_asym_09c_mask_2mm${refMaskSuffix};
	fnirtConfig=${FSLDIR}/etc/flirtsch/T1_2_MNI152_2mm.cnf;;
    MNI_FSL_2mm) bettedRefBrain=${stddir}/fsl_mni152/MNI152_T1_2mm_brain;
	unbettedRefBrain=${stddir}/fsl_mni152/MNI152_T1_2mm;
	refMask=${stddir}/fsl_mni152/MNI152_T1_2mm_brain_mask${refMaskSuffix};
	fnirtConfig=${FSLDIR}/etc/flirtsch/T1_2_MNI152_2mm.cnf;;
    SPM_2mm)  bettedRefBrain=${stddir}/spm8_mni/T1_brain;
	unbettedRefBrain=${stddir}/spm8_mni/T1;
	refMask=${stddir}/spm8_mni/brainmask_0.5thresh;
	fnirtConfig=${FSLDIR}/etc/flirtsch/T1_2_MNI152_2mm.cnf;;
    # 20210317WF for LC 
    1YO_2mm)  bettedRefBrain=${stddir}/UNCInfant/2mm/infant-1yr-brain
	unbettedRefBrain=${stddir}/UNCInfant/2mm/infant-1yr-withSkull
	refMask=${stddir}/UNCInfant/2mm/infant-1yr-brainmask;
	fnirtConfig=${FSLDIR}/etc/flirtsch/T1_2_MNI152_2mm.cnf;;
    # 20210416WF
    neo_2mm)  bettedRefBrain=${stddir}/UNCInfant/2mm_neo/infant-neo-brain
	unbettedRefBrain=${stddir}/UNCInfant/2mm_neo/infant-neo-withSkull
	refMask=${stddir}/UNCInfant/2mm_neo/infant-neo-brainmask;
	fnirtConfig=${FSLDIR}/etc/flirtsch/T1_2_MNI152_2mm.cnf;;
    *) echo -e "Reference brain '$reference' not recognized. Options are MNI_2mm, MNI_FSL_2mm, SPM_2mm, Tal_2mm, 1YO_2mm, neo_2mm."; exit 1 ;;
esac

echo "settings:
  preproc_override?  '$override_string'
  reference    = '$reference'
  ssmethod     = '$ssmethod'
  betOpts      = '$betOpts'
  strongBias   = '$strongBias'
  cleanup      = '$cleanup'
  bright_skull = '$bright_skull'
  unifize      = '$unifize'
  skullmask    = '$post_bet_skullmask'
  deface/deneck= '$deface'/'$deneck'
  brainmask    = '$brainmask'
  CUT_ZIDX     = '$CUT_ZIDX'
  "

[ -n "$DRYRUN" ] && exit 1


#define warped brain mask without dilation (used for creating a final betted image)
nodil=$( echo $refMask | perl -pe 's/(_dil[0-9]*)*(\.nii)*(\.gz)*$//' ) #remove any dil, as well as file extension

if [[ ! -f ${bettedRefBrain}.nii && ! -f ${bettedRefBrain}.nii.gz  ]]; then
    echo -e "Skull stripped reference brain not found: $bettedRefBrain\n"
    exit 1
fi

#if unbetted reference is set, but file does not exist, throw error
if [[ ! -f ${unbettedRefBrain}.nii && ! -f ${unbettedRefBrain}.nii.gz ]]; then
    echo -e "Reference brain not found: $unbettedRefBrain\n"
    exit 1
fi

#check for the fnirt config file
if [ ! -f $fnirtConfig ]; then
    echo -e "FNIRT config file not found: $fnirtConfig\n"
    exit 1
fi

#if nifti passed in (i.e., skip dicom to nifti, then verify its existence
if [[ -n $nifti && ! -f $nifti ]]; then
    echo -e "Instructed to start preprocessing with nifti (-n), but file not found.\nFile:${nifti}\n"
    exit 1
fi

if [ -n "$post_bet_skullmask" -a ! -r "$post_bet_skullmask" ]; then
   echo "set -post_bet_skullmask but '$post_bet_skullmask' does not exist (in $(pwd))"
   exit 1
fi

#figure out file extension for FSL programs
if [ -z $FSLOUTPUTTYPE ]; then
    export FSLOUTPUTTYPE=NIFTI_GZ
fi

if [ $FSLOUTPUTTYPE = NIFTI_GZ ]; then
    ext=".nii.gz"
elif [ $FSLOUTPUTTYPE = NIFTI ]; then
    ext=".nii"
else
    echo "Not setup to handle FSLOUTPUTTYPE: $FSLOUTPUTTYPE."
    exit 1
fi

if [ -n "${logFile}" ]; then
    if [ -f "${logFile}" ]; then
	#stat and date are not portable across Linux and BSD...
        [ $( uname ) = Darwin ] && mtime=$( stat -f "%Sm" -t "%Y%m%d_%H%M" "${logFile}" ) || mtime=$( date -r "${logFile}" +%Y%m%d_%H%M )
	mv "${logFile}" "${logFile}_${mtime}"
    fi
    
    #add absolute path to log file location
    logFile="$(pwd)/$( basename $logFile )" ##TODO: make path handling more robust to non-local directories

    echo "#!/bin/bash" > "${logFile}"
    echo "## Log of preprocessMprage commands" >> "${logFile}"
    echo -e "## Call: $thiscommandinfo" >> "${logFile}"
    [ -n "$override_string" ] && echo "## Override from .preproc_override: $override_string" >> "${logFile}"
fi

#setup qa_images directory
mpragedir=$(pwd)
qa_imgdir="${mpragedir}/qa_images"
[ ! -d "${qa_imgdir}" ] && mkdir "${qa_imgdir}"
qa_imglog="${qa_imgdir}/qa_images.log"

#if .mprage_nifti file exists, then DICOM -> NIFTI conversion completed and we should read nifti from file
if [[ -z "$nifti" && -f .mprage_nifti ]]; then
    read nifti < .mprage_nifti
fi

#check whether a .nii file was passed in, but FSL is using .nii.gz (or vice versa)
#in this case it is crucial to gzip or gunzip the file up front so that fslreorient2std overwrites
#the initial file and the reoriented file is used downstream

#for the moment, I haven't added the same check to preprocessFunctional because it prepends an underscore at the
#fslreorient2std step, which would avoid the ambiguity problem here.
if [ -n "$nifti" ]; then
    if [[ $nifti =~ .*.nii$ && $ext == .nii.gz ]]; then
	rel "Gzipping initial nifti input image to avoid ambiguous files downstream (FSL expects .nii.gz)" c
	rel "gzip $nifti"
    elif [[ $nifti =~ .*.nii.gz$ && $ext == .nii ]]; then
	rel "Gunzipping initial nifti input image to avoid ambiguous files downstream (FSL expects .nii)" c
	rel "gunzip $nifti"
    fi
fi

#if nifti undefined, assume the dicoms need to be converted
if [ -z "$nifti" ]; then
    ##############
    #convert dicom files to NIFTI
    rel "Converting DICOM files to NIFTI" c

    #check whether files exist
    #numFiles=$( ls | grep "$dicomPattern"  | wc -l )
    numFiles=$( ls $dicomPattern | wc -l )
    if [ $numFiles -eq 0 ]; then
	echo "No DICOM files found using pattern: $dicomPattern. If you have already converted DICOM to NIFTI and want to skip this step, pass in the unbetted structural image using the -n parameter. Example: preprocessMprage -n mprage.nii.gz"
	printHelp
	exit 1
    fi

    nifti="mprage.nii.gz" #used in bet step.

    #remove mprage.nii if it exists so that Dimon doesn't bomb out
    if [ -f mprage.nii ]; then
	rel "mv -f mprage.nii \"mprageBAK_$(date $datefmt).nii\""
    fi

    if [ -f mprage.nii.gz ]; then
	rel "mv -f mprage.nii.gz \"mprageBAK_$(date $datefmt).nii.gz\""
    fi

    #convert dicom to nifti using Dimon
    rel "Dimon \
	-infile_pattern \"${dicomPattern}\" \
	-GERT_Reco \
	-quit \
	-dicom_org \
	-sort_by_acq_time \
	-gert_write_as_nifti \
	-gert_create_dataset \
	-gert_to3d_prefix mprage"

    rm -f dimon.files*
    rm -f GERT_Reco_dicom*

    #if afnirc has compressor on, then above will already generate nii.gz
    if [ -f mprage.nii ]; then
	rel "gzip -f mprage.nii" #use -f to force overwrite in case where mprage.nii.gz exists, but we want to replace it.
    fi

    echo "mprage.nii.gz" > .mprage_nifti

    #Ask user what to do with original DICOM files unless passed on command line
    if [ -z $delDicom ]; then
	until [[ "$delDicom" = [AaNnYy] ]]; do
	    read -sn1 -p "Delete or archive original DICOM files? (y/n/a)" delDicom
	done
    fi

    case ${delDicom} in
	y|Y|yes|delete) echo -e "\nDeleting DICOM files"; rel "rm -f ${dicomPattern}" ;;
	n|N|no) echo -e "\nKeeping DICOM files" ;;
	a|A|archive) echo -e "\nArchiving DICOM files (mprage_dicom.tar.gz)"; rel "tar czf mprage_dicom.tar.gz ${dicomPattern}" && rel "rm -f ${dicomPattern}" ;;
    esac

fi

#strip off file extension to allow for suffix
T1=$( remove_ext $nifti )
[ -n "$brainmask" ] && brainmask=$( remove_ext $brainmask )

#need to set default after processing options and dicom import
#to account for reference choice and filename
if [ -z "$outputFile" ]; then outputFile="${T1}_nonlinear_warp_${reference}"; fi

#check for cleanup only
if [ $cleanup_only -eq 1 ]; then
    rel "Cleaning up intermediate files for preprocessMprage." c
    cleanup_preprocessMprage
    exit 0
fi

#### START preprocessing
backup_original "$T1" "$nifti" # get back to original or make a copy so we can later
reorient        # match templte LPI
[ -n "$CUT_ZIDX" ] && cut_zindex $CUT_ZIDX
[ $deneck  -eq 1 ] && deneck
[ $unifize -eq 1 ] && unifize
[ $deface  -eq 1 ] && deface

# downstream uses "$nifti"
last_prepre="$(current_step_file)"
rel "setting '$nifti' to the last pre-preproc step '$last_prepre'" c
rel "3dcopy -overwrite $last_prepre $nifti"

if [ $robustfov -eq 1 ]; then
    if [ $( imtest "${T1}_fullfov" ) -eq 0 ]; then
	rel "Copying original mprage before FOV cleanup to ${T1}_fullfov" c
	rel "imcp ${T1} ${T1}_fullfov" #shouldn't re-run this step on clipped data if fullfov already exists because it will keep chopping 1mm
    fi
    rel "robustfov -i \"${T1}_fullfov\" -r \"${T1}\""

    # 20210630 - DISABLED!
    #  use 3dresmample below to ensure mask matches anat
    if [ 0 -eq 1 -a  -n "$brainmask" ]; then
	[ $( imtest "${brainmask}_fullfov" ) -eq 0 ] && rel "imcp \"${brainmask}\" \"${brainmask}_fullfov\""
	rel "robustfov -i \"${brainmask}_fullfov\" -r \"${brainmask}\""
    fi
fi

if [ -n "$brainmask" ]; then 
   rel "fslreorient2std $brainmask ${brainmask}_std" # flip storage and header
   if 3dinfo -same_all_grid ${brainmask}_std$ext $T1${ext} |grep 0  -q; then
      rel "WARNING: mask doesn't match $T1 in dim/delta (maybe b/c robustfov, see -no_robustfov) resampling!" c
      rel "3dresample -master $T1$ext -inset ${brainmask}_std$ext -overwrite -rmode NN -prefix ${brainmask}_std$ext"
   fi

fi

#initial correction for any negative values in original image
fix_negatives "$T1"

#adapted from fsl_anat to create a blurred image by downsampling
function quick_smooth() {
  in=$1
  out=$2
  rel "fslmaths $in -subsamp2 -subsamp2 -subsamp2 -subsamp2 vol16"
  rel "flirt -in vol16 -ref $in -out $out -noresampblur -applyxfm -paddingsize 16"
  # possibly do a tiny extra smooth to $out here?
  rel "imrm vol16"
}

type=1 # For FAST: 1 = T1w, 2 = T2w, 3 = PD

#Jan2017: bet with -f 0.1 on the hpf smoothed data is giving some insane results. This appears to result from
#  air voxels becoming rather speckled as a function of the division by the smoothed image, resulting in the
#  skull-stripping algorithm breaking down. More sensible in my mind is to create the rough brain mask up front
#  and apply it within the bias correction steps.

#get a rough brain mask - it can be *VERY* rough (i.e. missing huge portions of the brain or including non-brain, but non-background)
#use -f 0.1 to err on being over inclusive
# Nov2016: add -R to the mix to improve initial skull strip that is more sensible (avoids ridiculous over-inclusiveness)
rel "bet ${T1} ${T1}_initial_brain -m -f 0.1 -R"
rel "fslmaths ${T1}_initial_brain_mask -dilF ${T1}_initial_brain_mask_dil1x" #dilate 1x to be inclusive

qa_image "${T1}" "${T1}_initial_brain" mprage_initial_brain.png "Unprocessed mprage overlaid with initial skull strip (bet -f 0.1 -R)"

#Adapted from fsl_anat
#### BIAS FIELD CORRECTION (main work, although also refined later on if segmentation run)
# required input: ${T1}
# output: ${T1}_biascorr  [ other intermediates to be cleaned up ]
if [ $biasCorrect -eq 1 ] ; then
    rel "Performing bias field correction." c
    if [ $strongBias -eq 1 ] ; then
	rel "Estimating and removing field (stage 1 - large-scale fields)" c
	niter=5 #FAST iterations
	smooth=10 # bias field smoothing extent (FWHM) in mm; -l in FAST

	#create a highly smoothed T1 image for getting a rough brain mask
	quick_smooth "${T1}" "${T1}_s20"

	#divide T1 image by smoothed image
	rel "fslmaths ${T1} -div ${T1}_s20 ${T1}_hpf"

        # get a smoothed version without the edge effects
	rel "fslmaths ${T1} -mas ${T1}_initial_brain_mask_dil1x ${T1}_hpf_s20" #mask by original by rough brain
	quick_smooth "${T1}_hpf_s20" "${T1}_hpf_s20" #smooth masked img
	quick_smooth "${T1}_initial_brain_mask_dil1x" "${T1}_initmask_s20"
	rel "fslmaths ${T1}_hpf_s20 -div ${T1}_initmask_s20 -mas ${T1}_initial_brain_mask_dil1x ${T1}_hpf2_s20"
	rel "fslmaths ${T1} -mas ${T1}_initial_brain_mask_dil1x -div ${T1}_hpf2_s20 ${T1}_hpf2_brain"

	# make sure the overall scaling doesn't change (equate medians)
	med0=$( fslstats ${T1} -k ${T1}_initial_brain_mask_dil1x -P 50 )
	med1=$( fslstats ${T1}_hpf2_brain -k ${T1}_initial_brain_mask_dil1x -P 50 )
	rel "fslmaths ${T1}_hpf2_brain -div $med1 -mul $med0 ${T1}_hpf2_brain"

	rel "Estimating and removing bias field (stage 2 - detailed fields)" c
	rel "fast -o ${T1}_initfast -l ${smooth} -b -B -t $type --iter=${niter} --nopve --fixed=0 -v ${T1}_hpf2_brain" #run on T1 with smoothed edge
	rel "fast -o ${T1}_initfast2 -l ${smooth} -b -B -t $type --iter=${niter} --nopve --fixed=0 -v ${T1}_initfast_restore" #run again on bias corrected image from above
    else
	niter=10  #FAST iterations
	smooth=20 #bias field smoothing extent (FWHM) in mm; -l in FAST

	#for weak bias field correction, do not iterate the FAST as above (initfast, initfast2)
	#copy across initial skull stripped brain from bet above
	rel "fslmaths ${T1}_initial_brain ${T1}_initfast2_restore"
    fi

    # Run FAST (again) to improve bias field
    rel "fast -o ${T1}_fast -l ${smooth} -b -B -t $type --iter=${niter} --nopve --fixed=0 -v ${T1}_initfast2_restore"
    
    rel "Extrapolating bias field from central region" c

    # use the latest fast output
    rel "fslmaths ${T1} -div ${T1}_fast_restore -mas ${T1}_initial_brain_mask_dil1x ${T1}_fast_totbias"
    rel "fslmaths ${T1}_initial_brain_mask_dil1x -ero -ero -ero -ero ${T1}_initfast2_brain_mask2"
    rel "fslmaths ${T1}_fast_totbias -sub 1 ${T1}_fast_totbias" #subtract 1 from total bias
    rel "fslsmoothfill -i ${T1}_fast_totbias -m ${T1}_initfast2_brain_mask2 -o ${T1}_fast_bias"
    rel "fslmaths ${T1}_fast_bias -add 1 ${T1}_fast_bias" #re-add 1
    rel "fslmaths ${T1}_fast_totbias -add 1 ${T1}_fast_totbias"
    rel "fslmaths ${T1} -div ${T1}_fast_bias ${T1}_biascorr" #bias-corrected structural
else
    #no bias correction -- just copy T1 as is
    rel "Skipping bias field correction." c
    rel "fslmaths ${T1} ${T1}_biascorr"
fi

#compute a liberal brain mask for the bias-corrected image (with skull) to use after gradient correction.
#this is mostly useful for reducing file sizes during gzip compression. Use the 98_2 method of FEAT
p_2=$( fslstats ${T1}_biascorr -p 2 )
p_98=$( fslstats ${T1}_biascorr -p 98 )
thresh=$( echo "scale=5; $p_2 + ($p_98 - $p_2)/10" | bc )

rel "fslmaths ${T1}_biascorr -thr $thresh -bin -fillh -dilF ${T1}_biascorr_mask"
rel "fslmaths ${T1}_biascorr -mas ${T1}_biascorr_mask ${T1}_biascorr"

command -v runROBEX.sh >/dev/null 2>&1 && init_ss_program=runROBEX.sh || init_ss_program=3dSkullStrip #use ROBEX by default if available

###############
#initial structural brain extraction for FLIRT
rel "Running brain extraction" c
if [ $ssmethod = bet ]; then
    [ -z "$betOpts" ] && betOpts="-R -f 0.5 -v" #default bet options if not set manually
    rel "bet ${T1}_biascorr ${T1}_skullstrip ${betOpts}"
    initstrip="bet"
elif [ -n "${brainmask}" ]; then
    if [ ! -f "${brainmask}${ext}" ]; then
	echo "Cannot locate -custom_brainmask: ${brainmask}${ext}"
	exit 1
    fi
    rel "fslmaths ${T1}_biascorr -mas ${brainmask}_std ${T1}_skullstrip" #apply custom brainmask
    ssmethod=custom #override ssmethod so that fnirt-inv doesn't fire below
    initstrip="${brainmask}"
elif [ $ssmethod = noisy ]; then
    rel "@NoisySkullStrip -input ${T1}_biascorr${ext}"
    rel "3dcopy ${T1}_biascorr${ext}.ns+orig. ${T1}_bet_initial${ext}" #default option for initial strip
    initstrip="@NoisySkullStrip"
elif [ $ssmethod = 3dSkullStrip ]; then
    [ -z "$betOpts" ] && betOpts="-touchup -orig_vol" #default 3dss options
    rel "3dSkullStrip -input ${T1}_biascorr${ext} -prefix ${T1}_skullstrip${ext} ${betOpts}" #default option for initial strip
    initstrip="3dSkullStrip"
elif [[ $ssmethod =~ ^freesurfer ]]; then
   fs_mask=${ssmethod##*=}
   fsbrain_to_local mprage.nii.gz $fs_mask ${T1}_skullstrip${ext}
   initstrip="freesurfer ($fs_mask)"
else
    if [ $init_ss_program = runROBEX.sh ]; then
	rel "runROBEX.sh ${T1}_biascorr${ext} ${T1}_skullstrip${ext}"
	initstrip="ROBEX"
    elif [ $init_ss_program = 3dSkullStrip ]; then
	#TODO: is there a less redundant way of using 3dSkullStrip in this circumstance versus the above?
	[ -z "$betOpts" ] && betOpts="-touchup -orig_vol" #default 3dss options
	rel "3dSkullStrip -input ${T1}_biascorr${ext} -prefix ${T1}_skullstrip${ext} ${betOpts}" #default option for initial strip
	initstrip="3dSkullStrip"
    fi
fi

# 20210510 apply $post_bet_skullmask
# overwrite _skullstrip (with backup to _skullstrip_pre-skullmask). 
# skillstrip used to create ${flirtAffCoef}  (i.e. this mask is also applied when using flirt-inv method)
if [ -n "$post_bet_skullmask" ]; then
   skullstripped=${T1}_skullstrip${ext}
   rel "3dcopy -overwrite $skullstripped ${T1}_skullstrip_pre-skullmask.nii.gz"
   qa_image "${T1}_biascorr" "${T1}_skullstrip_pre-skullmask" mprage_pre-skullstrip.png "Bias-corrected mprage overlaid with skull strip before masking: $initstrip"
   post_bet_skullmask="$(resample_or_keep $skullstripped $post_bet_skullmask)"
   rel "3dcalc -m '$post_bet_skullmask' -i $skullstripped -expr 'not(m)*i' -overwrite -prefix '$skullstripped'"
fi

qa_image "${T1}_biascorr" "${T1}_skullstrip" mprage_skullstrip.png "Bias-corrected mprage overlaid with skull strip: $initstrip"

#Handle gradient nonlinearity correction
#NB: relative to HCP pipeline, this step is executed after initial bias field correction.
#This was chosen because the bias field steps (based on fsl_anat) tended to be prone to small values outside the brain,
#leading the bias-corrected image to have rather low mean values and the skull strip step to under extract.
#There should be no adverse consequences conceptually to putting gradient unwarping just before the structural -> template
#warp since it will simply apply the gradient correction on the full image sans bias field.

#NB: also differing from HCP, the structural pipeline is generally performed with gradient-distorted images.
#Gradient undistortion is applied when computing the warp to MNI to maximize alignment to template.
#It is also applied at the end of the pipeline in files labeled _postgdc.
#The logic is that preprocessFunctional can take advantage of gradient unwarping with no major overhaul if we
# a) pass in the gradient-distorted mprage_bet
# b) pass in the concatenated structural -> undistorted structural -> MNI transformation

warpsuffix= #whether to include _postgdc in the flirt and fnirt steps (if gradient undistortion requested)
if [ -n "$gcoeffs" ]; then
    #check that command exists
    command -v gradient_unwarp.py >/dev/null 2>&1 || { echo "Cannot find gradient_unwarp.py. Aborting." >&2; exit 1; }

    #allow coeffs to be in scripts directory
    if [ ! -r "$gcoeffs" ]; then
	if [ ! -r "${scriptDir}/cfg_files/${gcoeffs}" ]; then
	    echo "Unable to locate -grad_unwarp file: $gcoeffs"
	    exit 1
	else
	    gcoeffs="${scriptDir}/cfg_files/${gcoeffs}"
	fi
    fi

    #initial run of gradient undistortion using python script
    rel "gradient_unwarp.py \"${T1}_biascorr${ext}\" \"${T1}_biascorr_gdc_orig${ext}\" siemens -g \"$gcoeffs\" -n"
    rel "convertwarp --abs --ref=\"${T1}_biascorr\" --warp1=fullWarp_abs --relout --out=gdc_warpfield" #--jacobian=gdc_jacobian
    #fslmaths gdc_jacobian -Tmean gdc_jacobian #omitting because jacobian is unused in the unwarping process    
    warpsuffix=_postgdc
    rel "applywarp --rel --interp=spline -i \"${T1}_biascorr\" -r \"${T1}_biascorr\" -w gdc_warpfield -o \"${T1}_biascorr${warpsuffix}\""
    rel "applywarp --rel --interp=spline -i \"${T1}_skullstrip\" -r \"${T1}_skullstrip\" -w gdc_warpfield -o \"${T1}_skullstrip${warpsuffix}\""

    #correct any negatives from GDC procedure
    fix_negatives "${T1}_biascorr${warpsuffix}"
    fix_negatives "${T1}_skullstrip${warpsuffix}"

    rel "remask after GDC to zero out non-tissue voxels" c
    rel "fslmaths ${T1}_skullstrip -bin -dilF ${T1}_skullstrip_mask_dil1x"
    rel "fslmaths ${T1}_skullstrip${warpsuffix} -mas ${T1}_skullstrip_mask_dil1x ${T1}_skullstrip${warpsuffix}"
    rel "imrm ${T1}_skullstrip_mask_dil1x"
    rel "fslmaths ${T1}_biascorr${warpsuffix} -mas ${T1}_biascorr_mask ${T1}_biascorr${warpsuffix}"
fi


##############
#warp structural to standard space (MNI or Talairach)
#first conduct affine (linear) warping to get linear warp coefficients
#note that flirt works with betted reference brain, but fnirt prefers unbetted
flirtAffCoef="${T1}_to_${reference}_affine.mat"

#use the postgdc files, if available, so that the flirt and fnirt steps are on postgdc files, then we convertwarp to bring together GDC + MNI warps
#if [ $( imtest "${T1}_warp_linear" ) -eq 0 ]; then
rel "Running affine (linear) transformation for structural -> template" c
rel "flirt -in ${T1}_skullstrip${warpsuffix} -ref $bettedRefBrain -omat ${flirtAffCoef} -out ${T1}\"_warp_linear\" -dof 12 -interp spline"
#fi

###############
#now nonlinear warp

rel "Running nonlinear transformation to warp mprage to: ${reference}" c
rel "Using warp resolution of $wr mm" c

#Note: fsl_anat generates a mask from the template brain_mask by -fillh -dilF.
#Previously FSL tools used the dil mask that had been dilated 2x -dilF -dilF.
#Not sure of the relevance (the fsl_anat way leads to a smaller mask).

function bright_skull_correction() {
    [[ -z $bright_skull || $bright_skull -eq 0 ]] && return 0 #do not run function
    
    rel "Applying bright skull intensity correction. Note that this will work best when ${T1}_skullstrip is good." c
    rel "I recommend checking that this step works as expected!" c
    rel "fsleyes bright_skull/${T1}_biascorr${warpsuffix}_origintensities ${T1}_biascorr${warpsuffix}" c
    
    rel "mkdir -p bright_skull"
    p90=$(fslstats ${T1}_biascorr${warpsuffix} -P 90) # 90th percentile for whole head image
    p70brain=$(fslstats ${T1}_skullstrip -P 70) # 70th percentile of brain voxels
    
    #Take the minimum of either 90th percentile of whole head, or 70th percentile of brain voxels
    #Should almost always be usually should be the latter
    #This is our normalization target for bright skull voxels
    
    renormtarget=$( echo $p90 $p70brain | awk '{if ($1 < $2) {print $1} else {print $2}}' )
    rel "fslmaths ${T1}_skullstrip -binv bright_skull/skull_mask" #use skull strip to determine what's inside and outside of brain
    rel "fslmaths ${T1}_skullstrip -bin bright_skull/brain_mask"

    #identify voxels that are outside of brain and exceed 90th percentile threshold (these will be renormalized)
    rel "fslmaths ${T1}_biascorr${warpsuffix} -thr $p90 -mas bright_skull/skull_mask bright_skull/${T1}_skull_highvals"

    #compute mean of high values and renormalization factor. Will shift intensities of high voxels to have mean = renorm
    mhigh=$(fslstats bright_skull/${T1}_skull_highvals -M)
    renorm=$( echo "$mhigh/$renormtarget" | bc -l)

    #divide all non-brain voxels by normalization factor
    rel "fslmaths ${T1}_biascorr${warpsuffix} -div $renorm -mas bright_skull/skull_mask bright_skull/${T1}_skull_renorm"

    #To avoid sharp transitions in intensity correction due to thresholding (where voxels are either in or out of 90th percentile mask),
    #create a smoothed (2mm FWHM) version of the mask, and an inverted (1 - value) smoothed mask.
    #Use these to compute a weighted combination of original and renormalized intensities such that voxels that are completely
    #in the high intensity sections are heavily altered, whereas there is a smooth falloff with mild modifications to adjacent voxels
    
    rel "fslmaths bright_skull/${T1}_skull_highvals -bin -s 1 bright_skull/${T1}_highvals_smooth_mask_s1"
    rel "fslmaths bright_skull/${T1}_highvals_smooth_mask_s1 -mul -1 -add 1 bright_skull/${T1}_highvals_smooth_mask_s1_invert"

    rel "fslmaths ${T1}_biascorr${warpsuffix} -mul bright_skull/${T1}_highvals_smooth_mask_s1_invert -mas bright_skull/skull_mask bright_skull/origweight"
    rel "fslmaths bright_skull/${T1}_skull_renorm -mul bright_skull/${T1}_highvals_smooth_mask_s1 -mas bright_skull/skull_mask bright_skull/renormweight"

    if [ $( echo "$renorm < 1" | bc ) -eq 1 ]; then
	rel "Bright skull correction would actually increase intensities of high intensity voxels." c
	rel "This suggests a problem in the skull strip, the overall pipeline, or this algorithm." c
	rel "Exiting the algorithm without changing intensities of ${T1}_biascorr${warpsuffix}" c
	return 0
    fi
    
    rel "imcp ${T1}_biascorr${warpsuffix} bright_skull/${T1}_biascorr${warpsuffix}_origintensities" #cache for comparison

    #add weighted combination of voxels outside of brain to the untouched brain voxels
    rel "fslmaths ${T1}_biascorr${warpsuffix} -mas bright_skull/brain_mask -add bright_skull/origweight -add bright_skull/renormweight ${T1}_biascorr${warpsuffix}"

    return 0
}

#run bright skull correction, if requested
bright_skull_correction

# 20210709. fnirt input is ${T1}_biascorr${warpsuffix}  (likely 'mprage_biascorr')
#           and is not skullstripped. fnirt template <-> mprage w/skull
#           problem for large FOV baby data

fnirt_in=${T1}_biascorr${warpsuffix}
maybe_dilate_bias_in

#re-run fnirt even if output already exists because files upstream may have been changed (e.g., betted mprage)
rel "fnirt --ref=$unbettedRefBrain --refmask=$refMask --in="$fnirt_in" --aff=$flirtAffCoef \
        --fout=\"$( remove_ext ${outputFile} )_warpfield\" --cout=\"${T1}_warpcoef\" --iout=\"$outputFile\" \
	--config=$fnirtConfig --logout=${T1}_to_${reference}_fnirt_settings.log --warpres=${wr},${wr},${wr} -v"

#compute MNI -> anat transformation
rel "invwarp --ref=${T1}_biascorr${warpsuffix} -w ${T1}_warpcoef -o template_to_subject_warpcoef"

#define warp files used for struct -> MNI and MNI -> struct
mni2anat="template_to_subject_warpcoef"
anat2mni="${T1}_warpcoef"

#Apply gradient correction and nonlinear warp in one interpolation (need to convertwarp to concatenate both)
if [ -n "$gcoeffs" ]; then
    rel "immv \"$outputFile\" \"${outputFile}_twostepinterp\""
    rel "convertwarp --ref=$unbettedRefBrain --warp1=gdc_warpfield --warp2=${T1}_warpcoef --out=${T1}_warpcoef_withgdc --relout"

    #not that ${T1}_biascorr is the pregdc file, so this is a one-step interpolation GDC + MNI
    rel "applywarp --rel --in=${T1}_biascorr --ref=$unbettedRefBrain --warp=${T1}_warpcoef_withgdc --interp=spline --out=\"$outputFile\""
    mni2anat="template_to_subject_warpcoef_gdistort"
    anat2mni="${T1}_warpcoef_withgdc"

    #compute warp of MNI template onto subject's brain, applying gradient disortion
    rel "invwarp --ref=${T1}_biascorr -w gdc_warpfield -o gdc_inv_warpfield"
    rel "convertwarp --ref=${T1}_biascorr --warp1=template_to_subject_warpcoef --warp2=gdc_inv_warpfield -o template_to_subject_warpcoef_gdistort"
fi

#warp template to subject (incl. gradient distortion if needed)
rel "applywarp --interp=spline --in=$bettedRefBrain --ref=${T1}_biascorr -w ${mni2anat} -o template_to_subject_brain"

#compute final skull-stripped image
if [ $ssmethod = fnirt-inv ]; then
    rel "Generating brain mask based on inverse warp from template" c
    rel "applywarp --interp=nn --in=$nodil --ref=${T1}_biascorr -w ${mni2anat} -o ${T1}_biascorr_brain_mask"
    rel "fslmaths ${T1}_biascorr_brain_mask -fillh ${T1}_biascorr_brain_mask"
    rel "fslmaths ${T1}_biascorr -mas ${T1}_biascorr_brain_mask ${T1}_bet"
else
    rel "fslmaths ${T1}_skullstrip ${T1}_bet" #straight copy
    rel "fslmaths ${T1}_skullstrip -bin ${T1}_biascorr_brain_mask -odt char" #create brainmask from non-zero voxels of betted image
fi

qa_image "${T1}_biascorr" "${T1}_bet" mprage_brain.png "Bias-corrected mprage overlaid with final skull strip: $ssmethod"

if [[ -f "${bettedRefBrain}.nii" && ! -h ./template_brain.nii ]]; then
    ln -s "${bettedRefBrain}.nii" ./template_brain.nii
    tbrain=template_brain.nii
elif [[ -f "${bettedRefBrain}.nii.gz" && ! -h ./template_brain.nii.gz ]]; then
    ln -s "${bettedRefBrain}.nii.gz" ./template_brain.nii.gz
    tbrain=template_brain.nii.gz
fi

#Perform tissue segmentation here to avoid running FAST multiple times with parallel preprocessFunctional jobs
#WM segmentation is used for bbr co-registration
#also improves bias field correction here by running again
if [ $( imtest ${T1}_bet_fast ) -eq 0 ]; then
    #adapted from fsl_anat
    # required input: ${T1}_biascorr ${T1}_biascorr_brain ${T1}_biascorr_brain_mask
    # output: ${T1}_biascorr ${T1}_biascorr_brain (modified) ${T1}_fast* (as normally output by fast) ${T1}_fast_bias (modified)
    rel "Performing tissue-type segmentation" c
    
    # setting the smoothing parameter and number of iterations for cases which bias field was not calculated. 
    niter=10 
    smooth=20
	
    rel "fast -o ${T1}_bet_fast -n 3 -g -H 0.1 -l ${smooth} -b -B -t $type --iter=${niter} ${T1}_bet"
    rel "immv ${T1}_biascorr ${T1}_bet_initbias" #bias field correction up to this point
    rel "fslmaths ${T1}_bet_fast_restore ${T1}_bet" #use the further refinement of the bias field as the mprage_bet output

    # extrapolate bias field and apply to the whole head image
    rel "fslmaths ${T1}_bet_initbias -div ${T1}_bet_fast_restore -mas ${T1}_biascorr_brain_mask ${T1}_fast_totbias"
    rel "fslmaths ${T1}_fast_totbias -sub 1 ${T1}_fast_totbias"
    rel "fslsmoothfill -i ${T1}_fast_totbias -m ${T1}_biascorr_brain_mask -o ${T1}_fast_bias"
    rel "fslmaths ${T1}_fast_bias -add 1 ${T1}_fast_bias"
    rel "fslmaths ${T1}_fast_totbias -add 1 ${T1}_fast_totbias"
    rel "fslmaths ${T1}_bet_initbias -div ${T1}_fast_bias ${T1}_biascorr"

    # regenerate the standard space structural with the updated bias field correction applied ${T1}_biascorr
    rel "applywarp --rel --in=${T1}_biascorr --ref=$unbettedRefBrain --warp=${anat2mni} --interp=spline --out=\"$outputFile\""
    rel "fslmaths ${T1}_bet_fast_pve_2 -thr 0.5 -bin ${T1}_bet_fast_wmseg" #create binary WM mask

    #generate skull-stripped version of final warped file (after final bias correction)
    rel "fslmaths ${outputFile} -mas $nodil $( remove_ext ${outputFile} )_bet"
fi

#reapply gradient distortion correction on final native-space files
if [ -n "$gcoeffs" ]; then
    rel "applywarp --rel --interp=spline -i \"${T1}_biascorr\" -r \"${T1}_biascorr\" -w gdc_warpfield -o \"${T1}_biascorr${warpsuffix}\""
    rel "applywarp --rel --interp=spline -i \"${T1}_bet\" -r \"${T1}_bet\" -w gdc_warpfield -o \"${T1}_bet${warpsuffix}\""

    #correct any negatives from GDC procedure
    fix_negatives "${T1}_biascorr${warpsuffix}"
    fix_negatives "${T1}_bet${warpsuffix}"

    #GDC procedure tends to create tiny values at many non-brain voxels since the gradient correction is applied wrt scanner coordinates
    #To reduce gray voxels outside the head and reduce file size (since a run of zeros can be compressed efficiently,
    #whereas a run of small values cannot), mask the postgdc files with dilated pre-gdc files
    rel "fslmaths \"${T1}_bet\" -bin -dilF -fillh bet_mask"
    
    rel "fslmaths \"${T1}_biascorr${warpsuffix}\" -mas \"${T1}_biascorr_mask\" \"${T1}_biascorr${warpsuffix}\""
    rel "fslmaths \"${T1}_bet${warpsuffix}\" -mas bet_mask \"${T1}_bet${warpsuffix}\""
    rel "imrm bet_mask"
fi

#create png to check registration
rel "applywarp --ref=\"$unbettedRefBrain\" --mask=\"$refMask\" --in=\"${T1}_bet\" --warp=\"${anat2mni}\" --interp=spline --out=subject_to_template_brain"

#create png of subject brain overlaid with template
qa_image subject_to_template_brain "$bettedRefBrain" highres2standard1.png "Subject brain with template outline overlaid"

#create png of template brain overlaid with subject
qa_image "$bettedRefBrain" subject_to_template_brain highres2standard2.png "Template with subject brain outline overlaid"

#combine into one image
rel "pngappend ${qa_imgdir}/highres2standard1.png - ${qa_imgdir}/highres2standard2.png ${qa_imgdir}/struct_to_template.png"
rel "rm -f highres2standard2.png  highres2standard1.png"
rel "imrm subject_to_template_brain"

#generate symbolic link to skull-stripped mprage file for compatibility with FSL tools
ln -sfn "${T1}_bet${ext}" "${T1}_brain${ext}"

if [ $cleanup -eq 1 ]; then
    rel "Cleaning up preprocessMprage intermediate files." c
    cleanup_preprocessMprage
fi

# write finish time, mv incomplete to complete
echo "# finished $(date $datefmt)" >> ${dotfile}_incomplete
mv ${dotfile}_incomplete ${dotfile}_complete

command -v fsleyes >/dev/null 2>&1 && viewcmd=fsleyes || viewcmd=fslview

#suggestions for QA
[ -r "qa.txt" ] && rm -f qa.txt
cat > qa.txt <<EOL
QA suggestions (from least to most important):
Initial skull strip. Primarily a worry if large parts of brain tissue are missing (e.g., cerebellum).
   $viewcmd ${T1}${ext} ${T1}_skullstrip${ext}

Linear (affine) alignment of skull-stripped subject brain onto template.
This will be negatively affected if skull strip is poor. But you should only expect overall global alignment
   $viewcmd $tbrain ${T1}_warp_linear${ext}

Nonlinear warp of template brain onto subject (inverse warp). Check especially for global shape problems.
Mismatch between template and subject brains indicates a failure of fnirt, possibly inflating the edge
of the brain too much, which then looks like a 'punched in' quality on this check
   $viewcmd ${T1}_bet${ext} template_to_subject_brain${ext}

Nonlinear warp of subject brain onto template. Look for precise alignment of sulci, especially along midline.
Focus especially on the _bet file compared to template, but also check the warped brain without skull strip
to ensure that no areas are inflated excessively beyond the bounds of the template. However, do not worry if the
skull looks oddly deformed on this check, just focus on the brain.
   $viewcmd $tbrain $( remove_ext ${outputFile} )_bet${ext} $( remove_ext ${outputFile} )${ext}
EOL
# vim: set tabstop=7: