Pipeline Documentation

Below are the scripts for the bash source code used in the pipeline.

Note

This page includes the documentation and source code for all of the scripts used.

lib.sh (Bash Library Functions )

Bash library functions that are referenced throughout the use of the pipeline.

#!/usr/bin/env bash
# -*- coding: utf-8 -*-
#
# DESCRIPTION:
#   Bash/shell function library for logging.
#
#
# NOTE:
#   Google shell style guide is used here for consistency. See the
#   style guide here: https://google.github.io/styleguide/shellguide.html
#


#######################################
# Prints message to the command line interface
#   in some arbitrary color.
# Args:
#   msg
#######################################
echo_color(){
  msg='\033[0;'"${@}"'\033[0m'
  echo -e ${msg}
}


#######################################
# Prints message to the command line interface
#   in red.
# Args:
#   msg
#######################################
echo_red(){
  echo_color '31m'"${@}"
}


#######################################
# Prints message to the command line interface
#   in green.
# Args:
#   msg
#######################################
echo_green(){
  echo_color '32m'"${@}"
}


#######################################
# Prints message to the command line interface
#   in blue.
# Args:
#   msg
#######################################
echo_blue(){
  echo_color '36m'"${@}"
}


#######################################
# Prints message to the command line interface
#   in red when an error is intened to be raised.
# Args:
#   msg
#######################################
exit_error(){
  echo_red "${@}"
  exit 1
}


#######################################
# Logs the command to file, and executes (runs) the command.
# Globals:
#   log
#   err
# Args:
#   Command to be logged and performed.
#######################################
run(){
  echo "${@}"
  "${@}" >>${log} 2>>${err}
  if [[ ! ${?} -eq 0 ]]; then
    echo "failed: see log files ${log} ${err} for details"
    exit 1
  fi
  echo "-----------------------"
}


#######################################
# Logs the command to file.
# Globals:
#   log
#   err
# Args:
#   Command to be logged and performed.
#######################################
log(){
  echo "${@}"
  echo "${@}" >>${log} 2>>${err}
  echo "-----------------------"
  echo "-----------------------" >>${log} 2>>${err}
}


#######################################
# Extracts and merges b0s in a dMRI volume.
# Globals:
#   log
#   err
# Required Args:
#   d, dwi: Input DWI file.
#   b, bval: Corresponding bval file.
#   e, bvec: Corresponding bvec file.
#   o, out: Output file name.
# Returns
#   0 if no errors, non-zero on error.
#######################################
extract_b0(){
  # Parse arguments
  while [[ ${#} -gt 0 ]]; do
    case "${1}" in
      -d|--dwi) shift; local dwi=${1} ;;
      -b|--bval) shift; local bval=${1} ;;
      -e|--bvec) shift; local bvec=${1} ;;
      -o|--out) shift; local out=${1} ;;
      -*) echo_red "$(basename ${0}) | extract_b0: Unrecognized option ${1}" >&2; Usage; ;;
      *) break ;;
    esac
    shift
  done

  # Create tmp dir
  local cwd=${PWD}
  local tmp_dir=$(remove_ext ${out})_tmp_${RANDOM}
  run mkdir -p ${tmp_dir}
  run cd ${tmp_dir}

  # Create mif file
  run mrconvert -fslgrad ${bvec} ${bval} ${dwi} dwi.mif

  # Extract b0s
  run dwiextract -bzero dwi.mif dwi.b0s.nii.gz

  # Merge b0s
  run fslmaths dwi.b0s.nii.gz -Tmean ${out}

  # Clean-up
  cd ${cwd}
  rm -rf ${tmp_dir}
}


#######################################
# N4 retrospective bias correction algorithm.
# Globals:
#   log
#   err
# Args:
#   Same arguments as N4BiasFieldCorrection.
# Returns
#   0 if no errors, non-zero on error.
#######################################
N4(){
  N4BiasFieldCorrection "${@}"
}


#######################################
# Creates mask from b0s of a dMRI.
#
# NOTE:
#   * Still a work in progress.
#   * Currently not used at the moment.
#
# Globals:
#   log
#   err
# Args:
#   Same arguments as N4BiasFieldCorrection.
# Returns
#   0 if no errors, non-zero on error.
#######################################
create_mask(){
  # Set defaults
  frac_int=0.25
  bias_correct="false"

  # Parse arguments
  while [[ ${#} -gt 0 ]]; do
    case "${1}" in
      -d|--dwi) shift; local dwi=${1} ;;
      -b|--bval) shift; local bval=${1} ;;
      -e|--bvec) shift; local bvec=${1} ;;
      -o|--outdir) shift; local out=${1} ;;
      --bias-correct) local bias_correct="true" ;;
      -f|-frac|--frac-int) shift; local frac_int=${1} ;;
      -*) echo_red "$(basename ${0}) | extract_b0: Unrecognized option ${1}" >&2; Usage; ;;
      *) break ;;
    esac
    shift
  done

  # Create tmp dir
  local cwd=${PWD}
  local tmp_dir=$(remove_ext ${out})_tmp_${RANDOM}
  run mkdir -p ${tmp_dir}
  run cd ${tmp_dir}

  # Extract b0s
  run extract_b0 --dwi ${dwi} --bval ${bval} --bvec ${bvec} --out b0s.nii.gz
  b0=$(realpath b0s.nii.gz)

  if [[ "${bias_correct}" = true ]]; then
    run bet ${b0} tmp -R -f 0.1 -m
    run N4 -i ${b0}.nii.gz \
    -x tmp_mask.nii.gz \
    -o "[restore.nii.gz,bias.nii.gz]" \
    -c "[50x50x50,0.001]" \
    -s 2 \
    -b "[100,3]" \
    -t "[0.15,0.01,200]"
  else
    run cp ${b0} restore.nii.gz
  fi

  # WORK IN PROGRESS
  # # Copy files with output preifx
  # run bet restore.nii.gz restore_brain -R -f ${frac_int} -m
  # run imcp restore.nii.gz ${out}/hifib0
  # run imcp restore_brain.nii.gz ${out}_hifi_brain.nii.gz
  # run imcp restore_brain_mask.nii.gz ${out}_hifi_brain_mask.nii.gz

  run cd ${cwd}

}


#######################################
# xfm_tck wrapper function for the
# tractography python CLI. CLI options
# are the same as the referenced command
# line tool.
# Globals:
#   log
#   err
# Args:
#   Same arguments as xfm_tck.py
# Returns
#   0 if no errors, non-zero on error.
#######################################
xfm_tck(){
  local scripts_dir=$(echo $(dirname $(realpath ${0})))
  # local cmd=$(realpath ${scripts_dir}/../pkgs/xfm_tck/xfm_tck.py)
  local cmd=$(realpath ${scripts_dir}/pkgs/xfm_tck/xfm_tck.py)
  run ${cmd} "${@}"
}


#######################################
# dwinfo wrapper function for the
# python CLI. CLI options
# are the same as the referenced command
# line tool.
# Globals:
#   log
#   err
# Args:
#   Same arguments as dwinfo.py
# Returns
#   0 if no errors, non-zero on error.
#######################################
dwinfo(){
  local scripts_dir=$(echo $(dirname $(realpath ${0})))
  local cmd=$(realpath ${scripts_dir}/../pkgs/dwinfo/dwinfo.py)
  run ${cmd} "${@}"
}


#######################################
# ``exists`` wrapper function. Checks
# if the input file or directory exists.
# If the input exists (as either a file
# or directory), then 'True' is
# returned/printed to the command line.
# If the input does not exist, then
# 'False' is printed to the command
# line.
#
# Globals:
#   log
#   err
# Args:
#   Input file or directory.
# Returns
#   True if directory/file exists, and
#     False otherwise
#######################################
exists(){
  # Define input
  input="${1}"
  input=$(realpath ${input})

  # Check if input exists as a file or directory
  if [[ -f ${input} ]] || [[ -d ${input} ]]; then
    echo "True"
  else
    echo "False"
  fi
}

dwi_preproc.sh

The main preprocessing script that gathers all the required information and runs each stage of dMR image preprocessing.

#!/usr/bin/env bash
# -*- coding: utf-8 -*-
#
# DESCRIPTION:
#
#
# NOTE:
#   Google shell style guide is used here for consistency. See the
#   style guide here: https://google.github.io/styleguide/shellguide.html
#

# SET SCRIPT GLOBALS
scripts_dir=$(echo $(dirname $(realpath ${0})))

# SOURCE LOGGING FUNCTIONS
. ${scripts_dir}/src/lib.sh
xfm_tck=$(realpath ${scripts_dir}/pkgs/xfm_tck/xfm_tck.py)

#######################################
# Prints usage to the command line interface.
# Args:
#   None
#######################################
Usage(){
  cat << USAGE
  Usage:

      $(basename ${0}) <--options> [--options]

  Performs similar preprocessing steps to that of the dHCP dMRI preprocessing pipeline.
  These preprocessing steps include:

    1. Topup (distortion estimation)
    2. Eddy (eddy current, motion, distortion, and slice-to-volume motion correction)
    3. QC
    4. DTI model fitting
    5. Tractography (Single-Shell CSD)

  Options marked as REPEATABLE may be specified more than once, however all such options
  must be specified the same number of times.

  Lastly, input data is assumed to be named in the BIDS v1.4.1+ convention, with '*_acq-' containing
  the shells (bvalues) of the acquisition. Other attributes in the filename should include:

    * subject ID (sub-<sub_id>_...)
    * run ID (..._run-<run_id>_...)

  Required arguments
    -d, --dwi                       Input 4D dMR/DW image file.
    -b, --bval                      Corresponding bval file.
    -e, --bvec                      Corresponding bvec file.
    -b0, --b0, --sbref              Reverse phase encoded b0 (single-band reference)
    --slspec                        Slice order specification file.
    --acqp                          Acquisition parameter file.
    --data-dir                      Output parent data directory.
    --template                      REPEATABLE: Standard whole-head template for registration and tractography.
    --template-brain                REPEATABLE: Standard brain template for registration and tractography.
    --labels                        REPEATABLE: Corrsponding template labels for tractography.
    --out-tract                     REPEATABLE: Corrsponding output directory basenames for tractography.

  Optional arguments
    --dwi-json                      Corresponding dMR/DW image JSON sidecar.
    --b0-json, --sbref-json         Corresponding b0/sbref JSON sidecar.
    --echo-spacing                  Echo-spacing parameter for the parameter acquisition file [default: 0.05].
    -mb, --multiband-factor         Multiband acceleration factor. NOTE: If this parameter is provided then
                                    '--slspec' does not need to be specified. Additionally, this parameter can
                                    also be specified via a JSON (sidecar) file.
    --idx                           Slice phase encoding index file.
    --mporder                       Number of discrete cosine functions used to model slice-to-volume motion.
                                    Set this parameter to 0 to disable slice-to-volume motion correction and
                                    distortion correction. Otherwise, this parameter is automatically computed.
                                    [default: automatically computed].
    --factor                        Factor to divide the mporder by (if necessary). A factor division of 4
                                    is recommended. [default: 0].
    -h, -help, --help               Prints the help menu, then exits.

USAGE
  exit 1
}


#######################################
# Uses shell's built-in hash function to check dependency.
# Globals:
#   log
#   err
# Args:
#   Dependency to check, e.g. shell command.
# Returns
#   0 if no errors, non-zero on error.
#######################################
dependency_check(){
  if ! hash ${1} 2>/dev/null; then
    exit_error "Dependency ${1} is not installed or added to the system path. Please check. Exiting..."
  fi
}


# SCRIPT BODY

# Set defaults
mporder=""
mb=""
echo_spacing=0.05
factor=0
b0_json=""
slspec=""
acqp=""

# Parse arguments
[[ ${#} -eq 0 ]] && Usage;
while [[ ${#} -gt 0 ]]; do
  case "${1}" in
    -d|--dwi) shift; dwi=${1} ;;
    -b|--bval) shift; bval=${1} ;;
    -e|--bvec) shift; bvec=${1} ;;
    -b0|--b0|--sbref) shift; sbref=${1} ;;
    --factor) shift; factor=${1} ;;
    --slspec) shift; slspec=${1} ;;
    -mb|--multiband-factor) shift; mb=${1} ;;
    --idx) shift; idx=${1} ;;
    --acqp) shift; acqp=${1} ;;
    --dwi-json) shift; dwi_json=${1} ;;
    --b0-json|--sbref-json) shift; b0_json=${1} ;;
    --data-dir) shift; data_dir=${1} ;;
    --template) shift; templates+=( ${1} ) ;;
    --template-brain) shift; template_brains+=( ${1} ) ;;
    --labels) shift; labels+=( ${1} ) ;;
    --out-tract) shift; out_tracts+=( ${1} ) ;;
    --mporder) shift; mporder=${1} ;;
    --echo-spacing) shift; echo_spacing=${1} ;;
    -h|-help|--help) shift; Usage; ;;
    -*) echo_red "$(basename ${0}): Unrecognized option ${1}" >&2; Usage; ;;
    *) break ;;
  esac
  shift
done

# Check args
if [[ ${#templates[@]} -eq ${#template_brains[@]} ]] && [[ ${#template_brains[@]} -eq ${#labels[@]} ]] && [[ ${#labels[@]} -eq ${#out_tracts[@]} ]]; then
  echo ""
else
  echo_red "Unequal number of mulit-input options."
  Usage;
fi

# Check dependencies
deps=( topup eddy mrconvert dwiextract ss3t_csd_beta1 )

for dep in ${deps[@]}; do
  dependency_check ${dep}
done

# variable info
sub_id=$(echo $(remove_ext $(basename ${dwi})) | sed "s@_@ @g" | awk '{print $1}' | sed "s@sub-@@g")
run_id=$(echo $(remove_ext $(basename ${dwi})) | sed "s@_@ @g" | awk '{print $4}' | sed "s@run-@@g")
bshell=$(echo $(remove_ext $(basename ${dwi})) | sed "s@_@ @g" | awk '{print $2}' | sed "s@acq-@@g")

outdir=${data_dir}/sub-${sub_id}/${bshell}/run-${run_id}
topup_dir=${outdir}/topup
eddy_dir=${outdir}/eddy
preproc_dir=${outdir}/preprocessed_data

log_dir=${outdir}/logs
log=${log_dir}/dwi.log
err=${log_dir}/dwi.err

mkdir -p ${log_dir}

# Preprocess data
${scripts_dir}/src/import.sh \
--b0 ${sbref} \
--dwi ${dwi} \
--bval ${bval} \
--bvec ${bvec} \
--data-dir ${data_dir} \
--dwi-json ${dwi_json}
# --slspec ${slspec} \
# --acqp ${acqp} \
# --multiband-factor ${mb} \
# --echo-spacing ${echo_spacing} \
# --b0-json ${b0_json}

${scripts_dir}/src/run_topup.sh \
--phase ${outdir}/import/phase \
--acqp ${outdir}/import/dwi.params.acqp \
--out-dir ${outdir}

${scripts_dir}/src/run_eddy.sh \
--dwi ${dwi} \
--bval ${bval} \
--bvec ${bvec} \
--outdir ${outdir} \
--acqp ${outdir}/import/dwi.params.acqp \
--slspec ${outdir}/import/dwi.slice_order \
--topup-dir ${topup_dir} \
--factor ${factor}
# --mporder ${mporder} \

${scripts_dir}/src/postproc.sh \
--dwi ${eddy_dir}/eddy_corrected.nii.gz \
--bval ${bval} \
--bvec ${eddy_dir}/eddy_corrected.eddy_rotated_bvecs \
--dwi-json ${dwi_json} \
--outdir ${outdir} \
--eddy-dir ${eddy_dir} \
--slspec ${outdir}/import/dwi.slice_order \
--idx ${outdir}/import/dwi.idx \
--acqp ${outdir}/import/dwi.params.acqp \
--topup-dir ${topup_dir}

# DISREGARD_CONDITIONS=0
DISREGARD_CONDITIONS=1

for (( i=0; i < ${#templates[@]}; i++)); do

  if [[ ! -d ${outdir}/tractography/${out_tracts[$i]} ]] || [[ ${DISREGARD_CONDITIONS} -eq 1 ]]; then
    # xfm_tck \
    ${xfm_tck} \
    --dwi=${preproc_dir}/dwi.nii.gz \
    --bval=${preproc_dir}/dwi.bval \
    --bvec=${preproc_dir}/dwi.bvec \
    --json=${preproc_dir}/dwi.json \
    --log=${log_dir}/tract.log \
    --template=${templates[$i]} \
    --template-brain=${template_brains[$i]} \
    --labels=${labels[$i]} \
    --out-dir=${outdir}/tractography/${out_tracts[$i]} \
    --frac-int=0.25 \
    --QIT \
    --symmetric \
    --zero-diagonal \
    --FA --MD --AD --RD # --no-cleanup
  fi
done

log "END: dMRI Preprocessing"

Import Data (Stage 0)

Import/copy the necessary data to a storage location with the required read/write access.

#!/usr/bin/env bash
# -*- coding: utf-8 -*-
#
# DESCRIPTION:
#
#
# NOTE:
#   Google shell style guide is used here for consistency. See the
#   style guide here: https://google.github.io/styleguide/shellguide.html
#

# SET SCRIPT GLOBALS
scripts_dir=$(echo $(dirname $(realpath ${0})))

# SOURCE LOGGING FUNCTIONS
. ${scripts_dir}/lib.sh
dwinfo=$(realpath ${scripts_dir}/../pkgs/dwinfo/dwinfo.py)

#######################################
# Prints usage to the command line interface.
# Args:
#   None
#######################################
Usage(){
  cat << USAGE
  Usage:

      $(basename ${0}) <--options> [--options]

  Imports the dMR (diffusion magnetic resonance) image data for preprocessing.

  Required arguments
    -d, --dwi                       Input DWI file.
    -b, --bval                      Corresponding bval file.
    -e, --bvec                      Corresponding bvec file.
    -b0, --b0                       Reversed phase encoded b0 (or single band reference).
    --data-dir                      Output parent data directory.
    --slspec                        Slice order specification file.
    --acqp                          Acquisition parameter file.

  Optional arguments
    --idx                           Slice phase encoding index file.
    -mb, --multiband-factor         Multiband acceleration factor.
    --dwi-json                      DWI/dMRI JSON sidecar.
    --b0-json                       b0 JSON sidecar.
    -h, -help, --help               Prints the help menu, then exits.

USAGE
  exit 1
}


#######################################
# Checks the dimensions of an input DWI file
# and a corresponding reverse phase encoded
# sbref/b0.
# Globals:
#   log
#   err
# Args:
#   dwi: Input DWI file.
#   b0: Reversed phase encoded b0 (or single band reference).
# Returns
#   0 if no errors, non-zero on error.
#######################################
check_dim(){
  # Parse arguments
  while [[ ${#} -gt 0 ]]; do
    case "${1}" in
      --dwi) shift; local dwi=${1} ;;
      --b0) shift; local b0=${1} ;;
      -*) echo_red "$(basename ${0}) | check_dim: Unrecognized option ${1}" >&2; Usage; ;;
      *) break ;;
    esac
    shift
  done

  # Set dimension variables
  local x1=$(fslval ${dwi} dim1)
  local y1=$(fslval ${dwi} dim2)
  local z1=$(fslval ${dwi} dim3)

  local x2=$(fslval ${b0} dim1)
  local y2=$(fslval ${b0} dim2)
  local z2=$(fslval ${b0} dim3)

  # Check dimensions
  [[ ${x1} -ne ${x2} ]] && log "ERROR | check_dim: Input DWI and sbref are of different dimensions - DWI x: ${x1} rpe_sbref x: ${x2}" && exit_error "ERROR | check_dim: Input DWI and sbref are of different dimensions - DWI x: ${x1} rpe_sbref x: ${x2}"
  [[ ${y1} -ne ${y2} ]] && log "ERROR | check_dim: Input DWI and sbref are of different dimensions - DWI y: ${y1} rpe_sbref y: ${y2}" && exit_error "ERROR | check_dim: Input DWI and sbref are of different dimensions - DWI y: ${y1} rpe_sbref y: ${y2}"
  [[ ${z1} -ne ${z2} ]] && log "ERROR | check_dim: Input DWI and sbref are of different dimensions - DWI z: ${z1} rpe_sbref z: ${z2}" && exit_error "ERROR | check_dim: Input DWI and sbref are of different dimensions - DWI z: ${z1} rpe_sbref z: ${z2}"
}


#######################################
# Writes the index file for use with FSL's
# topup, and eddy.
#
# NOTE: This function assumes uniform phase-encoding
#   for the DWI/dMRI. The output index file is just
#   a sequence 1's.
# Required Args:
#   dwi: Input DWI file.
#   out-idx: Slice phase encoding index.
# Returns
#   0 if no errors, non-zero on error.
#######################################
write_idx(){
  # Parse arguments
  while [[ ${#} -gt 0 ]]; do
    case "${1}" in
      --dwi) shift; local dwi=${1} ;;
      --out-idx) shift; local out_idx=${1} ;;
      -*) echo_red "$(basename ${0}) | write_idx: Unrecognized option ${1}" >&2; Usage; ;;
      *) break ;;
    esac
    shift
  done

  # Number of volumes/dynamics
  local nvols=( $( echo $(seq 1 1 $(fslval ${dwi} dim4)) ) )

  log "LOG | write_idx: Writing DWI/dMRI index file. NOTE: this assumes uniform phase encoding in the dMRI."

  for i in ${nvols[@]}; do
    echo "1" >> "${out_idx}"
  done
}


#######################################
# Imports relevant information such as slice
# order, phase encoding index, and acquisition
# parameters.
# Globals:
#   log
#   err
# Required Args:
#   slspec: Slice order specification file.
#   acqp: Acquisition parameter file.
#   out-slspec: Output filename of the slice order specification file.
#   out-acqp: Output filename of the acquisition parameter file.
# Optional Args:
#   idx: Slice phase encoding index.
#   out-idx: Output filename of the slice phase encoding index.
#   dwi: Input DWI file.
# Returns
#   0 if no errors, non-zero on error.
#######################################
import_info(){
  # Set defaults
  local idx=""
  local dwi=""

  # Parse arguments
  while [[ ${#} -gt 0 ]]; do
    case "${1}" in
      --slspec) shift; local slspec=${1} ;;
      --idx) shift; local idx=${1} ;;
      --acqp) shift; local acqp=${1} ;;
      --dwi) shift; local dwi=${1} ;;
      --out-slspec) shift; local out_slspec=${1} ;;
      --out-idx) shift; local out_idx=${1} ;;
      --out-acqp) shift; local out_acqp=${1} ;;
      -*) echo_red "$(basename ${0}) | import_info: Unrecognized option ${1}" >&2; Usage; ;;
      *) break ;;
    esac
    shift
  done

  if [[ ! -f ${idx} ]] && [[ -f ${dwi} ]] && [[ ! -z ${out_idx} ]]; then
    run write_idx --dwi ${dwi} --out-idx ${out_idx}
    local idx=${out_idx}
  else
    run cp ${idx} ${out_idx}
  fi

  # Check input arguments
  [[ -z ${slspec} ]] || [[ ! -f ${slspec} ]] && log "ERROR | import_info: Slice specification order file required." && exit_error "ERROR | import_info: Slice specification order file required."
  [[ -z ${idx} ]] || [[ ! -f ${idx} ]] && log "ERROR | import_info: Slice index file required." && exit_error "ERROR | import_info: Slice index file required."
  [[ -z ${acqp} ]] || [[ ! -f ${acqp} ]] && log "ERROR | import_info: Acquisition parameters file required." && exit_error "ERROR | import_info: Acquisition parameters file required."

  [[ -z ${out_slspec} ]] && log "ERROR | import_info: Slice specification order output filename required." && exit_error "ERROR | import_info: Slice specification order output filename required."
  [[ -z ${out_idx} ]] && log "ERROR | import_info: Slice index output filename required." && exit_error "ERROR | import_info: Slice index output filename required."
  [[ -z ${out_acqp} ]] && log "ERROR | import_info: Acquisition parameters output filename required." && exit_error "ERROR | import_info: Acquisition parameters output filename required."

  # Import info
  run cp ${slspec} ${out_slspec}
  run cp ${acqp} ${out_acqp}
}


# SCRIPT MAIN BODY

# Parse arguments
[[ ${#} -eq 0 ]] && Usage;
while [[ ${#} -gt 0 ]]; do
  case "${1}" in
    -d|--dwi) shift; dwi=${1} ;;
    -b|--bval) shift; bval=${1} ;;
    -e|--bvec) shift; bvec=${1} ;;
    -b0|--b0) shift; b0=${1} ;;
    --slspec) shift; slspec=${1} ;;
    -mb|--multiband-factor) shift; mb=${1} ;;
    --idx) shift; idx=${1} ;;
    --acqp) shift; acqp=${1} ;;
    --dwi-json) shift; dwi_json=${1} ;;
    --b0-json) shift; b0_json=${1} ;;
    --data-dir) shift; data_dir=${1} ;;
    --echo-spacing) shift; echo_spacing=${1} ;;
    -h|-help|--help) shift; Usage; ;;
    -*) echo_red "$(basename ${0}): Unrecognized option ${1}" >&2; Usage; ;;
    *) break ;;
  esac
  shift
done


# variable info
sub_id=$(echo $(remove_ext $(basename ${dwi})) | sed "s@_@ @g" | awk '{print $1}' | sed "s@sub-@@g")
run_id=$(echo $(remove_ext $(basename ${dwi})) | sed "s@_@ @g" | awk '{print $4}' | sed "s@run-@@g")
dwi_PE=$(echo $(remove_ext $(basename ${dwi})) | sed "s@_@ @g" | awk '{print $3}' | sed "s@dir-@@g")
bshell=$(echo $(remove_ext $(basename ${dwi})) | sed "s@_@ @g" | awk '{print $2}' | sed "s@acq-@@g")
outdir=${data_dir}/sub-${sub_id}/${bshell}/run-${run_id}

# Check (required) input arguments
if [[ -z ${data_dir} ]]; then
  log "ERROR | import_data: Output data directory was not specified."
  exit_error "ERROR | import_data: Output data directory was not specified."
fi

if [[ -z ${dwi} ]] || [[ ! -f ${dwi} ]]; then
  log "ERROR | import_data: Input DWI file was not specified or does not exist."
  exit_error "ERROR | import_data: Input DWI file was not specified or does not exist."
else
  dwi=$(realpath ${dwi})
fi

if [[ -z ${bval} ]] || [[ ! -f ${bval} ]]; then
  log "ERROR | import_data: Input bval file was not specified or does not exist."
  exit_error "ERROR | import_data: Input bval file was not specified or does not exist."
else
  bval=$(realpath ${bval})
fi

if [[ -z ${bvec} ]] || [[ ! -f ${bvec} ]]; then
  log "ERROR | import_data: Input bvec file was not specified or does not exist."
  exit_error "ERROR | import_data: Input bvec file was not specified or does not exist."
else
  bvec=$(realpath ${bvec})
fi

if [[ -z ${b0} ]] || [[ ! -f ${b0} ]]; then
  log "ERROR | import_data: Input reversed phase (rPE) b0 file was not specified or does not exist."
  exit_error "ERROR | import_data: Input reversed phase (rPE) b0 file was not specified or does not exist."
else
  b0=$(realpath ${b0})
fi

# Check (optional) arguments
if [[ ! -z ${dwi_json} ]]; then
  if [[ -f ${dwi_json} ]]; then
    dwi_json=$(realpath ${dwi_json})
  else
    log "ERROR | import_data: DWI JSON file specified, but it does not exist."
    exit_error "ERROR | import_data: DWI JSON file specified, but it does not exist."
  fi
fi

if [[ ! -z ${b0_json} ]]; then
  if [[ -f ${b0_json} ]]; then
    b0_json=$(realpath ${b0_json})
  else
    log "ERROR | import_data: Reversed phase (rPE) b0 JSON file specified, but it does not exist."
    exit_error "ERROR | import_data: Reversed phase (rPE) b0 JSON file specified, but it does not exist."
  fi
fi

# Declare global variables
cwd=${PWD}
log_dir=${outdir}/logs
log=${log_dir}/dwi.log
err=${log_dir}/dwi.err

# Import data
if [[ ! -d ${outdir}/import ]]; then

  log "START: Import Data"
  run mkdir -p ${log_dir}
  run mkdir -p ${outdir}/import

  check_dim --dwi ${dwi} --b0 ${b0}

  # Temporary directory (variable)
  tmp_dir=${outdir}/import/tmp_dir_${RANDOM}

  # Check if acqp file was passed
  if [[ -z ${acqp} ]] || [[ ! -f ${acqp} ]]; then
    [[ ! -d ${tmp_dir} ]] && run mkdir -p ${tmp_dir}

    if [[ -z ${echo_spacing} ]]; then
      _echo_spacing=$(${dwinfo} read-bids --bids-nifti=${dwi} --bids-label=EchoSpacing)
      echo_spacing=$(python -c "print(float('${_echo_spacing}')*100)")
    fi

    acqp=${tmp_dir}/params.acqp

    echo "0 1 0 ${echo_spacing}" > ${acqp}
    echo "0 -1 0 ${echo_spacing}" >> ${acqp}
  fi

  # Check if slice specification order file was passed
  if [[ -z ${slspec} ]] || [[ ! -f ${slspec} ]]; then

    [[ ! -d ${tmp_dir} ]] && run mkdir -p ${tmp_dir}
    [[ -z ${mb} ]] && mb=$(${dwinfo} read-bids --bids-nifti=${dwi} --bids-label=MultibandAccelerationFactor)

    slspec=${tmp_dir}/dwi.slice.order.txt

    ${dwinfo} sliceorder --bids-nifti=${dwi} --mb-factor=${mb} --output=${slspec} --interleaved
  fi

  run import_info --out-slspec ${outdir}/import/dwi.slice_order --out-acqp ${outdir}/import/dwi.params.acqp --dwi ${dwi} --out-idx ${outdir}/import/dwi.idx --slspec ${slspec} --idx ${idx} --acqp ${acqp}
  run extract_b0 --dwi ${dwi} --bval ${bval} --bvec ${bvec} --out ${outdir}/import/sbref_pa.nii.gz
  run fslmaths ${b0} -Tmean ${outdir}/import/sbref_ap &

  run imcp ${dwi} ${outdir}/import/dwi &
  # run imcp ${b0} ${outdir}/import/sbref_ap &
  run cp ${bval} ${outdir}/import/dwi.bval
  run cp ${bvec} ${outdir}/import/dwi.bvec

  [[ ! -z ${dwi_json} ]] && run cp ${dwi_json} ${outdir}/import/dwi.json
  [[ ! -z ${b0_json} ]] && run cp ${b0_json} ${outdir}/import/phase.json

  rm -rf ${tmp_dir} &

  wait

  # Merge b0s
  run cd ${outdir}
  run fslmerge -t ${outdir}/import/phase ${outdir}/import/sbref_pa ${outdir}/import/sbref_ap

  # Minor clean-up
  run imrm ${outdir}/import/sbref_pa ${outdir}/import/sbref_ap
  run cd ${cwd}

  log "END: Import Data"

else
  log "Import already completed."
fi

Topup (Stage 1)

Perform stage 1 of the preprocessing which includes FSL’s Topup.

#!/usr/bin/env bash
# -*- coding: utf-8 -*-
#
# DESCRIPTION:
#
#
# NOTE:
#   Google shell style guide is used here for consistency. See the
#   style guide here: https://google.github.io/styleguide/shellguide.html
#

# SET SCRIPT GLOBALS
scripts_dir=$(echo $(dirname $(realpath ${0})))

# SOURCE LOGGING FUNCTIONS
. ${scripts_dir}/lib.sh

#######################################
# Prints usage to the command line interface.
# Arguments:
#   None
#######################################
Usage(){
  cat << USAGE
  Usage:

      $(basename ${0}) <--options> [--options]

  Performs dMR image distortion estimation of 4D reverse phase-encoded b0s.

  Required arguments
    -p, --phase                     4D reverse phase-encoded b0s.
    -a, --acqp                      Acquisition parameter file.
    --out-dir                       Parent output directory.

  Optional arguments
    -c, --config                    Topup configuration file.
    -h, -help, --help               Prints the help menu, then exits.

USAGE
  exit 1
}


# SCRIPT MAIN BODY

# Set defaults
config=${scripts_dir}/../misc/b02b0.cnf

# Parse arguments
[[ ${#} -eq 0 ]] && Usage;
while [[ ${#} -gt 0 ]]; do
  case "${1}" in
    -p|--phase) shift; phase=${1} ;;
    -a|--acqp) shift; acqp=${1} ;;
    -c|--config) shift; config=${1} ;;
    --out-dir) shift; outdir=${1} ;;
    -*) echo_red "$(basename ${0}): Unrecognized option ${1}" >&2; Usage; ;;
    *) break ;;
  esac
  shift
done

# Log variabes
log_dir=${outdir}/logs
log=${log_dir}/dwi.log
err=${log_dir}/dwi.err

# Topup output dir
cwd=${PWD}
topup_dir=${outdir}/topup

if [[ ! -d ${topup_dir} ]]; then
  run mkdir -p ${topup_dir}

  cd ${topup_dir}

  run imcp ${phase} phase && phase=${PWD}/phase
  run cp ${acqp} params.acqp && acqp=params.acqp

  log "START: TOPUP"

  # Run topup
  run topup --imain=${phase} --datain=${acqp} --config=${config} --fout=${topup_dir}/fieldmap --iout=${topup_dir}/topup_b0s --out=${topup_dir}/topup_results -v

  # Run BET on topup output
  run fslmaths ${topup_dir}/topup_b0s -Tmean ${topup_dir}/topup_hifib0
  run bet ${topup_dir}/topup_hifib0 ${topup_dir}/nodif_brain -m -f 0.25 -R

  cd ${cwd}

  log "END: TOPUP"
else
  log "TOPUP already completed."
fi

EDDY (Stage 2)

Perform stage 2 of the preprocessing which includes FSL’s EDDY.

#!/usr/bin/env bash
# -*- coding: utf-8 -*-
#
# DESCRIPTION:
#
#
# NOTE:
#   Google shell style guide is used here for consistency. See the
#   style guide here: https://google.github.io/styleguide/shellguide.html
#

# SET SCRIPT GLOBALS
scripts_dir=$(echo $(dirname $(realpath ${0})))

# SOURCE LOGGING FUNCTIONS
. ${scripts_dir}/lib.sh
dwinfo=$(realpath ${scripts_dir}/../pkgs/dwinfo/dwinfo.py)

#######################################
# Prints usage to the command line interface.
# Arguments:
#   None
#######################################
Usage(){
  cat << USAGE
  Usage:

      $(basename ${0}) <--options> [--options]

  Performs eddy current correction, in addition to motion and distortion correction.
  Additionally, slice-to-volume motion correction is also performed.

  Required arguments
    -d, --dwi                       Input DWI file.
    -b, --bval                      Corresponding bval file.
    -e, --bvec                      Corresponding bvec file.
    --data-dir                      Output parent data directory.
    --slspec                        Slice order specification file.
    --acqp                          Acquisition parameter file.
    --idx                           Slice phase encoding index file.
    --outdir                        Parent output directory.
    --topup-dir                     TOPUP output directory.

  Optional arguments
    --mporder                       Number of discrete cosine functions used to model slice-to-volume motion.
                                    Set this parameter to 0 to disable slice-to-volume motion correction and
                                    distortion correction. Otherwise, this parameter is automatically computed.
                                    [default: automatically computed].
    --factor                        Factor to divide the mporder by (if necessary). A factor division of 4
                                    is recommended. [default: 0].
    -h, -help, --help               Prints the help menu, then exits.

USAGE
  exit 1
}


# SCRIPT MAIN BODY

# Set defaults
mporder=""
factor=0

# Parse arguments
[[ ${#} -eq 0 ]] && Usage;
while [[ ${#} -gt 0 ]]; do
  case "${1}" in
    -d|--dwi) shift; dwi=${1} ;;
    -b|--bval) shift; bval=${1} ;;
    -e|--bvec) shift; bvec=${1} ;;
    --slspec) shift; slspec=${1} ;;
    --factor) shift; factor=${1} ;;
    --idx) shift; idx=${1} ;;
    --acqp) shift; acqp=${1} ;;
    --outdir) shift; outdir=${1} ;;
    --topup-dir) shift; topup_dir=${1} ;;
    --mporder) shift; mporder=${1} ;;
    -h|-help|--help) shift; Usage; ;;
    -*) echo_red "$(basename ${0}): Unrecognized option ${1}" >&2; Usage; ;;
    *) break ;;
  esac
  shift
done

# Log variabes
log_dir=${outdir}/logs
log=${log_dir}/dwi.log
err=${log_dir}/dwi.err

# Eddy output dir
cwd=${PWD}
eddy_dir=${outdir}/eddy

# bsub -q gpu-v100 -gpu "num=1" -M 1000 -W 500 -n 1 -J "bash" -R "span[hosts=1]" -Is bash

# Compute mporder
if [[ -z ${mporder} ]]; then
  mporder=$(${dwinfo} mporder --bids-nifti ${dwi} --slice-order=${slspec} --factor-divide=${factor})
fi

if [[ ! -d ${eddy_dir} ]]; then
  run mkdir -p ${eddy_dir}

  log "START: EDDY"

  # Run eddy
  run eddy_cuda \
  --imain=${dwi} \
  --mask=${topup_dir}/nodif_brain_mask.nii.gz \
  --index=${outdir}/import/dwi.idx \
  --bvals=${bval} \
  --bvecs=${bvec} \
  --acqp=${acqp} \
  --out=${eddy_dir}/eddy_corrected \
  --very_verbose \
  --niter=5 \
  --fwhm=10,5,0,0,0 \
  --nvoxhp=5000 \
  --repol \
  --ol_type=both  \
  --ol_nstd=3 \
  --data_is_shelled \
  --cnr_maps \
  --residuals \
  --dont_mask_output \
  --slspec=${slspec} \
  --s2v_niter=10 \
  --mporder=${mporder} \
  --s2v_interp=trilinear \
  --s2v_lambda=1 \
  --topup=${topup_dir}/topup_results \
  --estimate_move_by_susceptibility \
  --mbs_niter=20 \
  --mbs_ksp=10 \
  --mbs_lambda=10

  dwi=${eddy_dir}/eddy_corrected.nii.gz
  bvec=${eddy_dir}/eddy_corrected.eddy_rotated_bvecs

  # Run BET on eddy output
  run extract_b0 --dwi ${dwi} --bval ${bval} --bvec ${bvec} --out ${eddy_dir}/hifib0.nii.gz
  run bet ${eddy_dir}/hifib0 ${eddy_dir}/nodif_brain -m -f 0.25 -R

  log "END: EDDY"
else
  log "EDDY already completed."
fi

Post-process (Stage 3)

Perform stage 3 of the preprocessing which includes:

• Computing DTI (diffusion tensor imaging) metrics (e.g. FA, MD, RD, etc.)

• QC (quality control, via EDDY QUAD)

#!/usr/bin/env bash
# -*- coding: utf-8 -*-
#
# DESCRIPTION:
#
#
# NOTE:
#   Google shell style guide is used here for consistency. See the
#   style guide here: https://google.github.io/styleguide/shellguide.html
#

# SET SCRIPT GLOBALS
scripts_dir=$(echo $(dirname $(realpath ${0})))

# SOURCE LOGGING FUNCTIONS
. ${scripts_dir}/lib.sh

#######################################
# Prints usage to the command line interface.
# Arguments:
#   None
#######################################
Usage(){
  cat << USAGE
  Usage:

      $(basename ${0}) <--options> [--options]

  Performs post-processing, which includes dMR image quality control, and the diffusion model tensor fit.

  Required arguments
    -d, --dwi                       Input DWI file.
    -b, --bval                      Corresponding bval file.
    -e, --bvec                      Corresponding bvec file.
    --dwi-json                      Corresponding dMRI JSON sidecar.
    --outdir                        Output parent data directory.
    --slspec                        Slice order specification file.
    --acqp                          Acquisition parameter file.
    --idx                           Slice phase encoding index file.
    --topup-dir                     TOPUP output directory.
    --eddy-dir                      EDDY output directory.

  Optional arguments
    -h, -help, --help               Prints the help menu, then exits.

USAGE
  exit 1
}


# SCRIPT MAIN BODY

# Parse arguments
[[ ${#} -eq 0 ]] && Usage;
while [[ ${#} -gt 0 ]]; do
  case "${1}" in
    -d|--dwi) shift; dwi=${1} ;;
    -b|--bval) shift; bval=${1} ;;
    -e|--bvec) shift; bvec=${1} ;;
    --dwi-json) shift; dwi_json=${1} ;;
    --outdir) shift; outdir=${1} ;;
    --eddy-dir) shift; eddy_dir=${1} ;;
    --slspec) shift; slspec=${1} ;;
    --idx) shift; idx=${1} ;;
    --acqp) shift; acqp=${1} ;;
    --topup-dir) shift; topup_dir=${1} ;;
    -h|-help|--help) shift; Usage; ;;
    -*) echo_red "$(basename ${0}): Unrecognized option ${1}" >&2; Usage; ;;
    *) break ;;
  esac
  shift
done

# Log variabes
log_dir=${outdir}/logs
log=${log_dir}/dwi.log
err=${log_dir}/dwi.err

# Post-process output dir
cwd=${PWD}
preproc_dir=${outdir}/preprocessed_data
qc_dir=${outdir}/eddy.qc

log "START: POST-PROCESS"

if [[ ! -d ${preproc_dir} ]]; then
  run mkdir -p ${preproc_dir}

  run cd ${preproc_dir}

  # Remove negative intensity values (caused by spline interpolation
  # during preprocessing).
  run fslmaths ${dwi} -thr 0 dwi
  run cp ${bval} dwi.bval
  run cp ${bvec} dwi.bvec

  [[ -f ${dwi_json} ]] && run cp ${dwi_json} dwi.json

  # Create brain mask
  tmp_dir=${preproc_dir}/tmp_${RANDOM}
  run mkdir -p ${tmp_dir}
  run extract_b0 --dwi ${preproc_dir}/dwi.nii.gz --bval ${preproc_dir}/dwi.bval --bvec ${preproc_dir}/dwi.bvec --out ${tmp_dir}/b0.nii.gz
  run bet ${tmp_dir}/b0.nii.gz ${preproc_dir}/nodif_brain -m -f 0.25 -R
  rm -rf ${tmp_dir}

  # DTIFIT
  [[ ! -d ${preproc_dir}/dtifit ]] && run mkdir -p ${preproc_dir}/dtifit
  run dtifit \
  -k dwi \
  -o ${preproc_dir}/dtifit/data \
  -m ${preproc_dir}/nodif_brain_mask \
  -r dwi.bvec \
  -b dwi.bval \
  --save_tensor
fi

if [[ ! -d ${qc_dir} ]]; then
  run cd ${preproc_dir}
  run eddy_quad ${eddy_dir}/eddy_corrected \
  -idx ${idx} \
  -par ${acqp} \
  -m ${preproc_dir}/nodif_brain \
  -b ${preproc_dir}/dwi.bval \
  -g ${preproc_dir}/dwi.bvec \
  -o ${qc_dir} \
  -f ${topup_dir}/fieldmap \
  -s ${slspec} \
  -v
fi
log "END: POST-PROCESS"