diff --git a/DiffusionPreprocessing/DiffPreprocPipeline.sh b/DiffusionPreprocessing/DiffPreprocPipeline.sh
index 501c89298..a33677069 100755
--- a/DiffusionPreprocessing/DiffPreprocPipeline.sh
+++ b/DiffusionPreprocessing/DiffPreprocPipeline.sh
@@ -91,14 +91,12 @@
 #
 #~ND~END~
 
-
 # --------------------------------------------------------------------------------
 #  Usage Description Function
 # --------------------------------------------------------------------------------
 
-show_usage()
-{
-	cat << EOF
+show_usage() {
+	cat <<EOF
 
 Perform the steps of the HCP Diffusion Preprocessing Pipeline
 
@@ -141,6 +139,12 @@ PARAMETERs are: [ ] = optional; < > = user supplied value
                           output the commands that would be executed instead of
                           actually running them. --printcom=echo is intended to
                           be used for testing purposes
+  [--select-best-b0]
+                          If set selects the best b0 for each phase encoding direction
+                          to pass on to topup rather than the default behaviour of
+                          using equally spaced b0's throughout the scan. The best b0
+                          is identified as the least distorted (i.e., most similar to
+                          the average b0 after registration).
   [--extra-eddy-arg=<value>]
                           Generic single token (no whitespace) argument to pass
                           to the DiffPreprocPipeline_Eddy.sh script and subsequently
@@ -238,8 +242,7 @@ EOF
 #  Support Functions
 # --------------------------------------------------------------------------------
 
-get_options()
-{
+get_options() {
 	local arguments=($@)
 
 	# initialize global output variables
@@ -250,6 +253,7 @@ get_options()
 	unset NegInputImages
 	unset echospacing
 	unset GdCoeffs
+	unset SelectBestB0
 	DWIName="Diffusion"
 	DegreesOfFreedom=${DEFAULT_DEGREES_OF_FREEDOM}
 	b0maxbval=${DEFAULT_B0_MAX_BVAL}
@@ -262,123 +266,127 @@ get_options()
 	local numArgs=${#arguments[@]}
 	local argument
 
-	while [ ${index} -lt ${numArgs} ] ; do
+	while [ ${index} -lt ${numArgs} ]; do
 		argument=${arguments[index]}
 
 		case ${argument} in
-			--help)
-				show_usage
-				exit 0
-				;;
-			--version)
-				version_show $@
-				exit 0
-				;;
-			--path=*)
-				StudyFolder=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--subject=*)
-				Subject=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--PEdir=*)
-				PEdir=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--posData=*)
-				PosInputImages=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--negData=*)
-				NegInputImages=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--echospacing=*)
-				echospacing=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--gdcoeffs=*)
-				GdCoeffs=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--dwiname=*)
-				DWIName=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--dof=*)
-				DegreesOfFreedom=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--b0maxbval=*)
-				b0maxbval=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--printcom=*)
-				runcmd=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--extra-eddy-arg=*)
-				extra_eddy_arg=${argument#*=}
-				extra_eddy_args+=" ${extra_eddy_arg} "
-				index=$(( index + 1 ))
-				;;
-			--combine-data-flag=*)
-				CombineDataFlag=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			*)
-				show_usage
-				echo "ERROR: Unrecognized Option: ${argument}"
-				exit 1
-				;;
+		--help)
+			show_usage
+			exit 0
+			;;
+		--version)
+			version_show "$@"
+			exit 0
+			;;
+		--path=*)
+			StudyFolder=${argument#*=}
+			index=$((index + 1))
+			;;
+		--subject=*)
+			Subject=${argument#*=}
+			index=$((index + 1))
+			;;
+		--PEdir=*)
+			PEdir=${argument#*=}
+			index=$((index + 1))
+			;;
+		--posData=*)
+			PosInputImages=${argument#*=}
+			index=$((index + 1))
+			;;
+		--negData=*)
+			NegInputImages=${argument#*=}
+			index=$((index + 1))
+			;;
+		--echospacing=*)
+			echospacing=${argument#*=}
+			index=$((index + 1))
+			;;
+		--gdcoeffs=*)
+			GdCoeffs=${argument#*=}
+			index=$((index + 1))
+			;;
+		--dwiname=*)
+			DWIName=${argument#*=}
+			index=$((index + 1))
+			;;
+		--dof=*)
+			DegreesOfFreedom=${argument#*=}
+			index=$((index + 1))
+			;;
+		--b0maxbval=*)
+			b0maxbval=${argument#*=}
+			index=$((index + 1))
+			;;
+		--printcom=*)
+			runcmd=${argument#*=}
+			index=$((index + 1))
+			;;
+		--select-best-b0)
+			SelectBestB0="true"
+			index=$((index + 1))
+			;;
+		--extra-eddy-arg=*)
+			extra_eddy_arg=${argument#*=}
+			extra_eddy_args+=" ${extra_eddy_arg} "
+			index=$((index + 1))
+			;;
+		--combine-data-flag=*)
+			CombineDataFlag=${argument#*=}
+			index=$((index + 1))
+			;;
+		*)
+			show_usage
+			echo "ERROR: Unrecognized Option: ${argument}"
+			exit 1
+			;;
 		esac
 	done
 
 	local error_msgs=""
 
 	# check required parameters
-	if [ -z ${StudyFolder} ] ; then
+	if [ -z ${StudyFolder} ]; then
 		error_msgs+="\nERROR: <study-path> not specified"
 	fi
 
-	if [ -z ${Subject} ] ; then
+	if [ -z ${Subject} ]; then
 		error_msgs+="\nERROR: <subject-id> not specified"
 	fi
 
-	if [ -z ${PEdir} ] ; then
+	if [ -z ${PEdir} ]; then
 		error_msgs+="\nERROR: <phase-encoding-dir> not specified"
 	fi
 
-	if [ -z ${PosInputImages} ] ; then
+	if [ -z ${PosInputImages} ]; then
 		error_msgs+="\nERROR: <positive-phase-encoded-data> not specified"
 	fi
 
-	if [ -z ${NegInputImages} ] ; then
+	if [ -z ${NegInputImages} ]; then
 		error_msgs+="\nERROR: <negative-phase-encoded-data> not specified"
 	fi
 
-	if [ -z ${echospacing} ] ; then
+	if [ -z ${echospacing} ]; then
 		error_msgs+="\nERROR: <echo-spacing> not specified"
 	fi
 
-	if [ -z ${GdCoeffs} ] ; then
+	if [ -z ${GdCoeffs} ]; then
 		error_msgs+="\nERROR: <path-to-gradients-coefficients-file> not specified"
 	fi
 
-	if [ -z ${b0maxbval} ] ; then
+	if [ -z ${b0maxbval} ]; then
 		error_msgs+="\nERROR: <b0-max-bval> not specified"
 	fi
 
-	if [ -z ${DWIName} ] ; then
+	if [ -z ${DWIName} ]; then
 		error_msgs+="\nERROR: <DWIName> not specified"
 	fi
 
-	if [ -z ${CombineDataFlag} ] ; then
+	if [ -z ${CombineDataFlag} ]; then
 		error_msgs+="\nERROR: CombineDataFlag not specified"
 	fi
 
-	if [ ! -z "${error_msgs}" ] ; then
+	if [ ! -z "${error_msgs}" ]; then
 		show_usage
 		echo -e ${error_msgs}
 		echo ""
@@ -399,31 +407,58 @@ get_options()
 	echo "   b0maxbval: ${b0maxbval}"
 	echo "   runcmd: ${runcmd}"
 	echo "   CombineDataFlag: ${CombineDataFlag}"
+	if [ ! -z ${SelectBestB0} ]; then
+		echo "   SelectBestB0: ${SelectBestB0}"
+	fi
 	echo "   extra_eddy_args: ${extra_eddy_args}"
 	echo "-- ${g_script_name}: Specified Command-Line Parameters - End --"
+
+	if [ ! -z "${SelectBestB0}" ]; then
+		dont_peas_set=false
+		fwhm_set=false
+		if [ ! -z "${extra_eddy_args}" ]; then
+			for extra_eddy_arg in ${extra_eddy_args}; do
+				if [[ ${extra_eddy_arg} == "--fwhm"* ]]; then
+					fwhm_set=true
+				fi
+				if [[ ${extra_eddy_arg} == "--dont_peas"* ]]; then
+					show_usage
+					log_Err "When using --select-best-b0, post-alignment of shells in eddy is required, "
+					log_Err "as the first b0 could be taken from anywhere within the diffusion data and "
+					log_Err "hence might not be aligned to the first diffusion-weighted image."
+					log_Err_Abort "Remove either the --extra_eddy_args=--dont_peas flag or the --select-best-b0 flag"
+				fi
+			done
+		fi
+		if [ ${fwhm_set} == false ]; then
+			log_Warn "Using --select-best-b0 prepends the best b0 to the start of the file passed into eddy."
+			log_Warn "To ensure eddy succesfully aligns this new first b0 with the actual first volume,"
+			log_Warn "we recommend to increase the FWHM for the first eddy iterations if using --select-best-b0"
+			log_Warn "This can be done by setting the --extra_eddy_args=--fwhm=... flag"
+		fi
+	fi
 }
 
 #
 # Function Description
 #  Validate necessary scripts exist
 #
-validate_scripts()
-{
+validate_scripts() {
 	local error_msgs=""
 
-	if [ ! -e ${HCPPIPEDIR}/DiffusionPreprocessing/DiffPreprocPipeline_PreEddy.sh ] ; then
+	if [ ! -e ${HCPPIPEDIR}/DiffusionPreprocessing/DiffPreprocPipeline_PreEddy.sh ]; then
 		error_msgs+="\nERROR: HCPPIPEDIR/DiffusionPreprocessing/DiffPreprocPipeline_PreEddy.sh not found"
 	fi
 
-	if [ ! -e ${HCPPIPEDIR}/DiffusionPreprocessing/DiffPreprocPipeline_Eddy.sh ] ; then
+	if [ ! -e ${HCPPIPEDIR}/DiffusionPreprocessing/DiffPreprocPipeline_Eddy.sh ]; then
 		error_msgs+="\nERROR: HCPPIPEDIR/DiffusionPreprocessing/DiffPreprocPipeline_Eddy.sh not found"
 	fi
 
-	if [ ! -e ${HCPPIPEDIR}/DiffusionPreprocessing/DiffPreprocPipeline_PostEddy.sh ] ; then
+	if [ ! -e ${HCPPIPEDIR}/DiffusionPreprocessing/DiffPreprocPipeline_PostEddy.sh ]; then
 		error_msgs+="\nERROR: HCPPIPEDIR/DiffusionPreprocessing/DiffPreprocPipeline_PostEddy.sh not found"
 	fi
 
-	if [ ! -z "${error_msgs}" ] ; then
+	if [ ! -z "${error_msgs}" ]; then
 		show_usage
 		echo -e ${error_msgs}
 		echo ""
@@ -435,8 +470,7 @@ validate_scripts()
 # Function Description
 #  Main processing of script
 #
-main()
-{
+main() {
 	# Get Command Line Options
 	get_options "$@"
 
@@ -455,6 +489,9 @@ main()
 	pre_eddy_cmd+=" --echospacing=${echospacing} "
 	pre_eddy_cmd+=" --b0maxbval=${b0maxbval} "
 	pre_eddy_cmd+=" --printcom=${runcmd} "
+	if [ ! -z "${SelectBestB0}" ]; then
+		pre_eddy_cmd+=" --select-best-b0 "
+	fi
 
 	log_Msg "pre_eddy_cmd: ${pre_eddy_cmd}"
 	${pre_eddy_cmd}
@@ -467,8 +504,8 @@ main()
 	eddy_cmd+=" --dwiname=${DWIName} "
 	eddy_cmd+=" --printcom=${runcmd} "
 
-	if [ ! -z "${extra_eddy_args}" ] ; then
-		for extra_eddy_arg in ${extra_eddy_args} ; do
+	if [ ! -z "${extra_eddy_args}" ]; then
+		for extra_eddy_arg in ${extra_eddy_args}; do
 			eddy_cmd+=" --extra-eddy-arg=${extra_eddy_arg} "
 		done
 	fi
@@ -486,6 +523,9 @@ main()
 	post_eddy_cmd+=" --dof=${DegreesOfFreedom} "
 	post_eddy_cmd+=" --combine-data-flag=${CombineDataFlag} "
 	post_eddy_cmd+=" --printcom=${runcmd} "
+	if [ ! -z "${SelectBestB0}" ]; then
+		post_eddy_cmd+=" --select-best-b0 "
+	fi
 
 	log_Msg "post_eddy_cmd: ${post_eddy_cmd}"
 	${post_eddy_cmd}
@@ -507,8 +547,8 @@ g_script_name=$(basename "${0}")
 
 # Allow script to return a Usage statement, before any other output
 if [ "$#" = "0" ]; then
-    show_usage
-    exit 1
+	show_usage
+	exit 1
 fi
 
 # Verify that HCPPIPEDIR Environment variable is set
@@ -518,13 +558,13 @@ if [ -z "${HCPPIPEDIR}" ]; then
 fi
 
 # Load function libraries
-source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@"         # Debugging functions; also sources log.shlib
-source ${HCPPIPEDIR}/global/scripts/opts.shlib                 # Command line option functions
-source ${HCPPIPEDIR}/global/scripts/version.shlib	           # version_ functions
+source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib
+source ${HCPPIPEDIR}/global/scripts/opts.shlib         # Command line option functions
+source ${HCPPIPEDIR}/global/scripts/version.shlib      # version_ functions
 
-opts_ShowVersionIfRequested $@
+opts_ShowVersionIfRequested "$@"
 
-if opts_CheckForHelpRequest $@; then
+if opts_CheckForHelpRequest "$@"; then
 	show_usage
 	exit 0
 fi
@@ -538,5 +578,4 @@ log_Check_Env_Var FSLDIR
 #
 # Invoke the 'main' function to get things started
 #
-main $@
-
+main "$@"
diff --git a/DiffusionPreprocessing/DiffPreprocPipeline_Eddy.sh b/DiffusionPreprocessing/DiffPreprocPipeline_Eddy.sh
index ff7573573..ebf23abcd 100755
--- a/DiffusionPreprocessing/DiffPreprocPipeline_Eddy.sh
+++ b/DiffusionPreprocessing/DiffPreprocPipeline_Eddy.sh
@@ -63,14 +63,12 @@
 #
 #~ND~END~
 
-
 # --------------------------------------------------------------------------------
 #  Usage Description Function
 # --------------------------------------------------------------------------------
 
-show_usage()
-{
-	cat << EOF
+show_usage() {
+	cat <<EOF
 
 Perform the Eddy step of the HCP Diffusion Preprocessing Pipeline
 
@@ -206,8 +204,7 @@ EOF
 #  Support Functions
 # --------------------------------------------------------------------------------
 
-get_options()
-{
+get_options() {
 	local arguments=($@)
 
 	# initialize global output variables
@@ -226,7 +223,7 @@ get_options()
 	resamp_value=""
 	unset ol_nstd_value
 	extra_eddy_args=""
-    no_gpu=""
+	no_gpu=""
 	cuda_version=""
 
 	# parse arguments
@@ -234,131 +231,131 @@ get_options()
 	local numArgs=${#arguments[@]}
 	local argument
 
-	while [ ${index} -lt ${numArgs} ] ; do
+	while [ ${index} -lt ${numArgs} ]; do
 		argument=${arguments[index]}
 
 		case ${argument} in
-			--help)
-				show_usage
-				exit 0
-				;;
-			--version)
-				version_show $@
-				exit 0
-				;;
-			--path=*)
-				StudyFolder=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--subject=*)
-				Subject=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--detailed-outlier-stats=*)
-				DetailedOutlierStats=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--detailed-outlier-stats)
-				DetailedOutlierStats="True"
-				index=$(( index + 1 ))
-				;;
-			--replace-outliers=*)
-				ReplaceOutliers=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--replace-outliers)
-				ReplaceOutliers="True"
-				index=$(( index + 1 ))
-				;;
-			--printcom=*)
-				runcmd=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--dwiname=*)
-				DWIName=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--nvoxhp=*)
-				nvoxhp=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--sep_offs_move=*)
-				sep_offs_move=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--sep_offs_move)
-				sep_offs_move="True"
-				index=$(( index + 1 ))
-				;;
-			--sep-offs-move)
-				sep_offs_move="True"
-				index=$(( index + 1 ))
-				;;
-			--rms=*)
-				rms=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--rms)
-				rms="True"
-				index=$(( index + 1 ))
-				;;
-			--ff=*)
-				ff_val=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--dont_peas)
-				dont_peas="--dont_peas"
-				index=$(( index + 1 ))
-				;;
-			--fwhm=*)
-				fwhm_value=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--resamp=*)
-				resamp_value=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--ol_nstd=*)
-				ol_nstd_value=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--extra-eddy-arg=*)
-				extra_eddy_arg=${argument#*=}
-				extra_eddy_args+=" ${extra_eddy_arg} "
-				index=$(( index + 1 ))
-				;;
-			--no-gpu)
-				no_gpu="TRUE"
-				index=$(( index + 1 ))
-				;;
-			--cuda-version=*)
-				cuda_version=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			*)
-				show_usage
-				echo "ERROR: Unrecognized Option: ${argument}"
-				exit 1
-				;;
+		--help)
+			show_usage
+			exit 0
+			;;
+		--version)
+			version_show "$@"
+			exit 0
+			;;
+		--path=*)
+			StudyFolder=${argument#*=}
+			index=$((index + 1))
+			;;
+		--subject=*)
+			Subject=${argument#*=}
+			index=$((index + 1))
+			;;
+		--detailed-outlier-stats=*)
+			DetailedOutlierStats=${argument#*=}
+			index=$((index + 1))
+			;;
+		--detailed-outlier-stats)
+			DetailedOutlierStats="True"
+			index=$((index + 1))
+			;;
+		--replace-outliers=*)
+			ReplaceOutliers=${argument#*=}
+			index=$((index + 1))
+			;;
+		--replace-outliers)
+			ReplaceOutliers="True"
+			index=$((index + 1))
+			;;
+		--printcom=*)
+			runcmd=${argument#*=}
+			index=$((index + 1))
+			;;
+		--dwiname=*)
+			DWIName=${argument#*=}
+			index=$((index + 1))
+			;;
+		--nvoxhp=*)
+			nvoxhp=${argument#*=}
+			index=$((index + 1))
+			;;
+		--sep_offs_move=*)
+			sep_offs_move=${argument#*=}
+			index=$((index + 1))
+			;;
+		--sep_offs_move)
+			sep_offs_move="True"
+			index=$((index + 1))
+			;;
+		--sep-offs-move)
+			sep_offs_move="True"
+			index=$((index + 1))
+			;;
+		--rms=*)
+			rms=${argument#*=}
+			index=$((index + 1))
+			;;
+		--rms)
+			rms="True"
+			index=$((index + 1))
+			;;
+		--ff=*)
+			ff_val=${argument#*=}
+			index=$((index + 1))
+			;;
+		--dont_peas)
+			dont_peas="--dont_peas"
+			index=$((index + 1))
+			;;
+		--fwhm=*)
+			fwhm_value=${argument#*=}
+			index=$((index + 1))
+			;;
+		--resamp=*)
+			resamp_value=${argument#*=}
+			index=$((index + 1))
+			;;
+		--ol_nstd=*)
+			ol_nstd_value=${argument#*=}
+			index=$((index + 1))
+			;;
+		--extra-eddy-arg=*)
+			extra_eddy_arg=${argument#*=}
+			extra_eddy_args+=" ${extra_eddy_arg} "
+			index=$((index + 1))
+			;;
+		--no-gpu)
+			no_gpu="TRUE"
+			index=$((index + 1))
+			;;
+		--cuda-version=*)
+			cuda_version=${argument#*=}
+			index=$((index + 1))
+			;;
+		*)
+			show_usage
+			echo "ERROR: Unrecognized Option: ${argument}"
+			exit 1
+			;;
 		esac
 	done
 
 	local error_msgs=""
 
 	# check required parameters
-	if [ -z ${StudyFolder} ] ; then
+	if [ -z ${StudyFolder} ]; then
 		error_msgs+="\nERROR: <study-path> not specified"
 	fi
 
-	if [ -z ${Subject} ] ; then
+	if [ -z ${Subject} ]; then
 		error_msgs+="\nERROR: <subject-id> not specified"
 	fi
 
-	if [ -z ${DWIName} ] ; then
+	if [ -z ${DWIName} ]; then
 		error_msgs+="\nERROR: <DWIName> not specified"
 	fi
 
-	if [ ! -z "${error_msgs}" ] ; then
+	if [ ! -z "${error_msgs}" ]; then
 		show_usage
 		echo -e ${error_msgs}
 		echo ""
@@ -391,15 +388,14 @@ get_options()
 # Function Description
 #  Validate necessary scripts exist
 #
-validate_scripts()
-{
+validate_scripts() {
 	local error_msgs=""
 
-	if [ ! -e ${HCPPIPEDIR_dMRI}/run_eddy.sh ] ; then
+	if [ ! -e ${HCPPIPEDIR_dMRI}/run_eddy.sh ]; then
 		error_msgs+="\nERROR: ${HCPPIPEDIR_dMRI}/run_eddy.sh not found"
 	fi
 
-	if [ ! -z "${error_msgs}" ] ; then
+	if [ ! -z "${error_msgs}" ]; then
 		show_usage
 		echo -e ${error_msgs}
 		echo ""
@@ -413,8 +409,7 @@ validate_scripts()
 #
 #  Gets user specified command line options, runs Eddy step of Diffusiong Preprocessing
 #
-main()
-{
+main() {
 	# Get Command Line Options
 	get_options "$@"
 
@@ -425,49 +420,49 @@ main()
 	outdir=${StudyFolder}/${Subject}/${DWIName}
 
 	# Determine stats_option value to pass to run_eddy.sh script
-	if [ "${DetailedOutlierStats}" = "True" ] ; then
+	if [ "${DetailedOutlierStats}" = "True" ]; then
 		stats_option="--wss"
 	else
 		stats_option=""
 	fi
 
 	# Determine replace_outliers_option value to pass to run_eddy.sh script
-	if [ "${ReplaceOutliers}" = "True" ] ; then
+	if [ "${ReplaceOutliers}" = "True" ]; then
 		replace_outliers_option="--repol"
 	else
 		replace_outliers_option=""
 	fi
 
 	# Determine nvoxhp_option value to pass to run_eddy.sh script
-	if [ "${nvoxhp}" != "" ] ; then
+	if [ "${nvoxhp}" != "" ]; then
 		nvoxhp_option="--nvoxhp=${nvoxhp}"
 	else
 		nvoxhp_option=""
 	fi
 
 	# Determine sep_offs_move_option value to pass to run_eddy.sh script
-	if [ "${sep_offs_move}" = "True" ] ; then
+	if [ "${sep_offs_move}" = "True" ]; then
 		sep_offs_move_option="--sep_offs_move"
 	else
 		sep_offs_move_option=""
 	fi
 
 	# Determine rms_option value to pass to run_eddy.sh script
-	if [ "${rms}" = "True" ] ; then
+	if [ "${rms}" = "True" ]; then
 		rms_option="--rms"
 	else
 		rms_option=""
 	fi
 
 	# Determine ff_option value to pass to run_eddy.sh script
-	if [ "${ff_val}" != "" ] ; then
+	if [ "${ff_val}" != "" ]; then
 		ff_option="--ff=${ff_val}"
 	else
 		ff_option=""
 	fi
 
 	# Determine ol_nstd value to pass to run_eddy.sh script
-	if [ -z ${ol_nstd_value} ] ; then
+	if [ -z ${ol_nstd_value} ]; then
 		ol_nstd_value_option=""
 	else
 		ol_nstd_value_option="--ol_nstd=${ol_nstd_value}"
@@ -494,18 +489,18 @@ main()
 
 	run_eddy_cmd+=" -w ${outdir}/eddy "
 
-	if [ ! -z "${dont_peas}" ] ; then
+	if [ ! -z "${dont_peas}" ]; then
 		run_eddy_cmd+=" --dont_peas "
 	fi
 
 	run_eddy_cmd+=" --fwhm=${fwhm_value} "
 
-	if [ ! -z "${resamp_value}" ] ; then
+	if [ ! -z "${resamp_value}" ]; then
 		run_eddy_cmd+=" --resamp=${resamp_value} "
 	fi
 
-	if [ ! -z "${extra_eddy_args}" ] ; then
-		for extra_eddy_arg in ${extra_eddy_args} ; do
+	if [ ! -z "${extra_eddy_args}" ]; then
+		for extra_eddy_arg in ${extra_eddy_args}; do
 			run_eddy_cmd+=" --extra-eddy-arg=${extra_eddy_arg} "
 		done
 	fi
@@ -529,8 +524,8 @@ g_script_name=$(basename "${0}")
 
 # Allow script to return a Usage statement, before any other output
 if [ "$#" = "0" ]; then
-    show_usage
-    exit 1
+	show_usage
+	exit 1
 fi
 
 # Verify that HCPPIPEDIR Environment variable is set
@@ -540,13 +535,13 @@ if [ -z "${HCPPIPEDIR}" ]; then
 fi
 
 # Load function libraries
-source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@"         # Debugging functions; also sources log.shlib
-source ${HCPPIPEDIR}/global/scripts/opts.shlib                 # Command line option functions
-source ${HCPPIPEDIR}/global/scripts/version.shlib	           # version_ functions
+source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib
+source ${HCPPIPEDIR}/global/scripts/opts.shlib         # Command line option functions
+source ${HCPPIPEDIR}/global/scripts/version.shlib      # version_ functions
 
-opts_ShowVersionIfRequested $@
+opts_ShowVersionIfRequested "$@"
 
-if opts_CheckForHelpRequest $@; then
+if opts_CheckForHelpRequest "$@"; then
 	show_usage
 	exit 0
 fi
@@ -563,4 +558,4 @@ HCPPIPEDIR_dMRI=${HCPPIPEDIR}/DiffusionPreprocessing/scripts
 #
 # Invoke the 'main' function to get things started
 #
-main $@
+main "$@"
diff --git a/DiffusionPreprocessing/DiffPreprocPipeline_PostEddy.sh b/DiffusionPreprocessing/DiffPreprocPipeline_PostEddy.sh
index 9c645e608..1abbfe0f2 100755
--- a/DiffusionPreprocessing/DiffPreprocPipeline_PostEddy.sh
+++ b/DiffusionPreprocessing/DiffPreprocPipeline_PostEddy.sh
@@ -63,14 +63,12 @@
 #
 #~ND~END~
 
-
 # --------------------------------------------------------------------------------
 #  Usage Description Function
 # --------------------------------------------------------------------------------
 
-show_usage()
-{
-	cat << EOF
+show_usage() {
+	cat <<EOF
 
 Perform the Post-Eddy steps of the HCP Diffusion Preprocessing Pipeline
 
@@ -95,6 +93,9 @@ PARAMETERs are: [ ] = optional; < > = user supplied value
                           output the commands that would be executed instead of
                           actually running them. --printcom=echo is intended to
                           be used for testing purposes
+  [--select_best_b0]
+                          Whether the --select-best-b0 flag was used in
+                          DiffPreprocPipeline_PreEddy.
   [--combine-data-flag=<value>]
                           Specified value is passed as the CombineDataFlag value
                           for the eddy_postproc.sh script.
@@ -145,14 +146,14 @@ EOF
 #  ${CombineDataFlag}     CombineDataFlag value to pass to the eddy_postproc.sh
 #                         script.
 #
-get_options()
-{
+get_options() {
 	local arguments=($@)
 
 	# initialize global output variables
 	unset StudyFolder
 	unset Subject
 	unset GdCoeffs
+	unset SelectBestB0
 	DWIName="Diffusion"
 	DegreesOfFreedom=${DEFAULT_DEGREES_OF_FREEDOM}
 	runcmd=""
@@ -163,82 +164,86 @@ get_options()
 	local numArgs=${#arguments[@]}
 	local argument
 
-	while [ ${index} -lt ${numArgs} ] ; do
+	while [ ${index} -lt ${numArgs} ]; do
 		argument=${arguments[index]}
 
 		case ${argument} in
-			--help)
-				show_usage
-				exit 0
-				;;
-			--version)
-				version_show $@
-				exit 0
-				;;
-			--path=*)
-				StudyFolder=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--subject=*)
-				Subject=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--gdcoeffs=*)
-				GdCoeffs=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--dof=*)
-				DegreesOfFreedom=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--printcom=*)
-				runcmd=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--dwiname=*)
-				DWIName=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--combine-data-flag=*)
-				CombineDataFlag=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			*)
-				show_usage
-				echo "ERROR: Unrecognized Option: ${argument}"
-				exit 1
-				;;
+		--help)
+			show_usage
+			exit 0
+			;;
+		--version)
+			version_show "$@"
+			exit 0
+			;;
+		--path=*)
+			StudyFolder=${argument#*=}
+			index=$((index + 1))
+			;;
+		--subject=*)
+			Subject=${argument#*=}
+			index=$((index + 1))
+			;;
+		--gdcoeffs=*)
+			GdCoeffs=${argument#*=}
+			index=$((index + 1))
+			;;
+		--dof=*)
+			DegreesOfFreedom=${argument#*=}
+			index=$((index + 1))
+			;;
+		--printcom=*)
+			runcmd=${argument#*=}
+			index=$((index + 1))
+			;;
+		--dwiname=*)
+			DWIName=${argument#*=}
+			index=$((index + 1))
+			;;
+		--combine-data-flag=*)
+			CombineDataFlag=${argument#*=}
+			index=$((index + 1))
+			;;
+		--select-best-b0)
+			SelectBestB0="true"
+			index=$((index + 1))
+			;;
+		*)
+			show_usage
+			echo "ERROR: Unrecognized Option: ${argument}"
+			exit 1
+			;;
 		esac
 	done
 
 	local error_msgs=""
 
 	# check required parameters
-	if [ -z ${StudyFolder} ] ; then
+	if [ -z ${StudyFolder} ]; then
 		error_msgs+="\nERROR: <study-path> not specified"
 	fi
 
-	if [ -z ${Subject} ] ; then
+	if [ -z ${Subject} ]; then
 		error_msgs+="\nERROR: <subject-id> not specified"
 	fi
 
-	if [ -z ${GdCoeffs} ] ; then
+	if [ -z ${GdCoeffs} ]; then
 		error_msgs+="\nERROR: <path-to-gradients-coefficients-file> not specified"
 	fi
 
-	if [ -z ${DWIName} ] ; then
+	if [ -z ${DWIName} ]; then
 		error_msgs+="\nERROR: <DWIName> not specified"
 	fi
 
-	if [ -z ${DegreesOfFreedom} ] ; then
+	if [ -z ${DegreesOfFreedom} ]; then
 		error_msgs+="\nERROR: DegreesOfFreedom not specified"
 	fi
 
-	if [ -z ${CombineDataFlag} ] ; then
+	if [ -z ${CombineDataFlag} ]; then
 		error_msgs+="\nERROR: CombineDataFlag not specified"
 	fi
 
-	if [ ! -z "${error_msgs}" ] ; then
+	if [ ! -z "${error_msgs}" ]; then
 		show_usage
 		echo -e ${error_msgs}
 		echo ""
@@ -254,6 +259,9 @@ get_options()
 	echo "   DegreesOfFreedom: ${DegreesOfFreedom}"
 	echo "   runcmd: ${runcmd}"
 	echo "   CombineDataFlag: ${CombineDataFlag}"
+	if [ -z "${SelectBestB0}" ]; then
+		echo "   SelectBestB0: ${SelectBestB0}"
+	fi
 	echo "-- ${g_script_name}: Specified Command-Line Parameters - End --"
 }
 
@@ -261,19 +269,18 @@ get_options()
 # Function Description
 #  Validate necessary scripts exist
 #
-validate_scripts()
-{
+validate_scripts() {
 	local error_msgs=""
 
-	if [ ! -e ${HCPPIPEDIR_dMRI}/eddy_postproc.sh ] ; then
+	if [ ! -e ${HCPPIPEDIR_dMRI}/eddy_postproc.sh ]; then
 		error_msgs+="\nERROR: ${HCPPIPEDIR_dMRI}/eddy_postproc.sh not found"
 	fi
 
-	if [ ! -e ${HCPPIPEDIR_dMRI}/DiffusionToStructural.sh ] ; then
+	if [ ! -e ${HCPPIPEDIR_dMRI}/DiffusionToStructural.sh ]; then
 		error_msgs+="\nERROR: ${HCPPIPEDIR_dMRI}/DiffusionToStructural.sh not found"
 	fi
 
-	if [ ! -z "${error_msgs}" ] ; then
+	if [ ! -z "${error_msgs}" ]; then
 		show_usage
 		echo -e ${error_msgs}
 		echo ""
@@ -287,8 +294,7 @@ validate_scripts()
 #
 #  Gets user specified command line options, runs Post-Eddy steps of Diffusion Preprocessing
 #
-main()
-{
+main() {
 	# Get Command Line Options
 	get_options "$@"
 
@@ -301,14 +307,18 @@ main()
 
 	# Determine whether Gradient Nonlinearity Distortion coefficients are supplied
 	GdFlag=0
-	if [ ! ${GdCoeffs} = "NONE" ] ; then
+	if [ ! ${GdCoeffs} = "NONE" ]; then
 		log_Msg "Gradient nonlinearity distortion correction coefficients found!"
 		GdFlag=1
 	fi
 
 	log_Msg "Running Eddy PostProcessing"
 	# Note that gradient distortion correction is applied after 'eddy' in the dMRI Pipeline
-	${runcmd} ${HCPPIPEDIR_dMRI}/eddy_postproc.sh ${outdir} ${GdCoeffs} ${CombineDataFlag}
+	select_flag="0"
+	if [ ! -z "${SelectBestB0}" ]; then
+		select_flag="1"
+	fi
+	${runcmd} ${HCPPIPEDIR_dMRI}/eddy_postproc.sh ${outdir} ${GdCoeffs} ${CombineDataFlag} ${select_flag}
 
 	# Establish variables that follow naming conventions
 	T1wFolder="${StudyFolder}/${Subject}/T1w" #Location of T1w images
@@ -319,8 +329,8 @@ main()
 	FreeSurferBrainMask="${T1wFolder}/brainmask_fs"
 	RegOutput="${outdir}"/reg/"Scout2T1w"
 	QAImage="${outdir}"/reg/"T1wMulEPI"
-	DiffRes=`${FSLDIR}/bin/fslval ${outdir}/data/data pixdim1`
-	DiffRes=`printf "%0.2f" ${DiffRes}`
+	DiffRes=$(${FSLDIR}/bin/fslval ${outdir}/data/data pixdim1)
+	DiffRes=$(printf "%0.2f" ${DiffRes})
 
 	log_Msg "Running Diffusion to Structural Registration"
 	${runcmd} ${HCPPIPEDIR_dMRI}/DiffusionToStructural.sh \
@@ -340,7 +350,6 @@ main()
 		--gdflag=${GdFlag} \
 		--diffresol=${DiffRes}
 
-
 	to_location="${outdirT1w}/eddylogs"
 	from_directory="${outdir}/eddy"
 	log_Msg "Copying eddy log files to package location: ${to_location}"
@@ -349,13 +358,13 @@ main()
 	from_files=$(ls ${from_directory}/eddy_unwarped_images.* | grep -v .nii)
 
 	mkdir -p ${to_location}
-	for filename in ${from_files} ; do
+	for filename in ${from_files}; do
 		cp -p ${filename} ${to_location}
 	done
 
-  mkdir -p ${outdirT1w}/QC
-  cp -p ${outdir}/QC/* ${outdirT1w}/QC
-  immv ${outdirT1w}/cnr_maps ${outdirT1w}/QC/cnr_maps
+	mkdir -p ${outdirT1w}/QC
+	cp -p ${outdir}/QC/* ${outdirT1w}/QC
+	immv ${outdirT1w}/cnr_maps ${outdirT1w}/QC/cnr_maps
 
 	log_Msg "Completed!"
 	exit 0
@@ -373,8 +382,8 @@ g_script_name=$(basename "${0}")
 
 # Allow script to return a Usage statement, before any other output
 if [ "$#" = "0" ]; then
-    show_usage
-    exit 1
+	show_usage
+	exit 1
 fi
 
 # Verify that HCPPIPEDIR Environment variable is set
@@ -384,13 +393,13 @@ if [ -z "${HCPPIPEDIR}" ]; then
 fi
 
 # Load function libraries
-source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@"         # Debugging functions; also sources log.shlib
-source ${HCPPIPEDIR}/global/scripts/opts.shlib                 # Command line option functions
-source ${HCPPIPEDIR}/global/scripts/version.shlib	           # version_ functions
+source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib
+source ${HCPPIPEDIR}/global/scripts/opts.shlib         # Command line option functions
+source ${HCPPIPEDIR}/global/scripts/version.shlib      # version_ functions
 
-opts_ShowVersionIfRequested $@
+opts_ShowVersionIfRequested "$@"
 
-if opts_CheckForHelpRequest $@; then
+if opts_CheckForHelpRequest "$@"; then
 	show_usage
 	exit 0
 fi
@@ -400,7 +409,7 @@ ${HCPPIPEDIR}/show_version
 # Verify required environment variables are set and log value
 log_Check_Env_Var HCPPIPEDIR
 log_Check_Env_Var FSLDIR
-log_Check_Env_Var HCPPIPEDIR_Global  # Needed in eddy_postproc.sh and DiffusionToStructural.sh
+log_Check_Env_Var HCPPIPEDIR_Global # Needed in eddy_postproc.sh and DiffusionToStructural.sh
 
 # Set other necessary variables, contingent on HCPPIPEDIR
 HCPPIPEDIR_dMRI=${HCPPIPEDIR}/DiffusionPreprocessing/scripts
@@ -408,4 +417,4 @@ HCPPIPEDIR_dMRI=${HCPPIPEDIR}/DiffusionPreprocessing/scripts
 #
 # Invoke the 'main' function to get things started
 #
-main $@
+main "$@"
diff --git a/DiffusionPreprocessing/DiffPreprocPipeline_PreEddy.sh b/DiffusionPreprocessing/DiffPreprocPipeline_PreEddy.sh
index e753d0a7e..1217be0e4 100755
--- a/DiffusionPreprocessing/DiffPreprocPipeline_PreEddy.sh
+++ b/DiffusionPreprocessing/DiffPreprocPipeline_PreEddy.sh
@@ -62,14 +62,12 @@
 #
 #~ND~END~
 
-
 # --------------------------------------------------------------------------------
 #  Usage Description Function
 # --------------------------------------------------------------------------------
 
-show_usage()
-{
-	cat << EOF
+show_usage() {
+	cat <<EOF
 
 Perform the Pre-Eddy steps of the HCP Diffusion Preprocessing Pipeline
 
@@ -99,6 +97,12 @@ PARAMETERs are: [ ] = optional; < > = user supplied value
   [--b0maxbval=<b0-max-bval>]
                           Volumes with a bvalue smaller than this value will be
                           considered as b0s. Defaults to ${DEFAULT_B0_MAX_BVAL}
+  [--select-best-b0]
+                          If set selects the best b0 for each phase encoding direction
+                          to pass on to topup rather than the default behaviour of
+                          using equally spaced b0's throughout the scan. The best b0
+                          is identified as the least distorted (i.e., most similar to
+                          the average b0 after registration).
   [--printcom=<print-command>]
                           Use the specified <print-command> to echo or otherwise
                           output the commands that would be executed instead of
@@ -146,8 +150,7 @@ EOF
 #  Support Functions
 # --------------------------------------------------------------------------------
 
-get_options()
-{
+get_options() {
 	local arguments=($@)
 
 	# initialize global output variables
@@ -157,6 +160,7 @@ get_options()
 	unset PosInputImages
 	unset NegInputImages
 	unset echospacing
+	unset SelectBestB0
 	DWIName="Diffusion"
 	b0maxbval=${DEFAULT_B0_MAX_BVAL}
 	runcmd=""
@@ -166,98 +170,102 @@ get_options()
 	local numArgs=${#arguments[@]}
 	local argument
 
-	while [ ${index} -lt ${numArgs} ] ; do
+	while [ ${index} -lt ${numArgs} ]; do
 		argument=${arguments[index]}
 
 		case ${argument} in
-			--help)
-				show_usage
-				exit 0
-				;;
-			--version)
-				version_show $@
-				exit 0
-				;;
-			--path=*)
-				StudyFolder=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--subject=*)
-				Subject=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--PEdir=*)
-				PEdir=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--posData=*)
-				PosInputImages=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--negData=*)
-				NegInputImages=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--dwiname=*)
-				DWIName=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--echospacing=*)
-				echospacing=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--b0maxbval=*)
-				b0maxbval=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--printcom=*)
-				runcmd=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			*)
-				show_usage
-				echo "ERROR: Unrecognized Option: ${argument}"
-				exit 1
-				;;
+		--help)
+			show_usage
+			exit 0
+			;;
+		--version)
+			version_show "$@"
+			exit 0
+			;;
+		--path=*)
+			StudyFolder=${argument#*=}
+			index=$((index + 1))
+			;;
+		--subject=*)
+			Subject=${argument#*=}
+			index=$((index + 1))
+			;;
+		--PEdir=*)
+			PEdir=${argument#*=}
+			index=$((index + 1))
+			;;
+		--posData=*)
+			PosInputImages=${argument#*=}
+			index=$((index + 1))
+			;;
+		--negData=*)
+			NegInputImages=${argument#*=}
+			index=$((index + 1))
+			;;
+		--dwiname=*)
+			DWIName=${argument#*=}
+			index=$((index + 1))
+			;;
+		--echospacing=*)
+			echospacing=${argument#*=}
+			index=$((index + 1))
+			;;
+		--b0maxbval=*)
+			b0maxbval=${argument#*=}
+			index=$((index + 1))
+			;;
+		--printcom=*)
+			runcmd=${argument#*=}
+			index=$((index + 1))
+			;;
+		--select-best-b0)
+			SelectBestB0="true"
+			index=$((index + 1))
+			;;
+		*)
+			show_usage
+			echo "ERROR: Unrecognized Option: ${argument}"
+			exit 1
+			;;
 		esac
 	done
 
 	local error_msgs=""
 
 	# check required parameters
-	if [ -z ${StudyFolder} ] ; then
+	if [ -z ${StudyFolder} ]; then
 		error_msgs+="\nERROR: <study-path> not specified"
 	fi
 
-	if [ -z ${Subject} ] ; then
+	if [ -z ${Subject} ]; then
 		error_msgs+="\nERROR: <subject-id> not specified"
 	fi
 
-	if [ -z ${PEdir} ] ; then
+	if [ -z ${PEdir} ]; then
 		error_msgs+="\nERROR: <phase-encoding-dir> not specified"
 	fi
 
-	if [ -z ${PosInputImages} ] ; then
+	if [ -z ${PosInputImages} ]; then
 		error_msgs+="\nERROR: <positive-phase-encoded-data> not specified"
 	fi
 
-	if [ -z ${NegInputImages} ] ; then
+	if [ -z ${NegInputImages} ]; then
 		error_msgs+="\nERROR: <negative-phase-encoded-data> not specified"
 	fi
 
-	if [ -z ${echospacing} ] ; then
+	if [ -z ${echospacing} ]; then
 		error_msgs+="\nERROR: <echo-spacing> not specified"
 	fi
 
-	if [ -z ${b0maxbval} ] ; then
+	if [ -z ${b0maxbval} ]; then
 		error_msgs+="\nERROR: <b0-max-bval> not specified"
 	fi
 
-	if [ -z ${DWIName} ] ; then
+	if [ -z ${DWIName} ]; then
 		error_msgs+="\nERROR: <DWIName> not specified"
 	fi
 
-	if [ ! -z "${error_msgs}" ] ; then
+	if [ ! -z "${error_msgs}" ]; then
 		show_usage
 		echo -e ${error_msgs}
 		echo ""
@@ -275,6 +283,9 @@ get_options()
 	echo "   DWIName: ${DWIName}"
 	echo "   b0maxbval: ${b0maxbval}"
 	echo "   runcmd: ${runcmd}"
+	if [ ! -z "${SelectBestB0}" ]; then
+		echo "   SelectBestB0: ${SelectBestB0}"
+	fi
 	echo "-- ${g_script_name}: Specified Command-Line Parameters - End --"
 }
 
@@ -282,19 +293,20 @@ get_options()
 # Function Description
 #  Validate necessary scripts exist
 #
-validate_scripts()
-{
+validate_scripts() {
 	local error_msgs=""
 
-	if [ ! -e ${HCPPIPEDIR_dMRI}/basic_preproc.sh ] ; then
-		error_msgs+="\nERROR: ${HCPPIPEDIR_dMRI}/basic_preproc.sh not found"
-	fi
+	for extension in norm_intensity sequence best_b0; do
+		if [ ! -e ${HCPPIPEDIR_dMRI}/basic_preproc_${extension}.sh ]; then
+			error_msgs+="\nERROR: ${HCPPIPEDIR_dMRI}/basic_preproc_${extension}.sh not found"
+		fi
+	done
 
-	if [ ! -e ${HCPPIPEDIR_dMRI}/run_topup.sh ] ; then
+	if [ ! -e ${HCPPIPEDIR_dMRI}/run_topup.sh ]; then
 		error_msgs+="\nERROR: ${HCPPIPEDIR_dMRI}/run_topup.sh not found"
 	fi
 
-	if [ ! -z "${error_msgs}" ] ; then
+	if [ ! -z "${error_msgs}" ]; then
 		show_usage
 		echo -e ${error_msgs}
 		echo ""
@@ -306,10 +318,8 @@ validate_scripts()
 # Function Description
 #  find the minimum of two specified numbers
 #
-min()
-{
-	if [ $1 -le $2 ]
-	then
+min() {
+	if [ $1 -le $2 ]; then
 		echo $1
 	else
 		echo $2
@@ -322,11 +332,10 @@ min()
 #
 #  Gets user specified command line options, runs Pre-Eddy steps of Diffusion Preprocessing
 #
-main()
-{
+main() {
 	# Hard-Coded variables for the pipeline
-	MissingFileFlag="EMPTY"  # String used in the input arguments to indicate that a complete series is missing
-	b0dist=45                # Minimum distance in volumes between b0s considered for preprocessing
+	MissingFileFlag="EMPTY" # String used in the input arguments to indicate that a complete series is missing
+	b0dist=45               # Minimum distance in volumes between b0s considered for preprocessing
 
 	# Get Command Line Options
 	get_options "$@"
@@ -339,7 +348,7 @@ main()
 	outdirT1w=${StudyFolder}/${Subject}/T1w/${DWIName}
 
 	# Delete any existing output sub-directories
-	if [ -d ${outdir} ] ; then
+	if [ -d ${outdir} ]; then
 		${runcmd} rm -rf ${outdir}/rawdata
 		${runcmd} rm -rf ${outdir}/topup
 		${runcmd} rm -rf ${outdir}/eddy
@@ -358,7 +367,7 @@ main()
 	${runcmd} mkdir ${outdir}/data
 	${runcmd} mkdir ${outdir}/reg
 
-	if [[ ${PEdir} -ne 1 && ${PEdir} -ne 2 ]] ; then
+	if [[ ${PEdir} -ne 1 && ${PEdir} -ne 2 ]]; then
 		log_Msg "ERROR: Invalid Phase Encoding Directory (PEdir} specified: ${PEdir}"
 		exit 1
 	fi
@@ -370,20 +379,20 @@ main()
 
 	# copy positive raw data
 	log_Msg "Copying positive raw data to working directory"
-	PosInputImages=`echo ${PosInputImages} | sed 's/@/ /g'`
+	PosInputImages=$(echo ${PosInputImages} | sed 's/@/ /g')
 	log_Msg "PosInputImages: ${PosInputImages}"
 
 	Pos_count=1
-	for Image in ${PosInputImages} ; do
-		if [[ ${Image} =~ ^.*EMPTY.*$  ]] ; then
+	for Image in ${PosInputImages}; do
+		if [[ ${Image} =~ ^.*EMPTY.*$ ]]; then
 			Image=EMPTY
 		fi
 
-		if [ ${Image} = ${MissingFileFlag} ] ; then
+		if [ ${Image} = ${MissingFileFlag} ]; then
 			PosVols[${Pos_count}]=0
 		else
-			PosVols[${Pos_count}]=`${FSLDIR}/bin/fslval ${Image} dim4`
-			absname=`${FSLDIR}/bin/imglob ${Image}`
+			PosVols[${Pos_count}]=$(${FSLDIR}/bin/fslval ${Image} dim4)
+			absname=$(${FSLDIR}/bin/imglob ${Image})
 			${runcmd} ${FSLDIR}/bin/imcp ${absname} ${outdir}/rawdata/${basePos}_${Pos_count}
 			${runcmd} cp ${absname}.bval ${outdir}/rawdata/${basePos}_${Pos_count}.bval
 			${runcmd} cp ${absname}.bvec ${outdir}/rawdata/${basePos}_${Pos_count}.bvec
@@ -393,20 +402,20 @@ main()
 
 	# copy negative raw data
 	log_Msg "Copying negative raw data to working directory"
-	NegInputImages=`echo ${NegInputImages} | sed 's/@/ /g'`
+	NegInputImages=$(echo ${NegInputImages} | sed 's/@/ /g')
 	log_Msg "NegInputImages: ${NegInputImages}"
 
 	Neg_count=1
-	for Image in ${NegInputImages} ; do
-		if [[ ${Image} =~ ^.*EMPTY.*$  ]] ; then
+	for Image in ${NegInputImages}; do
+		if [[ ${Image} =~ ^.*EMPTY.*$ ]]; then
 			Image=EMPTY
 		fi
 
-		if [ ${Image} = ${MissingFileFlag} ] ; then
+		if [ ${Image} = ${MissingFileFlag} ]; then
 			NegVols[${Neg_count}]=0
 		else
-			NegVols[${Neg_count}]=`${FSLDIR}/bin/fslval ${Image} dim4`
-			absname=`${FSLDIR}/bin/imglob ${Image}`
+			NegVols[${Neg_count}]=$(${FSLDIR}/bin/fslval ${Image} dim4)
+			absname=$(${FSLDIR}/bin/imglob ${Image})
 			${runcmd} ${FSLDIR}/bin/imcp ${absname} ${outdir}/rawdata/${baseNeg}_${Neg_count}
 			${runcmd} cp ${absname}.bval ${outdir}/rawdata/${baseNeg}_${Neg_count}.bval
 			${runcmd} cp ${absname}.bvec ${outdir}/rawdata/${baseNeg}_${Neg_count}.bvec
@@ -414,8 +423,23 @@ main()
 		Neg_count=$((${Neg_count} + 1))
 	done
 
+	#Compute Total_readout in secs with up to 6 decimal places
+	any=$(ls ${outdir}/rawdata/${basePos}*.nii* | head -n 1)
+	if [ ${PEdir} -eq 1 ]; then #RL/LR phase encoding
+		dimP=$(${FSLDIR}/bin/fslval ${any} dim1)
+	elif [ ${PEdir} -eq 2 ]; then #PA/AP phase encoding
+		dimP=$(${FSLDIR}/bin/fslval ${any} dim2)
+	fi
+	dimPminus1=$(($dimP - 1))
+	#Total_readout=EffectiveEchoSpacing*(ReconMatrixPE-1)
+	# Factors such as in-plane acceleration, phase oversampling, phase resolution, phase field-of-view, and interpolation
+	# must already be accounted for as part of the "EffectiveEchoSpacing"
+	ro_time=$(echo "${echospacing} * ${dimPminus1}" | bc -l)
+	ro_time=$(echo "scale=6; ${ro_time} / 1000" | bc -l) # Convert from ms to sec
+	log_Msg "Total readout time is $ro_time secs"
+
 	# verify positive and negative datasets are provided in pairs
-	if [ ${Pos_count} -ne ${Neg_count} ] ; then
+	if [ ${Pos_count} -ne ${Neg_count} ]; then
 		log_Msg "Wrong number of input datasets! Make sure that you provide pairs of input filenames."
 		log_Msg "If the respective file does not exist, use EMPTY in the input arguments."
 		exit 1
@@ -432,36 +456,40 @@ main()
 	log_Msg "Create two files for each phase encoding direction"
 
 	Paired_flag=0
-	for (( j=1; j<${Pos_count}; j++ )) ; do
-		CorrVols=`min ${NegVols[${j}]} ${PosVols[${j}]}`
-		${runcmd} echo ${CorrVols} ${PosVols[${j}]} >> ${outdir}/eddy/Pos_SeriesVolNum.txt
-		if [ ${PosVols[${j}]} -ne 0 ]
-		then
-			${runcmd} echo ${CorrVols} >> ${outdir}/rawdata/${basePos}_SeriesCorrespVolNum.txt
-			if [ ${CorrVols} -ne 0 ]
-			then
+	for ((j = 1; j < ${Pos_count}; j++)); do
+		CorrVols=$(min ${NegVols[${j}]} ${PosVols[${j}]})
+		${runcmd} echo ${CorrVols} ${PosVols[${j}]} >>${outdir}/eddy/Pos_SeriesVolNum.txt
+		if [ ${PosVols[${j}]} -ne 0 ]; then
+			${runcmd} echo ${CorrVols} >>${outdir}/rawdata/${basePos}_SeriesCorrespVolNum.txt
+			if [ ${CorrVols} -ne 0 ]; then
 				Paired_flag=1
 			fi
 		fi
 	done
 
-	for (( j=1; j<${Neg_count}; j++ )) ; do
-		CorrVols=`min ${NegVols[${j}]} ${PosVols[${j}]}`
-		${runcmd} echo ${CorrVols} ${NegVols[${j}]} >> ${outdir}/eddy/Neg_SeriesVolNum.txt
-		if [ ${NegVols[${j}]} -ne 0 ]
-		then
-			${runcmd} echo ${CorrVols} >> ${outdir}/rawdata/${baseNeg}_SeriesCorrespVolNum.txt
+	for ((j = 1; j < ${Neg_count}; j++)); do
+		CorrVols=$(min ${NegVols[${j}]} ${PosVols[${j}]})
+		${runcmd} echo ${CorrVols} ${NegVols[${j}]} >>${outdir}/eddy/Neg_SeriesVolNum.txt
+		if [ ${NegVols[${j}]} -ne 0 ]; then
+			${runcmd} echo ${CorrVols} >>${outdir}/rawdata/${baseNeg}_SeriesCorrespVolNum.txt
 		fi
 	done
 
-	if [ ${Paired_flag} -eq 0 ] ; then
-		log_Msg "Wrong Input! No pairs of phase encoding directions have been found!"
-		log_Msg "At least one pair is needed!"
-		exit 1
+	if [ ${Paired_flag} -eq 0 ]; then
+		log_Err "Wrong Input! No pairs of phase encoding directions have been found!"
+		log_Err_Abort "At least one pair is needed!"
 	fi
 
-	log_Msg "Running Basic Preprocessing"
-	${runcmd} ${HCPPIPEDIR_dMRI}/basic_preproc.sh ${outdir} ${echospacing} ${PEdir} ${b0dist} ${b0maxbval}
+	log_Msg "Running Intensity Normalisation"
+	${runcmd} ${HCPPIPEDIR_dMRI}/basic_preproc_norm_intensity.sh ${outdir} ${b0maxbval}
+
+	if [ ! -z "${SelectBestB0}" ]; then
+		log_Msg "Running basic preprocessing in preparation of topup (using least distorted b0's)"
+		${runcmd} ${HCPPIPEDIR_dMRI}/basic_preproc_best_b0.sh ${outdir} ${ro_time} ${PEdir} ${b0maxbval}
+	else
+		log_Msg "Running basic preprocessing in preparation of topup (using uniformly interspaced b0)"
+		${runcmd} ${HCPPIPEDIR_dMRI}/basic_preproc_sequence.sh ${outdir} ${ro_time} ${PEdir} ${b0dist} ${b0maxbval}
+	fi
 
 	log_Msg "Running Topup"
 	${runcmd} ${HCPPIPEDIR_dMRI}/run_topup.sh ${outdir}/topup
@@ -482,8 +510,8 @@ g_script_name=$(basename "${0}")
 
 # Allow script to return a Usage statement, before any other output
 if [ "$#" = "0" ]; then
-    show_usage
-    exit 1
+	show_usage
+	exit 1
 fi
 
 # Verify that HCPPIPEDIR Environment variable is set
@@ -493,13 +521,13 @@ if [ -z "${HCPPIPEDIR}" ]; then
 fi
 
 # Load function libraries
-source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@"         # Debugging functions; also sources log.shlib
-source ${HCPPIPEDIR}/global/scripts/opts.shlib                 # Command line option functions
-source ${HCPPIPEDIR}/global/scripts/version.shlib	           # version_ functions
+source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib
+source ${HCPPIPEDIR}/global/scripts/opts.shlib         # Command line option functions
+source ${HCPPIPEDIR}/global/scripts/version.shlib      # version_ functions
 
-opts_ShowVersionIfRequested $@
+opts_ShowVersionIfRequested "$@"
 
-if opts_CheckForHelpRequest $@; then
+if opts_CheckForHelpRequest "$@"; then
 	show_usage
 	exit 0
 fi
@@ -509,7 +537,7 @@ ${HCPPIPEDIR}/show_version
 # Verify required environment variables are set and log value
 log_Check_Env_Var HCPPIPEDIR
 log_Check_Env_Var FSLDIR
-log_Check_Env_Var HCPPIPEDIR_Config  # Needed in run_topup.sh
+log_Check_Env_Var HCPPIPEDIR_Config # Needed in run_topup.sh
 
 # Set other necessary variables, contingent on HCPPIPEDIR
 HCPPIPEDIR_dMRI=${HCPPIPEDIR}/DiffusionPreprocessing/scripts
@@ -517,4 +545,4 @@ HCPPIPEDIR_dMRI=${HCPPIPEDIR}/DiffusionPreprocessing/scripts
 #
 # Invoke the 'main' function to get things started
 #
-main $@
+main "$@"
diff --git a/DiffusionPreprocessing/scripts/DiffusionToStructural.sh b/DiffusionPreprocessing/scripts/DiffusionToStructural.sh
index 21a9feec9..bc6b55aaa 100755
--- a/DiffusionPreprocessing/scripts/DiffusionToStructural.sh
+++ b/DiffusionPreprocessing/scripts/DiffusionToStructural.sh
@@ -3,53 +3,52 @@
 source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib
 echo -e "\n START: DiffusionToStructural"
 
-
-########################################## SUPPORT FUNCTIONS ########################################## 
+########################################## SUPPORT FUNCTIONS ##########################################
 
 # function for parsing options
 getopt1() {
-    sopt="$1"
-    shift 1
-    for fn in $@ ; do
-	if [ `echo $fn | grep -- "^${sopt}=" | wc -w` -gt 0 ] ; then
-	    echo $fn | sed "s/^${sopt}=//"
-	    return 0
-	fi
-    done
+	sopt="$1"
+	shift 1
+	for fn in "$@"; do
+		if [ $(echo $fn | grep -- "^${sopt}=" | wc -w) -gt 0 ]; then
+			echo $fn | sed "s/^${sopt}=//"
+			return 0
+		fi
+	done
 }
 
 defaultopt() {
-    echo $1
+	echo $1
 }
 
 ################################################## OPTION PARSING #####################################################
 # Input Variables
 
-FreeSurferSubjectFolder=`getopt1 "--t1folder" $@`  # "$1" #${StudyFolder}/${Subject}/T1w
-FreeSurferSubjectID=`getopt1 "--subject" $@`       # "$2" #Subject ID
-WorkingDirectory=`getopt1 "--workingdir" $@`       # "$3" #Path to registration working dir, e.g. ${StudyFolder}/${Subject}/Diffusion/reg
-DataDirectory=`getopt1 "--datadiffdir" $@`         # "$4" #Path to diffusion space diffusion data, e.g. ${StudyFolder}/${Subject}/Diffusion/data
-T1wImage=`getopt1 "--t1" $@`                       # "$5" #T1w_acpc_dc image
-T1wRestoreImage=`getopt1 "--t1restore" $@`         # "$6" #T1w_acpc_dc_restore image
-T1wBrainImage=`getopt1 "--t1restorebrain" $@`      # "$7" #T1w_acpc_dc_restore_brain image
-BiasField=`getopt1 "--biasfield" $@`               # "$8" #Bias_Field_acpc_dc
-InputBrainMask=`getopt1 "--brainmask" $@`          # "$9" #Freesurfer Brain Mask, e.g. brainmask_fs
-GdcorrectionFlag=`getopt1 "--gdflag" $@`           # "$10"#Flag for gradient nonlinearity correction (0/1 for Off/On) 
-DiffRes=`getopt1 "--diffresol" $@`                 # "$11"#Diffusion resolution in mm (assume isotropic)
-dof=`getopt1 "--dof" $@`                           # Degrees of freedom for registration to T1w (defaults to 6)
+FreeSurferSubjectFolder=$(getopt1 "--t1folder" "$@")  # "$1" #${StudyFolder}/${Subject}/T1w
+FreeSurferSubjectID=$(getopt1 "--subject" "$@")       # "$2" #Subject ID
+WorkingDirectory=$(getopt1 "--workingdir" "$@")       # "$3" #Path to registration working dir, e.g. ${StudyFolder}/${Subject}/Diffusion/reg
+DataDirectory=$(getopt1 "--datadiffdir" "$@")         # "$4" #Path to diffusion space diffusion data, e.g. ${StudyFolder}/${Subject}/Diffusion/data
+T1wImage=$(getopt1 "--t1" "$@")                       # "$5" #T1w_acpc_dc image
+T1wRestoreImage=$(getopt1 "--t1restore" "$@")         # "$6" #T1w_acpc_dc_restore image
+T1wBrainImage=$(getopt1 "--t1restorebrain" "$@")      # "$7" #T1w_acpc_dc_restore_brain image
+BiasField=$(getopt1 "--biasfield" "$@")               # "$8" #Bias_Field_acpc_dc
+InputBrainMask=$(getopt1 "--brainmask" "$@")          # "$9" #Freesurfer Brain Mask, e.g. brainmask_fs
+GdcorrectionFlag=$(getopt1 "--gdflag" "$@")           # "$10"#Flag for gradient nonlinearity correction (0/1 for Off/On)
+DiffRes=$(getopt1 "--diffresol" "$@")                 # "$11"#Diffusion resolution in mm (assume isotropic)
+dof=$(getopt1 "--dof" "$@")                           # Degrees of freedom for registration to T1w (defaults to 6)
 
 # Output Variables
-T1wOutputDirectory=`getopt1 "--datadiffT1wdir" $@` # "$12" #Path to T1w space diffusion data (for producing output)
-RegOutput=`getopt1 "--regoutput" $@`               # "$13" #Temporary file for sanity checks 
-QAImage=`getopt1 "--QAimage" $@`                   # "$14" #Temporary file for sanity checks 
+T1wOutputDirectory=$(getopt1 "--datadiffT1wdir" "$@") # "$12" #Path to T1w space diffusion data (for producing output)
+RegOutput=$(getopt1 "--regoutput" "$@")               # "$13" #Temporary file for sanity checks
+QAImage=$(getopt1 "--QAimage" "$@")                   # "$14" #Temporary file for sanity checks
 
 # Set default option values
-dof=`defaultopt $dof 6`
+dof=$(defaultopt $dof 6)
 
 # Paths for scripts etc (uses variables defined in SetUpHCPPipeline.sh)
 GlobalScripts=${HCPPIPEDIR_Global}
 
-T1wBrainImageFile=`basename $T1wBrainImage`
+T1wBrainImageFile=$(basename $T1wBrainImage)
 regimg="nodif"
 
 ${FSLDIR}/bin/imcp "$T1wBrainImage" "$WorkingDirectory"/"$T1wBrainImageFile"
@@ -87,54 +86,53 @@ ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/nodif_brain_mask -kernel 3D -dilM "
 
 DilationsNum=6 #Dilated mask for masking the final data and grad_dev
 ${FSLDIR}/bin/imcp "$T1wOutputDirectory"/nodif_brain_mask "$T1wOutputDirectory"/nodif_brain_mask_temp
-for (( j=0; j<${DilationsNum}; j++ ))
-do
-    ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/nodif_brain_mask_temp -kernel 3D -dilM "$T1wOutputDirectory"/nodif_brain_mask_temp
+for ((j = 0; j < ${DilationsNum}; j++)); do
+	${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/nodif_brain_mask_temp -kernel 3D -dilM "$T1wOutputDirectory"/nodif_brain_mask_temp
 done
 
 #Rotate bvecs from diffusion to structural space
 ${GlobalScripts}/Rotate_bvecs.sh "$DataDirectory"/bvecs "$WorkingDirectory"/diff2str.mat "$T1wOutputDirectory"/bvecs
 cp "$DataDirectory"/bvals "$T1wOutputDirectory"/bvals
 
-DilateDistance=`echo "$DiffRes * 4" | bc`  # Extrapolates the diffusion data up to 4 voxels outside of the FOV
+DilateDistance=$(echo "$DiffRes * 4" | bc) # Extrapolates the diffusion data up to 4 voxels outside of the FOV
 
 # Register diffusion data to T1w space. Account for gradient nonlinearities (using one-step resampling)  if requested
 if [ ${GdcorrectionFlag} -eq 1 ]; then
-    echo "Correcting Diffusion data for gradient nonlinearities and registering to structural space"
+	echo "Correcting Diffusion data for gradient nonlinearities and registering to structural space"
 	# We combine the GDC warp with the diff2str affine, and apply to warped/data_warped (i.e., the post-eddy, pre-GDC data)
 	# so that we can have a one-step resampling following 'eddy'
-    ${FSLDIR}/bin/convertwarp --rel --relout --warp1="$DataDirectory"/warped/fullWarp --postmat="$WorkingDirectory"/diff2str.mat --ref="$T1wRestoreImage"_${DiffRes} --out="$WorkingDirectory"/grad_unwarp_diff2str
-    # Dilation outside of the field of view to minimise the effect of the hard field of view edge on the interpolation
-    ${CARET7DIR}/wb_command -volume-dilate $DataDirectory/warped/data_warped.nii.gz $DilateDistance NEAREST $DataDirectory/warped/data_warped_dilated.nii.gz
-    ${FSLDIR}/bin/applywarp --rel -i "$DataDirectory"/warped/data_warped_dilated -r "$T1wRestoreImage"_${DiffRes} -w "$WorkingDirectory"/grad_unwarp_diff2str --interp=spline -o "$T1wOutputDirectory"/data
-    ${FSLDIR}/bin/imrm $DataDirectory/warped/data_warped_dilated
-
-    # Transforms field of view mask to T1-weighted space
-    # (Be sure to use the fov_mask derived prior to application of GDC)
-    ${FSLDIR}/bin/applywarp --rel -i "$DataDirectory"/warped/fov_mask_warped -r "$T1wRestoreImage"_${DiffRes} -w "$WorkingDirectory"/grad_unwarp_diff2str --interp=trilinear -o "$T1wOutputDirectory"/fov_mask
-
-    # Transform CNR maps
-    ${CARET7DIR}/wb_command -volume-dilate $DataDirectory/warped/cnr_maps_warped.nii.gz $DilateDistance NEAREST $DataDirectory/warped/cnr_maps_dilated.nii.gz
-    ${FSLDIR}/bin/applywarp --rel -i "$DataDirectory"/warped/cnr_maps_dilated -r "$T1wRestoreImage"_${DiffRes} -w "$WorkingDirectory"/grad_unwarp_diff2str --interp=spline -o "$T1wOutputDirectory"/cnr_maps
-    ${FSLDIR}/bin/imrm $DataDirectory/warped/cnr_maps_dilated
-
-    # Now register the grad_dev tensor
-    ${FSLDIR}/bin/vecreg -i "$DataDirectory"/grad_dev -o "$T1wOutputDirectory"/grad_dev -r "$T1wRestoreImage"_${DiffRes} -t "$WorkingDirectory"/diff2str.mat --interp=spline
-    ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/grad_dev -mas "$T1wOutputDirectory"/nodif_brain_mask_temp "$T1wOutputDirectory"/grad_dev  #Mask-out values outside the brain 
+	${FSLDIR}/bin/convertwarp --rel --relout --warp1="$DataDirectory"/warped/fullWarp --postmat="$WorkingDirectory"/diff2str.mat --ref="$T1wRestoreImage"_${DiffRes} --out="$WorkingDirectory"/grad_unwarp_diff2str
+	# Dilation outside of the field of view to minimise the effect of the hard field of view edge on the interpolation
+	${CARET7DIR}/wb_command -volume-dilate $DataDirectory/warped/data_warped.nii.gz $DilateDistance NEAREST $DataDirectory/warped/data_warped_dilated.nii.gz
+	${FSLDIR}/bin/applywarp --rel -i "$DataDirectory"/warped/data_warped_dilated -r "$T1wRestoreImage"_${DiffRes} -w "$WorkingDirectory"/grad_unwarp_diff2str --interp=spline -o "$T1wOutputDirectory"/data
+	${FSLDIR}/bin/imrm $DataDirectory/warped/data_warped_dilated
+
+	# Transforms field of view mask to T1-weighted space
+	# (Be sure to use the fov_mask derived prior to application of GDC)
+	${FSLDIR}/bin/applywarp --rel -i "$DataDirectory"/warped/fov_mask_warped -r "$T1wRestoreImage"_${DiffRes} -w "$WorkingDirectory"/grad_unwarp_diff2str --interp=trilinear -o "$T1wOutputDirectory"/fov_mask
+
+	# Transform CNR maps
+	${CARET7DIR}/wb_command -volume-dilate $DataDirectory/warped/cnr_maps_warped.nii.gz $DilateDistance NEAREST $DataDirectory/warped/cnr_maps_dilated.nii.gz
+	${FSLDIR}/bin/applywarp --rel -i "$DataDirectory"/warped/cnr_maps_dilated -r "$T1wRestoreImage"_${DiffRes} -w "$WorkingDirectory"/grad_unwarp_diff2str --interp=spline -o "$T1wOutputDirectory"/cnr_maps
+	${FSLDIR}/bin/imrm $DataDirectory/warped/cnr_maps_dilated
+
+	# Now register the grad_dev tensor
+	${FSLDIR}/bin/vecreg -i "$DataDirectory"/grad_dev -o "$T1wOutputDirectory"/grad_dev -r "$T1wRestoreImage"_${DiffRes} -t "$WorkingDirectory"/diff2str.mat --interp=spline
+	${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/grad_dev -mas "$T1wOutputDirectory"/nodif_brain_mask_temp "$T1wOutputDirectory"/grad_dev #Mask-out values outside the brain
 else
-    # Dilation outside of the field of view to minimise the effect of the hard field of view edge on the interpolation
-    ${CARET7DIR}/wb_command -volume-dilate $DataDirectory/data.nii.gz $DilateDistance NEAREST $DataDirectory/data_dilated.nii.gz
-    # Register diffusion data to T1w space without considering gradient nonlinearities
-    ${FSLDIR}/bin/flirt -in "$DataDirectory"/data_dilated -ref "$T1wRestoreImage"_${DiffRes} -applyxfm -init "$WorkingDirectory"/diff2str.mat -interp spline -out "$T1wOutputDirectory"/data
-    ${FSLDIR}/bin/imrm $DataDirectory/data_dilated
-
-    # Transform CNR maps
-    ${CARET7DIR}/wb_command -volume-dilate $DataDirectory/cnr_maps.nii.gz $DilateDistance NEAREST $DataDirectory/cnr_maps_dilated.nii.gz
-    ${FSLDIR}/bin/flirt -in "$DataDirectory"/cnr_maps_dilated -ref "$T1wRestoreImage"_${DiffRes} -applyxfm -init "$WorkingDirectory"/diff2str.mat -interp spline -out "$T1wOutputDirectory"/cnr_maps
-    ${FSLDIR}/bin/imrm $DataDirectory/cnr_maps_dilated
-
-    # Transforms field of view mask to T1-weighted space
-    ${FSLDIR}/bin/flirt -in "$DataDirectory"/fov_mask -ref "$T1wRestoreImage"_${DiffRes} -applyxfm -init "$WorkingDirectory"/diff2str.mat -interp trilinear -out "$T1wOutputDirectory"/fov_mask
+	# Dilation outside of the field of view to minimise the effect of the hard field of view edge on the interpolation
+	${CARET7DIR}/wb_command -volume-dilate $DataDirectory/data.nii.gz $DilateDistance NEAREST $DataDirectory/data_dilated.nii.gz
+	# Register diffusion data to T1w space without considering gradient nonlinearities
+	${FSLDIR}/bin/flirt -in "$DataDirectory"/data_dilated -ref "$T1wRestoreImage"_${DiffRes} -applyxfm -init "$WorkingDirectory"/diff2str.mat -interp spline -out "$T1wOutputDirectory"/data
+	${FSLDIR}/bin/imrm $DataDirectory/data_dilated
+
+	# Transform CNR maps
+	${CARET7DIR}/wb_command -volume-dilate $DataDirectory/cnr_maps.nii.gz $DilateDistance NEAREST $DataDirectory/cnr_maps_dilated.nii.gz
+	${FSLDIR}/bin/flirt -in "$DataDirectory"/cnr_maps_dilated -ref "$T1wRestoreImage"_${DiffRes} -applyxfm -init "$WorkingDirectory"/diff2str.mat -interp spline -out "$T1wOutputDirectory"/cnr_maps
+	${FSLDIR}/bin/imrm $DataDirectory/cnr_maps_dilated
+
+	# Transforms field of view mask to T1-weighted space
+	${FSLDIR}/bin/flirt -in "$DataDirectory"/fov_mask -ref "$T1wRestoreImage"_${DiffRes} -applyxfm -init "$WorkingDirectory"/diff2str.mat -interp trilinear -out "$T1wOutputDirectory"/fov_mask
 fi
 
 # only include voxels fully(!) within the field of view for every volume
@@ -142,7 +140,7 @@ ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/fov_mask -thr 0.999 -bin "$T1wOutpu
 
 # Mask out data outside the brain and outside the fov
 ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/data -mas "$T1wOutputDirectory"/nodif_brain_mask_temp -mas "$T1wOutputDirectory"/fov_mask "$T1wOutputDirectory"/data
-${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/data -thr 0 "$T1wOutputDirectory"/data      #Remove negative intensity values (from eddy) from final data
+${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/data -thr 0 "$T1wOutputDirectory"/data #Remove negative intensity values (from eddy) from final data
 ${FSLDIR}/bin/imrm "$T1wOutputDirectory"/nodif_brain_mask_temp
 
 ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/cnr_maps -mas "$T1wOutputDirectory"/fov_mask "$T1wOutputDirectory"/cnr_maps
@@ -153,11 +151,11 @@ ${FSLDIR}/bin/immv "$T1wOutputDirectory"/nodif_brain_mask.nii.gz "$T1wOutputDire
 ${FSLDIR}/bin/fslmaths "$T1wOutputDirectory"/nodif_brain_mask_old.nii.gz -mas "$T1wOutputDirectory"/temp "$T1wOutputDirectory"/nodif_brain_mask
 
 # Create a simple summary text file of the percentage of spatial coverage of the dMRI data inside the FS-derived brain mask
-NvoxBrainMask=`fslstats "$T1wOutputDirectory"/nodif_brain_mask_old -V | awk '{print $1}'`
-NvoxFinalMask=`fslstats "$T1wOutputDirectory"/nodif_brain_mask -V | awk '{print $1}'`
-PctCoverage=`echo "scale=4; 100 * ${NvoxFinalMask} / ${NvoxBrainMask}" | bc -l`
-echo "PctCoverage, NvoxFinalMask, NvoxBrainMask" >| "$T1wOutputDirectory"/nodif_brain_mask.stats.txt
-echo "${PctCoverage}, ${NvoxFinalMask}, ${NvoxBrainMask}" >> "$T1wOutputDirectory"/nodif_brain_mask.stats.txt
+NvoxBrainMask=$(fslstats "$T1wOutputDirectory"/nodif_brain_mask_old -V | awk '{print $1}')
+NvoxFinalMask=$(fslstats "$T1wOutputDirectory"/nodif_brain_mask -V | awk '{print $1}')
+PctCoverage=$(echo "scale=4; 100 * ${NvoxFinalMask} / ${NvoxBrainMask}" | bc -l)
+echo "PctCoverage, NvoxFinalMask, NvoxBrainMask" >|"$T1wOutputDirectory"/nodif_brain_mask.stats.txt
+echo "${PctCoverage}, ${NvoxFinalMask}, ${NvoxBrainMask}" >>"$T1wOutputDirectory"/nodif_brain_mask.stats.txt
 
 ${FSLDIR}/bin/imrm "$T1wOutputDirectory"/temp
 ${FSLDIR}/bin/imrm "$T1wOutputDirectory"/nodif_brain_mask_old
diff --git a/DiffusionPreprocessing/scripts/basic_preproc.sh b/DiffusionPreprocessing/scripts/basic_preproc.sh
deleted file mode 100755
index be23d4367..000000000
--- a/DiffusionPreprocessing/scripts/basic_preproc.sh
+++ /dev/null
@@ -1,281 +0,0 @@
-#!/bin/bash
-source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib
-scriptName="basic_preproc.sh"
-echo -e "\n START: ${scriptName}"
-
-workingdir=$1
-echo_spacing=$2  #in msec
-PEdir=$3
-b0dist=$4
-b0maxbval=$5
-
-echo "${scriptName}: Input Parameter: workingdir: ${workingdir}"
-echo "${scriptName}: Input Parameter: echo_spacing: ${echo_spacing}"  # *Effective* Echo Spacing, in msec
-echo "${scriptName}: Input Parameter: PEdir: ${PEdir}"
-echo "${scriptName}: Input Parameter: b0dist: ${b0dist}"
-echo "${scriptName}: Input Parameter: b0maxbval: ${b0maxbval}"
-
-isodd(){
-	echo "$(( $1 % 2 ))"
-}
-
-rawdir=${workingdir}/rawdata
-topupdir=${workingdir}/topup
-eddydir=${workingdir}/eddy
-
-if [[ ${PEdir} -ne 1 && ${PEdir} -ne 2 ]] ; then
-	echo -e "\n ${scriptName}: ERROR: basic_preproc: Unrecognized PEdir: ${PEdir}"
-	exit 1
-fi
-
-# Use same convention for basePos and baseNeg names as in DiffPreprocPipeline_PreEddy.sh
-basePos="Pos"
-baseNeg="Neg"
-
-#Compute Total_readout in secs with up to 6 decimal places
-any=`ls ${rawdir}/${basePos}*.nii* | head -n 1`
-if [ ${PEdir} -eq 1 ]; then    #RL/LR phase encoding
-	dimP=`${FSLDIR}/bin/fslval ${any} dim1`
-elif [ ${PEdir} -eq 2 ]; then  #PA/AP phase encoding
-	dimP=`${FSLDIR}/bin/fslval ${any} dim2`
-fi
-dimPminus1=$(($dimP - 1))
-#Total_readout=EffectiveEchoSpacing*(ReconMatrixPE-1)
-# Factors such as in-plane acceleration, phase oversampling, phase resolution, phase field-of-view, and interpolation
-# must already be accounted for as part of the "EffectiveEchoSpacing"
-ro_time=`echo "${echo_spacing} * ${dimPminus1}" | bc -l`
-ro_time=`echo "scale=6; ${ro_time} / 1000" | bc -l`  # Convert from ms to sec
-echo "${scriptName}: Total readout time is $ro_time secs"
-
-
-################################################################################################
-## Intensity Normalisation across Series 
-################################################################################################
-echo "${scriptName}: Rescaling series to ensure consistency across baseline intensities"
-entry_cnt=0
-for entry in ${rawdir}/${basePos}_[0-9]*.nii* ${rawdir}/${baseNeg}_[0-9]*.nii*  #For each series, get the mean b0 and rescale to match the first series baseline
-do
-	basename=`imglob ${entry}`
-	echo "${scriptName}: Processing $basename"
-	
-	echo "${scriptName}: About to fslmaths ${entry} -Xmean -Ymean -Zmean ${basename}_mean"
-	${FSLDIR}/bin/fslmaths ${entry} -Xmean -Ymean -Zmean ${basename}_mean
-	if [ ! -e ${basename}_mean.nii.gz ] ; then
-		echo "${scriptName}: ERROR: Mean file: ${basename}_mean.nii.gz not created"
-		exit 1
-	fi
-	
-	echo "${scriptName}: Getting Posbvals from ${basename}.bval"
-	Posbvals=`cat ${basename}.bval`
-	echo "${scriptName}: Posbvals: ${Posbvals}"
-	
-	mcnt=0
-	for i in ${Posbvals} #extract all b0s for the series
-	do
-		echo "${scriptName}: Posbvals i: ${i}"
-		cnt=`$FSLDIR/bin/zeropad $mcnt 4`
-		echo "${scriptName}: cnt: ${cnt}"
-		if [ $i -lt ${b0maxbval} ]; then
-			echo "${scriptName}: About to fslroi ${basename}_mean ${basename}_b0_${cnt} ${mcnt} 1"
-			$FSLDIR/bin/fslroi ${basename}_mean ${basename}_b0_${cnt} ${mcnt} 1
-		fi
-		mcnt=$((${mcnt} + 1))
-	done
-	
-	echo "${scriptName}: About to fslmerge -t ${basename}_mean `echo ${basename}_b0_????.nii*`"
-	${FSLDIR}/bin/fslmerge -t ${basename}_mean `echo ${basename}_b0_????.nii*`
-	
-	echo "${scriptName}: About to fslmaths ${basename}_mean -Tmean ${basename}_mean"
-	${FSLDIR}/bin/fslmaths ${basename}_mean -Tmean ${basename}_mean #This is the mean baseline b0 intensity for the series
-	${FSLDIR}/bin/imrm ${basename}_b0_????
-	if [ ${entry_cnt} -eq 0 ]; then      #Do not rescale the first series
-		rescale=`fslmeants -i ${basename}_mean`
-	else
-		scaleS=`fslmeants -i ${basename}_mean`
-		${FSLDIR}/bin/fslmaths ${basename} -mul ${rescale} -div ${scaleS} ${basename}_new
-		${FSLDIR}/bin/imrm ${basename}   #For the rest, replace the original dataseries with the rescaled one
-		${FSLDIR}/bin/immv ${basename}_new ${basename}
-	fi
-	entry_cnt=$((${entry_cnt} + 1))
-	${FSLDIR}/bin/imrm ${basename}_mean
-done
-
-
-################################################################################################
-## b0 extraction and Creation of Index files for topup/eddy 
-################################################################################################
-echo "Extracting b0s from PE_Positive volumes and creating index and series files"
-declare -i sesdimt #declare sesdimt as integer
-tmp_indx=1
-while read line ; do  #Read SeriesCorrespVolNum.txt file
-	PCorVolNum[${tmp_indx}]=`echo $line | awk {'print $1'}`
-	tmp_indx=$((${tmp_indx}+1))
-done < ${rawdir}/${basePos}_SeriesCorrespVolNum.txt
-
-scount=1
-scount2=1
-indcount=0
-for entry in ${rawdir}/${basePos}_[0-9]*.nii*  #For each Pos volume
-do
-	#Extract b0s and create index file
-	basename=`imglob ${entry}`
-	Posbvals=`cat ${basename}.bval`
-	count=0  #Within series counter
-	count3=$((${b0dist} + 1))
-	for i in ${Posbvals}
-	do
-		if [ $count -ge ${PCorVolNum[${scount2}]} ]; then
-			tmp_ind=${indcount}
-			if [ $[tmp_ind] -eq 0 ]; then
-				tmp_ind=$((${indcount}+1))
-			fi
-			echo ${tmp_ind} >>${rawdir}/index.txt
-		else  #Consider a b=0 a volume that has a bvalue<${b0maxbval} and is at least ${b0dist} volumes away from the previous
-			if [ $i -lt ${b0maxbval} ] && [ ${count3} -gt ${b0dist} ]; then
-				cnt=`$FSLDIR/bin/zeropad $indcount 4`
-				echo "Extracting Pos Volume $count from ${entry} as a b=0. Measured b=$i" >>${rawdir}/extractedb0.txt
-				$FSLDIR/bin/fslroi ${entry} ${rawdir}/Pos_b0_${cnt} ${count} 1
-				if [ ${PEdir} -eq 1 ]; then    #RL/LR phase encoding
-					echo 1 0 0 ${ro_time} >> ${rawdir}/acqparams.txt
-				elif [ ${PEdir} -eq 2 ]; then  #AP/PA phase encoding
-					echo 0 1 0 ${ro_time} >> ${rawdir}/acqparams.txt
-				fi
-				indcount=$((${indcount} + 1))
-				count3=0
-			fi
-			echo ${indcount} >>${rawdir}/index.txt
-			count3=$((${count3} + 1))
-		fi
-		count=$((${count} + 1))
-	done
-	
-	#Create series file
-	sesdimt=`${FSLDIR}/bin/fslval ${entry} dim4` #Number of data points per Pos series
-	for (( j=0; j<${sesdimt}; j++ ))
-	do
-		echo ${scount} >> ${rawdir}/series_index.txt
-	done
-	scount=$((${scount} + 1))
-	scount2=$((${scount2} + 1))
-done
-
-echo "Extracting b0s from PE_Negative volumes and creating index and series files"
-tmp_indx=1
-while read line ; do  #Read SeriesCorrespVolNum.txt file
-	NCorVolNum[${tmp_indx}]=`echo $line | awk {'print $1'}`
-	tmp_indx=$((${tmp_indx}+1))
-done < ${rawdir}/${baseNeg}_SeriesCorrespVolNum.txt
-
-Poscount=${indcount}
-indcount=0
-scount2=1
-for entry in ${rawdir}/${baseNeg}_[0-9]*.nii* #For each Neg volume
-do
-	#Extract b0s and create index file
-	basename=`imglob ${entry}`
-	Negbvals=`cat ${basename}.bval`
-	count=0
-	count3=$((${b0dist} + 1))
-	for i in ${Negbvals}
-	do
-		if [ $count -ge ${NCorVolNum[${scount2}]} ]; then
-			tmp_ind=${indcount}
-			if [ $[tmp_ind] -eq 0 ]; then
-				tmp_ind=$((${indcount}+1))
-			fi
-			echo $((${tmp_ind} + ${Poscount})) >>${rawdir}/index.txt
-		else #Consider a b=0 a volume that has a bvalue<${b0maxbval} and is at least ${b0dist} volumes away from the previous
-			if [ $i -lt ${b0maxbval} ] && [ ${count3} -gt ${b0dist} ]; then
-				cnt=`$FSLDIR/bin/zeropad $indcount 4`
-				echo "Extracting Neg Volume $count from ${entry} as a b=0. Measured b=$i" >>${rawdir}/extractedb0.txt
-				$FSLDIR/bin/fslroi ${entry} ${rawdir}/Neg_b0_${cnt} ${count} 1
-				if [ ${PEdir} -eq 1 ]; then    #RL/LR phase encoding
-					echo -1 0 0 ${ro_time} >> ${rawdir}/acqparams.txt
-				elif [ ${PEdir} -eq 2 ]; then  #AP/PA phase encoding
-					echo 0 -1 0 ${ro_time} >> ${rawdir}/acqparams.txt
-				fi
-				indcount=$((${indcount} + 1))
-				count3=0
-			fi
-			echo $((${indcount} + ${Poscount})) >>${rawdir}/index.txt
-			count3=$((${count3} + 1))
-		fi
-		count=$((${count} + 1))
-	done
-	
-	#Create series file
-	sesdimt=`${FSLDIR}/bin/fslval ${entry} dim4`
-	for (( j=0; j<${sesdimt}; j++ ))
-	do
-		echo ${scount} >> ${rawdir}/series_index.txt #Create series file
-	done
-	scount=$((${scount} + 1))
-	scount2=$((${scount2} + 1))
-done
-
-################################################################################################
-## Merging Files and correct number of slices 
-################################################################################################
-echo "Merging Pos and Neg images"
-${FSLDIR}/bin/fslmerge -t ${rawdir}/Pos_b0 `${FSLDIR}/bin/imglob ${rawdir}/Pos_b0_????.*`
-${FSLDIR}/bin/fslmerge -t ${rawdir}/Neg_b0 `${FSLDIR}/bin/imglob ${rawdir}/Neg_b0_????.*`
-${FSLDIR}/bin/imrm ${rawdir}/Pos_b0_????
-${FSLDIR}/bin/imrm ${rawdir}/Neg_b0_????
-${FSLDIR}/bin/fslmerge -t ${rawdir}/Pos `echo ${rawdir}/${basePos}_[0-9]*.nii*`
-${FSLDIR}/bin/fslmerge -t ${rawdir}/Neg `echo ${rawdir}/${baseNeg}_[0-9]*.nii*`
-
-paste `echo ${rawdir}/${basePos}*.bval` >${rawdir}/Pos.bval
-paste `echo ${rawdir}/${basePos}*.bvec` >${rawdir}/Pos.bvec
-paste `echo ${rawdir}/${baseNeg}*.bval` >${rawdir}/Neg.bval
-paste `echo ${rawdir}/${baseNeg}*.bvec` >${rawdir}/Neg.bvec
-
-
-dimz=`${FSLDIR}/bin/fslval ${rawdir}/Pos dim3`
-if [ `isodd $dimz` -eq 1 ];then
-	echo "Remove one slice from data to get even number of slices"
-	${FSLDIR}/bin/fslroi ${rawdir}/Pos ${rawdir}/Posn 0 -1 0 -1 1 -1
-	${FSLDIR}/bin/fslroi ${rawdir}/Neg ${rawdir}/Negn 0 -1 0 -1 1 -1
-	${FSLDIR}/bin/fslroi ${rawdir}/Pos_b0 ${rawdir}/Pos_b0n 0 -1 0 -1 1 -1
-	${FSLDIR}/bin/fslroi ${rawdir}/Neg_b0 ${rawdir}/Neg_b0n 0 -1 0 -1 1 -1
-	${FSLDIR}/bin/imrm ${rawdir}/Pos
-	${FSLDIR}/bin/imrm ${rawdir}/Neg
-	${FSLDIR}/bin/imrm ${rawdir}/Pos_b0
-	${FSLDIR}/bin/imrm ${rawdir}/Neg_b0
-	${FSLDIR}/bin/immv ${rawdir}/Posn ${rawdir}/Pos
-	${FSLDIR}/bin/immv ${rawdir}/Negn ${rawdir}/Neg
-	${FSLDIR}/bin/immv ${rawdir}/Pos_b0n ${rawdir}/Pos_b0
-	${FSLDIR}/bin/immv ${rawdir}/Neg_b0n ${rawdir}/Neg_b0
-fi
-
-echo "Perform final merge"
-${FSLDIR}/bin/fslmerge -t ${rawdir}/Pos_Neg_b0 ${rawdir}/Pos_b0 ${rawdir}/Neg_b0 
-${FSLDIR}/bin/fslmerge -t ${rawdir}/Pos_Neg ${rawdir}/Pos ${rawdir}/Neg
-paste ${rawdir}/Pos.bval ${rawdir}/Neg.bval >${rawdir}/Pos_Neg.bvals
-paste ${rawdir}/Pos.bvec ${rawdir}/Neg.bvec >${rawdir}/Pos_Neg.bvecs
-
-${FSLDIR}/bin/imrm ${rawdir}/Pos
-${FSLDIR}/bin/imrm ${rawdir}/Neg
-
-
-################################################################################################
-## Move files to appropriate directories 
-################################################################################################
-echo "Move files to appropriate directories"
-mv ${rawdir}/extractedb0.txt ${topupdir}
-mv ${rawdir}/acqparams.txt ${topupdir}
-${FSLDIR}/bin/immv ${rawdir}/Pos_Neg_b0 ${topupdir}
-${FSLDIR}/bin/immv ${rawdir}/Pos_b0 ${topupdir}
-${FSLDIR}/bin/immv ${rawdir}/Neg_b0 ${topupdir}
-
-cp ${topupdir}/acqparams.txt ${eddydir}
-mv ${rawdir}/index.txt ${eddydir}
-mv ${rawdir}/series_index.txt ${eddydir}
-${FSLDIR}/bin/immv ${rawdir}/Pos_Neg ${eddydir}
-mv ${rawdir}/Pos_Neg.bvals ${eddydir}
-mv ${rawdir}/Pos_Neg.bvecs ${eddydir}
-mv ${rawdir}/Pos.bv?? ${eddydir}
-mv ${rawdir}/Neg.bv?? ${eddydir}
-
-echo -e "\n END: basic_preproc"
-
-
diff --git a/DiffusionPreprocessing/scripts/basic_preproc_best_b0.sh b/DiffusionPreprocessing/scripts/basic_preproc_best_b0.sh
new file mode 100755
index 000000000..a2803786c
--- /dev/null
+++ b/DiffusionPreprocessing/scripts/basic_preproc_best_b0.sh
@@ -0,0 +1,292 @@
+#!/bin/bash
+source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib
+scriptName="basic_preproc_best_b0.sh"
+echo -e "\n START: ${scriptName}"
+
+workingdir=$1
+ro_time=$2 #in sec
+PEdir=$3
+b0maxbval=$4
+
+echo "${scriptName}: Input Parameter: workingdir: ${workingdir}"
+echo "${scriptName}: Input Parameter: ro_time: ${ro_time}" # Readout time in sec
+echo "${scriptName}: Input Parameter: PEdir: ${PEdir}"
+echo "${scriptName}: Input Parameter: b0maxbval: ${b0maxbval}"
+
+isodd() {
+	echo "$(($1 % 2))"
+}
+
+rawdir=${workingdir}/rawdata
+topupdir=${workingdir}/topup
+eddydir=${workingdir}/eddy
+
+if [[ ${PEdir} -ne 1 && ${PEdir} -ne 2 ]]; then
+	echo -e "\n ERROR: ${scriptName}: Unrecognized PEdir: ${PEdir}"
+	exit 1
+fi
+
+# Use same convention for basePos and baseNeg names as in DiffPreprocPipeline_PreEddy.sh
+basePos="Pos"
+baseNeg="Neg"
+
+################################################################################################
+## Identifying the best b0's to use in topup
+################################################################################################
+# This code in this section was adapted from a script written for the developing HCP (https://git.fmrib.ox.ac.uk/matteob/dHCP_neo_dMRI_pipeline_release/blob/master/utils/pickBestB0s.sh)
+# The original script was released under the Apache license 2.0 (https://git.fmrib.ox.ac.uk/matteob/dHCP_neo_dMRI_pipeline_release/blob/master/LICENSE)
+
+# Merge b0's for both phase encoding directions
+for pe_sign in ${basePos} ${baseNeg}; do
+	merge_command=("${FSLDIR}/bin/fslmerge" -t "${rawdir}/all_${pe_sign}_b0s")
+	for entry in ${rawdir}/${pe_sign}_[0-9]*.nii*; do
+		basename=$(imglob ${entry})
+		# TODO: replace with FSL built-in version of select_dwi_vols once -db flag is supported (should be in 6.0.4)
+		${HCPPIPEDIR_Global}/select_dwi_vols ${basename} ${basename}.bval ${basename}_b0s 0 -db ${b0maxbval}
+		merge_command+=("${basename}_b0s")
+	done
+	echo about to "${merge_command[@]}"
+	"${merge_command[@]}"
+	for entry in ${rawdir}/${pe_sign}_[0-9]*_b0s.nii*; do
+		${FSLDIR}/bin/imrm ${entry}
+	done
+done
+
+# Here we identify the b0's that are least affected by motion artefacts to pass them on to topup.
+# These b0's are identified by being most similar to a reference b0, which is determined in one of two ways:
+# 1. If there are enough b0's with a given phase encoding direction (>= 5) we adopt the average b0 as a reference.
+# 2. If there are fewer b0's, then the average b0 might be contaminated by any motion artefacts in one or two b0's.
+#    In that case we use topup to combine the b0's with the b0's with opposite phase encoding to get a reference b0.
+# To further reduce the chance that our reference b0 is contaminated by motion we compute the average multiple times,
+# each time only using those b0's that were most similar to the previous average b0 (and hence least likely to be
+# affected by motion)
+for pe_sign in ${basePos} ${baseNeg}; do
+	echo "Identifying best b0's for ${pe_sign} phase encoding"
+	if [ ${pe_sign} = ${basePos} ]; then
+		pe_other=${baseNeg}
+	else
+		pe_other=${basePos}
+	fi
+	select_b0_dir=${rawdir}/select_b0_${pe_sign}
+	mkdir -p ${select_b0_dir}
+
+	N_b0s=$(${FSLDIR}/bin/fslval ${rawdir}/all_${pe_sign}_b0s dim4)
+
+	# if there are less than 5 B0's for a specific phase encoding use topup to find the best ones
+	# otherwise simply register the B0's of the same phase encoding to each other
+	if [[ ${N_b0s} -lt 5 ]]; then
+		N_other=$(${FSLDIR}/bin/fslval ${rawdir}/all_${pe_other}_b0s dim4)
+		if [[ ${N_other} -gt 4 ]]; then N_other=4; fi
+		echo "Score all ${N_b0s} ${pe_sign} B0's based on alignment with mean B0 after topup with ${N_other} ${pe_other} B0's"
+
+		# Select sub-set of other B0's to run topup on; we use the first ${N_other}
+		${FSLDIR}/bin/fslroi ${rawdir}/all_${pe_other}_b0s ${select_b0_dir}/opposite_b0s 0 ${N_other}
+
+		# Merge all b0's to do a rough initial aligment using topup
+		${FSLDIR}/bin/fslmerge -t ${select_b0_dir}/all_b0s ${rawdir}/all_${pe_sign}_b0s ${select_b0_dir}/opposite_b0s
+
+		# Create the acqparams file for the initial alignment
+		for idx in $(seq 1 ${N_b0s}); do
+			if [ ${PEdir} -eq 1 ]; then #RL/LR phase encoding
+				echo 1 0 0 ${ro_time} >>${select_b0_dir}/acqparams.txt
+			elif [ ${PEdir} -eq 2 ]; then #AP/PA phase encoding
+				echo 0 1 0 ${ro_time} >>${select_b0_dir}/acqparams.txt
+			fi
+		done
+
+		for idx in $(seq 1 ${N_other}); do
+			if [ ${PEdir} -eq 1 ]; then #RL/LR phase encoding
+				echo -1 0 0 ${ro_time} >>${select_b0_dir}/acqparams.txt
+			elif [ ${PEdir} -eq 2 ]; then #AP/PA phase encoding
+				echo 0 -1 0 ${ro_time} >>${select_b0_dir}/acqparams.txt
+			fi
+		done
+
+		dimz=$(${FSLDIR}/bin/fslval ${select_b0_dir}/all_b0s dim3)
+		if [ $(isodd $dimz) -eq 1 ]; then
+			${FSLDIR}/bin/fslroi ${select_b0_dir}/all_b0s ${select_b0_dir}/all_b0s_even 0 -1 0 -1 1 -1
+		else
+			${FSLDIR}/bin/imcp ${select_b0_dir}/all_b0s ${select_b0_dir}/all_b0s_even
+		fi
+		# run topup to roughly align the b0's
+		configdir=${HCPPIPEDIR_Config}
+		topup_config_file=${configdir}/best_b0.cnf
+		${FSLDIR}/bin/topup --imain=${select_b0_dir}/all_b0s_even \
+			--datain=${select_b0_dir}/acqparams.txt \
+			--config=${topup_config_file} \
+			--fout=${select_b0_dir}/fieldmap \
+			--iout=${select_b0_dir}/topup_b0s \
+			--out=${select_b0_dir}/topup_results \
+			-v
+
+		# compute squared residual from the mean b0
+		# once "bad" b0's are identified, the mean b0 will be recomputed without them
+		# and the scores recomputed; this process is iterated 3 times
+		${FSLDIR}/bin/fslmaths ${select_b0_dir}/topup_b0s -Tmean ${select_b0_dir}/topup_b0s_avg
+		for ((i = 1; i <= 3; i++)); do
+			${FSLDIR}/bin/fslmaths ${select_b0_dir}/topup_b0s -sub ${select_b0_dir}/topup_b0s_avg -sqr ${select_b0_dir}/topup_b0s_res
+
+			# Get brain mask from averaged results
+			${FSLDIR}/bin/bet ${select_b0_dir}/topup_b0s_avg ${select_b0_dir}/nodif_brain -m -R -f 0.3
+
+			# compute average squared residual over brain mask
+			scores=($(${FSLDIR}/bin/fslstats -t ${select_b0_dir}/topup_b0s_res -k ${select_b0_dir}/nodif_brain_mask -M))
+			scores_str="${scores[@]}"
+
+			echo "Iteration ${i} in finding best b0 for ${pe_sign}"
+			echo "Current scores: ${scores_str}"
+
+			# Recomputes the average using only the b0's with scores below (median(score) + 2 * mad(score)),
+			# where mad is the median absolute deviation.
+			idx=$(fslpython -c "from numpy import median, where; sc = [float(s) for s in '${scores_str}'.split()]; print(','.join(str(idx) for idx in where(sc < median(abs(sc - median(sc)) * 2 + median(sc)))[0]))")
+			echo "Recomputing average b0 using indices (counting from zero): ${idx}"
+			${FSLDIR}/bin/fslselectvols -i ${select_b0_dir}/topup_b0s -o ${select_b0_dir}/topup_b0s_avg --vols=${idx} -m
+		done
+		# recompute the squared residuals and brainmask using the final average
+		${FSLDIR}/bin/fslmaths ${select_b0_dir}/topup_b0s -sub ${select_b0_dir}/topup_b0s_avg -sqr ${select_b0_dir}/topup_b0s_res
+		${FSLDIR}/bin/bet ${select_b0_dir}/topup_b0s_avg ${select_b0_dir}/nodif_brain -m -R -f 0.3
+
+		# select only the polarity of interest and compute the final scores (average squared residual over the brainmask)
+		${FSLDIR}/bin/fslroi ${select_b0_dir}/topup_b0s_res ${select_b0_dir}/topup_b0s_res_${pe_sign} 0 ${N_b0s}
+		scores=($(${FSLDIR}/bin/fslstats -t ${select_b0_dir}/topup_b0s_res_${pe_sign} -k ${select_b0_dir}/nodif_brain_mask -M))
+	else # Number of b0's with this phase encoding (${N_b0s}) >= 5
+		echo "Score all ${pe_sign} B0's based on similarity with the mean ${pe_sign} B0"
+		echo ${FSLDIR}/bin/mcflirt -in ${rawdir}/all_${pe_sign}_b0s -out ${select_b0_dir}/all_b0s_mcf
+		${FSLDIR}/bin/mcflirt -in ${rawdir}/all_${pe_sign}_b0s -out ${select_b0_dir}/all_b0s_mcf
+		${FSLDIR}/bin/fslmaths ${select_b0_dir}/all_b0s_mcf -Tmean ${select_b0_dir}/all_b0s_mcf_avg
+
+		# compute squared residual from the mean b0
+		# once "bad" b0's are identified, the mean b0 will be recomputed without them
+		# and the scores recomputed; this process is iterated 3 times
+		for ((i = 1; i <= 3; i++)); do
+			${FSLDIR}/bin/fslmaths ${select_b0_dir}/all_b0s_mcf -sub ${select_b0_dir}/all_b0s_mcf_avg -sqr ${select_b0_dir}/all_b0s_mcf_res
+
+			# Get brain mask from averaged results
+			${FSLDIR}/bin/bet ${select_b0_dir}/all_b0s_mcf_avg ${select_b0_dir}/nodif_brain -m -R -f 0.3
+			scores=($(${FSLDIR}/bin/fslstats -t ${select_b0_dir}/all_b0s_mcf_res -k ${select_b0_dir}/nodif_brain_mask -M))
+			scores_str="${scores[@]}"
+
+			echo "Iteration ${i} in finding best b0 for ${pe_sign}"
+			echo "Current scores: ${scores_str}"
+
+			# Recomputes the average using only the b0's with scores below (median(score) + 2 * mad(score)),
+			# where mad is the median absolute deviation.
+			idx=$(fslpython -c "from numpy import median, where; sc = [float(s) for s in '${scores_str}'.split()]; print(','.join(str(idx) for idx in where(sc < median(abs(sc - median(sc)) * 2 + median(sc)))[0]))")
+			echo "Recomputing average b0 using indices (counting from zero): ${idx}"
+			${FSLDIR}/bin/fslselectvols -i ${select_b0_dir}/all_b0s_mcf -o ${select_b0_dir}/all_b0s_mcf_avg --vols=${idx} -m
+		done
+	fi
+	echo "Final b0 scores for ${pe_sign}: " "${scores[@]}"
+
+	printf "%s\n" "${scores[@]}" >${select_b0_dir}/scores.txt
+done
+
+################################################################################################
+## b0 extraction and Creation of acquisition paramater file for topup/eddy
+################################################################################################
+echo "Find the best B0 in the positive and negative volumes"
+
+rm -f ${rawdir}/index_best_b0s.txt
+
+for pe_sign in ${basePos} ${baseNeg}; do
+	# find index of minimum score
+	scores=()
+	while read line; do scores+=("$line"); done <${rawdir}/select_b0_${pe_sign}/scores.txt
+	min_idx=0
+	for idx in $(seq 0 $(($(${FSLDIR}/bin/fslval ${rawdir}/all_${pe_sign}_b0s dim4) - 1))); do
+		if [ $(echo "${scores[${idx}]} < ${scores[${min_idx}]}" | bc -l) -eq 1 ]; then min_idx=$idx; fi
+	done
+	echo "Selecting ${pe_sign} B0 with index ${min_idx} (counting from zero)"
+	echo "${pe_sign} ${min_idx}" >>${rawdir}/index_best_b0s.txt
+	${FSLDIR}/bin/fslroi ${rawdir}/all_${pe_sign}_b0s ${rawdir}/best_${pe_sign}_b0 ${min_idx} 1
+
+	if [ ${pe_sign} = ${basePos} ]; then
+		# store this to extract the b-value later
+		best_pos_b0_index=${min_idx}
+	fi
+done
+
+# producing acqparams.txt
+if [ ${PEdir} -eq 1 ]; then #RL/LR phase encoding
+	echo 1 0 0 ${ro_time} >${rawdir}/acqparams.txt
+	echo -1 0 0 ${ro_time} >>${rawdir}/acqparams.txt
+elif [ ${PEdir} -eq 2 ]; then #AP/PA phase encoding
+	echo 0 1 0 ${ro_time} >${rawdir}/acqparams.txt
+	echo 0 -1 0 ${ro_time} >>${rawdir}/acqparams.txt
+fi
+
+################################################################################################
+## Merging Files and correct number of slices
+################################################################################################
+echo "Merging Pos and Neg images"
+${FSLDIR}/bin/fslmerge -t ${rawdir}/Pos $(echo ${rawdir}/${basePos}_[0-9]*.nii*)
+${FSLDIR}/bin/fslmerge -t ${rawdir}/Neg $(echo ${rawdir}/${baseNeg}_[0-9]*.nii*)
+
+paste -d' ' $(echo ${rawdir}/${basePos}_[0-9]*.bval) >${rawdir}/Pos.bval
+paste -d' ' $(echo ${rawdir}/${basePos}_[0-9]*.bvec) >${rawdir}/Pos.bvec
+paste -d' ' $(echo ${rawdir}/${baseNeg}_[0-9]*.bval) >${rawdir}/Neg.bval
+paste -d' ' $(echo ${rawdir}/${baseNeg}_[0-9]*.bvec) >${rawdir}/Neg.bvec
+
+# start index file with a 1 to indicate the reference B0 image
+echo 1 >${rawdir}/index.txt
+for idx in $(seq 1 $(${FSLDIR}/bin/fslval ${rawdir}/Pos dim4)); do
+	echo 1 >>${rawdir}/index.txt
+done
+for idx in $(seq 1 $(${FSLDIR}/bin/fslval ${rawdir}/Neg dim4)); do
+	echo 2 >>${rawdir}/index.txt
+done
+
+dimz=$(${FSLDIR}/bin/fslval ${rawdir}/Pos dim3)
+if [ $(isodd $dimz) -eq 1 ]; then
+	echo "Remove one slice from data to get even number of slices"
+	${FSLDIR}/bin/fslroi ${rawdir}/Pos ${rawdir}/Posn 0 -1 0 -1 1 -1
+	${FSLDIR}/bin/fslroi ${rawdir}/Neg ${rawdir}/Negn 0 -1 0 -1 1 -1
+	${FSLDIR}/bin/fslroi ${rawdir}/best_Pos_b0 ${rawdir}/best_Pos_b0n 0 -1 0 -1 1 -1
+	${FSLDIR}/bin/fslroi ${rawdir}/best_Neg_b0 ${rawdir}/best_Neg_b0n 0 -1 0 -1 1 -1
+	${FSLDIR}/bin/imrm ${rawdir}/Pos
+	${FSLDIR}/bin/imrm ${rawdir}/Neg
+	${FSLDIR}/bin/imrm ${rawdir}/best_Pos_b0
+	${FSLDIR}/bin/imrm ${rawdir}/best_Neg_b0
+	${FSLDIR}/bin/immv ${rawdir}/Posn ${rawdir}/Pos
+	${FSLDIR}/bin/immv ${rawdir}/Negn ${rawdir}/Neg
+	${FSLDIR}/bin/immv ${rawdir}/best_Pos_b0n ${rawdir}/best_Pos_b0
+	${FSLDIR}/bin/immv ${rawdir}/best_Neg_b0n ${rawdir}/best_Neg_b0
+fi
+
+echo "Perform final merge"
+${FSLDIR}/bin/fslmerge -t ${rawdir}/Pos_Neg_b0 ${rawdir}/best_Pos_b0 ${rawdir}/best_Neg_b0
+# include Pos_b0 as the first volume of Pos_Neg, so that eddy will use it as reference
+# it will be removed after eddy is completed
+${FSLDIR}/bin/fslmerge -t ${rawdir}/Pos_Neg ${rawdir}/best_Pos_b0 ${rawdir}/Pos ${rawdir}/Neg
+
+# extract b-value and bvec of best b0 in Pos set
+best_pos_bval=$(cat ${rawdir}/Pos.bval | awk "{print \$$((${best_pos_b0_index} + 1))"})
+cat ${rawdir}/Pos.bvec | awk "{print \$$((${best_pos_b0_index} + 1))}" >${rawdir}/zero.bvecs
+
+# merge all b-values and bvecs
+echo $best_pos_bval $(paste -d' ' ${rawdir}/Pos.bval ${rawdir}/Neg.bval) >${rawdir}/Pos_Neg.bvals
+paste -d' ' ${rawdir}/zero.bvecs ${rawdir}/Pos.bvec ${rawdir}/Neg.bvec >${rawdir}/Pos_Neg.bvecs
+
+rm ${rawdir}/zero.bvecs
+${FSLDIR}/bin/imrm ${rawdir}/Pos
+${FSLDIR}/bin/imrm ${rawdir}/Neg
+
+################################################################################################
+## Move files to appropriate directories
+################################################################################################
+echo "Move files to appropriate directories"
+mv ${rawdir}/index_best_b0s.txt ${topupdir}
+mv ${rawdir}/acqparams.txt ${topupdir}
+${FSLDIR}/bin/immv ${rawdir}/Pos_Neg_b0 ${topupdir}
+${FSLDIR}/bin/immv ${rawdir}/best_Pos_b0 ${topupdir}/Pos_b0
+${FSLDIR}/bin/immv ${rawdir}/best_Neg_b0 ${topupdir}/Neg_b0
+
+cp ${topupdir}/acqparams.txt ${eddydir}
+mv ${rawdir}/index.txt ${eddydir}
+${FSLDIR}/bin/immv ${rawdir}/Pos_Neg ${eddydir}
+mv ${rawdir}/Pos_Neg.bvals ${eddydir}
+mv ${rawdir}/Pos_Neg.bvecs ${eddydir}
+mv ${rawdir}/Pos.bv?? ${eddydir}
+mv ${rawdir}/Neg.bv?? ${eddydir}
+
+echo -e "\n END: ${scriptName}"
diff --git a/DiffusionPreprocessing/scripts/basic_preproc_norm_intensity.sh b/DiffusionPreprocessing/scripts/basic_preproc_norm_intensity.sh
new file mode 100755
index 000000000..9887fd416
--- /dev/null
+++ b/DiffusionPreprocessing/scripts/basic_preproc_norm_intensity.sh
@@ -0,0 +1,70 @@
+#!/bin/bash
+source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib
+scriptName="basic_preproc_norm_intensity.sh"
+echo -e "\n START: ${scriptName}"
+
+workingdir=$1
+b0maxbval=$2
+
+echo "${scriptName}: Input Parameter: workingdir: ${workingdir}"
+echo "${scriptName}: Input Parameter: b0maxbval: ${b0maxbval}"
+
+# Use same convention for basePos and baseNeg names as in DiffPreprocPipeline_PreEddy.sh
+basePos="Pos"
+baseNeg="Neg"
+
+rawdir=${workingdir}/rawdata
+topupdir=${workingdir}/topup
+eddydir=${workingdir}/eddy
+
+################################################################################################
+## Intensity Normalisation across Series
+################################################################################################
+echo "${scriptName}: Rescaling series to ensure consistency across baseline intensities"
+entry_cnt=0
+for entry in ${rawdir}/${basePos}_[0-9]*.nii* ${rawdir}/${baseNeg}_[0-9]*.nii*; do #For each series, get the mean b0 and rescale to match the first series baseline
+	basename=$(imglob ${entry})
+	echo "${scriptName}: Processing $basename"
+
+	echo "${scriptName}: About to fslmaths ${entry} -Xmean -Ymean -Zmean ${basename}_mean"
+	${FSLDIR}/bin/fslmaths ${entry} -Xmean -Ymean -Zmean ${basename}_mean
+	if [ ! -e ${basename}_mean.nii.gz ]; then
+		echo "${scriptName}: ERROR: Mean file: ${basename}_mean.nii.gz not created"
+		exit 1
+	fi
+
+	echo "${scriptName}: Getting Posbvals from ${basename}.bval"
+	Posbvals=$(cat ${basename}.bval)
+	echo "${scriptName}: Posbvals: ${Posbvals}"
+
+	mcnt=0
+	for i in ${Posbvals}; do #extract all b0s for the series
+		echo "${scriptName}: Posbvals i: ${i}"
+		cnt=$($FSLDIR/bin/zeropad $mcnt 4)
+		echo "${scriptName}: cnt: ${cnt}"
+		if [ $i -lt ${b0maxbval} ]; then
+			echo "${scriptName}: About to fslroi ${basename}_mean ${basename}_b0_${cnt} ${mcnt} 1"
+			$FSLDIR/bin/fslroi ${basename}_mean ${basename}_b0_${cnt} ${mcnt} 1
+		fi
+		mcnt=$((${mcnt} + 1))
+	done
+
+	echo "${scriptName}: About to fslmerge -t ${basename}_mean $(echo ${basename}_b0_????.nii*)"
+	${FSLDIR}/bin/fslmerge -t ${basename}_mean $(echo ${basename}_b0_????.nii*)
+
+	echo "${scriptName}: About to fslmaths ${basename}_mean -Tmean ${basename}_mean"
+	${FSLDIR}/bin/fslmaths ${basename}_mean -Tmean ${basename}_mean #This is the mean baseline b0 intensity for the series
+	${FSLDIR}/bin/imrm ${basename}_b0_????
+	if [ ${entry_cnt} -eq 0 ]; then #Do not rescale the first series
+		rescale=$(fslmeants -i ${basename}_mean)
+	else
+		scaleS=$(fslmeants -i ${basename}_mean)
+		${FSLDIR}/bin/fslmaths ${basename} -mul ${rescale} -div ${scaleS} ${basename}_new
+		${FSLDIR}/bin/imrm ${basename} #For the rest, replace the original dataseries with the rescaled one
+		${FSLDIR}/bin/immv ${basename}_new ${basename}
+	fi
+	entry_cnt=$((${entry_cnt} + 1))
+	${FSLDIR}/bin/imrm ${basename}_mean
+done
+
+echo -e "\n END: ${scriptName}"
diff --git a/DiffusionPreprocessing/scripts/basic_preproc_sequence.sh b/DiffusionPreprocessing/scripts/basic_preproc_sequence.sh
new file mode 100755
index 000000000..ec8da810f
--- /dev/null
+++ b/DiffusionPreprocessing/scripts/basic_preproc_sequence.sh
@@ -0,0 +1,201 @@
+#!/bin/bash
+source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib
+scriptName="basic_preproc_sequence.sh"
+echo -e "\n START: ${scriptName}"
+
+workingdir=$1
+ro_time=$2 #in sec
+PEdir=$3
+b0dist=$4
+b0maxbval=$5
+
+echo "${scriptName}: Input Parameter: workingdir: ${workingdir}"
+echo "${scriptName}: Input Parameter: ro_time: ${ro_time}" # Readout time in sec
+echo "${scriptName}: Input Parameter: PEdir: ${PEdir}"
+echo "${scriptName}: Input Parameter: b0dist: ${b0dist}"
+echo "${scriptName}: Input Parameter: b0maxbval: ${b0maxbval}"
+isodd() {
+	echo "$(($1 % 2))"
+}
+
+rawdir=${workingdir}/rawdata
+topupdir=${workingdir}/topup
+eddydir=${workingdir}/eddy
+
+if [[ ${PEdir} -ne 1 && ${PEdir} -ne 2 ]]; then
+	echo -e "\n ${scriptName}: ERROR: basic_preproc: Unrecognized PEdir: ${PEdir}"
+	exit 1
+fi
+
+# Use same convention for basePos and baseNeg names as in DiffPreprocPipeline_PreEddy.sh
+basePos="Pos"
+baseNeg="Neg"
+
+################################################################################################
+## b0 extraction and Creation of Index files for topup/eddy
+################################################################################################
+echo "Extracting b0s from PE_Positive volumes and creating index and series files"
+declare -i sesdimt #declare sesdimt as integer
+tmp_indx=1
+while read line; do #Read SeriesCorrespVolNum.txt file
+	PCorVolNum[${tmp_indx}]=$(echo $line | awk {'print $1'})
+	tmp_indx=$((${tmp_indx} + 1))
+done <${rawdir}/${basePos}_SeriesCorrespVolNum.txt
+
+scount=1
+scount2=1
+indcount=0
+for entry in ${rawdir}/${basePos}_[0-9]*.nii*; do #For each Pos volume
+	#Extract b0s and create index file
+	basename=$(imglob ${entry})
+	Posbvals=$(cat ${basename}.bval)
+	count=0 #Within series counter
+	count3=$((${b0dist} + 1))
+	for i in ${Posbvals}; do
+		if [ $count -ge ${PCorVolNum[${scount2}]} ]; then
+			tmp_ind=${indcount}
+			if [ $((tmp_ind)) -eq 0 ]; then
+				tmp_ind=$((${indcount} + 1))
+			fi
+			echo ${tmp_ind} >>${rawdir}/index.txt
+		else #Consider a b=0 a volume that has a bvalue<${b0maxbval} and is at least ${b0dist} volumes away from the previous
+			if [ $i -lt ${b0maxbval} ] && [ ${count3} -gt ${b0dist} ]; then
+				cnt=$($FSLDIR/bin/zeropad $indcount 4)
+				echo "Extracting Pos Volume $count from ${entry} as a b=0. Measured b=$i" >>${rawdir}/extractedb0.txt
+				$FSLDIR/bin/fslroi ${entry} ${rawdir}/Pos_b0_${cnt} ${count} 1
+				if [ ${PEdir} -eq 1 ]; then #RL/LR phase encoding
+					echo 1 0 0 ${ro_time} >>${rawdir}/acqparams.txt
+				elif [ ${PEdir} -eq 2 ]; then #AP/PA phase encoding
+					echo 0 1 0 ${ro_time} >>${rawdir}/acqparams.txt
+				fi
+				indcount=$((${indcount} + 1))
+				count3=0
+			fi
+			echo ${indcount} >>${rawdir}/index.txt
+			count3=$((${count3} + 1))
+		fi
+		count=$((${count} + 1))
+	done
+
+	#Create series file
+	sesdimt=$(${FSLDIR}/bin/fslval ${entry} dim4) #Number of data points per Pos series
+	for ((j = 0; j < ${sesdimt}; j++)); do
+		echo ${scount} >>${rawdir}/series_index.txt
+	done
+	scount=$((${scount} + 1))
+	scount2=$((${scount2} + 1))
+done
+
+echo "Extracting b0s from PE_Negative volumes and creating index and series files"
+tmp_indx=1
+while read line; do #Read SeriesCorrespVolNum.txt file
+	NCorVolNum[${tmp_indx}]=$(echo $line | awk {'print $1'})
+	tmp_indx=$((${tmp_indx} + 1))
+done <${rawdir}/${baseNeg}_SeriesCorrespVolNum.txt
+
+Poscount=${indcount}
+indcount=0
+scount2=1
+for entry in ${rawdir}/${baseNeg}_[0-9]*.nii*; do #For each Neg volume
+	#Extract b0s and create index file
+	basename=$(imglob ${entry})
+	Negbvals=$(cat ${basename}.bval)
+	count=0
+	count3=$((${b0dist} + 1))
+	for i in ${Negbvals}; do
+		if [ $count -ge ${NCorVolNum[${scount2}]} ]; then
+			tmp_ind=${indcount}
+			if [ $((tmp_ind)) -eq 0 ]; then
+				tmp_ind=$((${indcount} + 1))
+			fi
+			echo $((${tmp_ind} + ${Poscount})) >>${rawdir}/index.txt
+		else #Consider a b=0 a volume that has a bvalue<${b0maxbval} and is at least ${b0dist} volumes away from the previous
+			if [ $i -lt ${b0maxbval} ] && [ ${count3} -gt ${b0dist} ]; then
+				cnt=$($FSLDIR/bin/zeropad $indcount 4)
+				echo "Extracting Neg Volume $count from ${entry} as a b=0. Measured b=$i" >>${rawdir}/extractedb0.txt
+				$FSLDIR/bin/fslroi ${entry} ${rawdir}/Neg_b0_${cnt} ${count} 1
+				if [ ${PEdir} -eq 1 ]; then #RL/LR phase encoding
+					echo -1 0 0 ${ro_time} >>${rawdir}/acqparams.txt
+				elif [ ${PEdir} -eq 2 ]; then #AP/PA phase encoding
+					echo 0 -1 0 ${ro_time} >>${rawdir}/acqparams.txt
+				fi
+				indcount=$((${indcount} + 1))
+				count3=0
+			fi
+			echo $((${indcount} + ${Poscount})) >>${rawdir}/index.txt
+			count3=$((${count3} + 1))
+		fi
+		count=$((${count} + 1))
+	done
+
+	#Create series file
+	sesdimt=$(${FSLDIR}/bin/fslval ${entry} dim4)
+	for ((j = 0; j < ${sesdimt}; j++)); do
+		echo ${scount} >>${rawdir}/series_index.txt #Create series file
+	done
+	scount=$((${scount} + 1))
+	scount2=$((${scount2} + 1))
+done
+
+################################################################################################
+## Merging Files and correct number of slices
+################################################################################################
+echo "Merging Pos and Neg images"
+${FSLDIR}/bin/fslmerge -t ${rawdir}/Pos_b0 $(${FSLDIR}/bin/imglob ${rawdir}/Pos_b0_????.*)
+${FSLDIR}/bin/fslmerge -t ${rawdir}/Neg_b0 $(${FSLDIR}/bin/imglob ${rawdir}/Neg_b0_????.*)
+${FSLDIR}/bin/imrm ${rawdir}/Pos_b0_????
+${FSLDIR}/bin/imrm ${rawdir}/Neg_b0_????
+${FSLDIR}/bin/fslmerge -t ${rawdir}/Pos $(echo ${rawdir}/${basePos}_[0-9]*.nii*)
+${FSLDIR}/bin/fslmerge -t ${rawdir}/Neg $(echo ${rawdir}/${baseNeg}_[0-9]*.nii*)
+
+paste $(echo ${rawdir}/${basePos}_[0-9]*.bval) >${rawdir}/Pos.bval
+paste $(echo ${rawdir}/${basePos}_[0-9]*.bvec) >${rawdir}/Pos.bvec
+paste $(echo ${rawdir}/${baseNeg}_[0-9]*.bval) >${rawdir}/Neg.bval
+paste $(echo ${rawdir}/${baseNeg}_[0-9]*.bvec) >${rawdir}/Neg.bvec
+
+dimz=$(${FSLDIR}/bin/fslval ${rawdir}/Pos dim3)
+if [ $(isodd $dimz) -eq 1 ]; then
+	echo "Remove one slice from data to get even number of slices"
+	${FSLDIR}/bin/fslroi ${rawdir}/Pos ${rawdir}/Posn 0 -1 0 -1 1 -1
+	${FSLDIR}/bin/fslroi ${rawdir}/Neg ${rawdir}/Negn 0 -1 0 -1 1 -1
+	${FSLDIR}/bin/fslroi ${rawdir}/Pos_b0 ${rawdir}/Pos_b0n 0 -1 0 -1 1 -1
+	${FSLDIR}/bin/fslroi ${rawdir}/Neg_b0 ${rawdir}/Neg_b0n 0 -1 0 -1 1 -1
+	${FSLDIR}/bin/imrm ${rawdir}/Pos
+	${FSLDIR}/bin/imrm ${rawdir}/Neg
+	${FSLDIR}/bin/imrm ${rawdir}/Pos_b0
+	${FSLDIR}/bin/imrm ${rawdir}/Neg_b0
+	${FSLDIR}/bin/immv ${rawdir}/Posn ${rawdir}/Pos
+	${FSLDIR}/bin/immv ${rawdir}/Negn ${rawdir}/Neg
+	${FSLDIR}/bin/immv ${rawdir}/Pos_b0n ${rawdir}/Pos_b0
+	${FSLDIR}/bin/immv ${rawdir}/Neg_b0n ${rawdir}/Neg_b0
+fi
+
+echo "Perform final merge"
+${FSLDIR}/bin/fslmerge -t ${rawdir}/Pos_Neg_b0 ${rawdir}/Pos_b0 ${rawdir}/Neg_b0
+${FSLDIR}/bin/fslmerge -t ${rawdir}/Pos_Neg ${rawdir}/Pos ${rawdir}/Neg
+paste ${rawdir}/Pos.bval ${rawdir}/Neg.bval >${rawdir}/Pos_Neg.bvals
+paste ${rawdir}/Pos.bvec ${rawdir}/Neg.bvec >${rawdir}/Pos_Neg.bvecs
+
+${FSLDIR}/bin/imrm ${rawdir}/Pos
+${FSLDIR}/bin/imrm ${rawdir}/Neg
+
+################################################################################################
+## Move files to appropriate directories
+################################################################################################
+echo "Move files to appropriate directories"
+mv ${rawdir}/extractedb0.txt ${topupdir}
+mv ${rawdir}/acqparams.txt ${topupdir}
+${FSLDIR}/bin/immv ${rawdir}/Pos_Neg_b0 ${topupdir}
+${FSLDIR}/bin/immv ${rawdir}/Pos_b0 ${topupdir}
+${FSLDIR}/bin/immv ${rawdir}/Neg_b0 ${topupdir}
+
+cp ${topupdir}/acqparams.txt ${eddydir}
+mv ${rawdir}/index.txt ${eddydir}
+mv ${rawdir}/series_index.txt ${eddydir}
+${FSLDIR}/bin/immv ${rawdir}/Pos_Neg ${eddydir}
+mv ${rawdir}/Pos_Neg.bvals ${eddydir}
+mv ${rawdir}/Pos_Neg.bvecs ${eddydir}
+mv ${rawdir}/Pos.bv?? ${eddydir}
+mv ${rawdir}/Neg.bv?? ${eddydir}
+
+echo -e "\n END: ${scriptName}"
diff --git a/DiffusionPreprocessing/scripts/eddy_postproc.sh b/DiffusionPreprocessing/scripts/eddy_postproc.sh
index 1c7893217..a88cfaeb7 100755
--- a/DiffusionPreprocessing/scripts/eddy_postproc.sh
+++ b/DiffusionPreprocessing/scripts/eddy_postproc.sh
@@ -3,13 +3,17 @@ source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; al
 echo -e "\n START: eddy_postproc"
 
 #Hard-Coded filename. Flag from eddy to indicate that the jac method has been used for resampling
-EddyJacFlag="JacobianResampling" 
+EddyJacFlag="JacobianResampling"
 
 workingdir=$1
-GdCoeffs=$2  #Coefficients for gradient nonlinearity distortion correction. If "NONE" this corrections is turned off
-CombineDataFlag=$3   #2 for including in the ouput all volumes uncombined (i.e. output file of eddy)
-                     #1 for including in the ouput and combine only volumes where both LR/RL (or AP/PA) pairs have been acquired
-                     #0 As 1, but also include uncombined single volumes"
+GdCoeffs=$2        #Coefficients for gradient nonlinearity distortion correction. If "NONE" this corrections is turned off
+CombineDataFlag=$3 #2 for including in the ouput all volumes uncombined (i.e. output file of eddy)
+                   #1 for including in the ouput and combine only volumes where both LR/RL (or AP/PA) pairs have been acquired
+                   #0 As 1, but also include uncombined single volumes"
+SelectBestB0=$4 #0 only the actual diffusion data was fed into eddy
+                #1 least distorted b0 was prepended to the eddy input
+                # Note: This numeric value is used within the script as a numeric that controls
+                # the number of volumes to skip, so it isn't just used as 0/1 "boolean".
 
 globalscriptsdir=${HCPPIPEDIR_Global}
 
@@ -31,70 +35,80 @@ qc_command+=(-v)
 "${qc_command[@]}"
 
 #Prepare for next eddy Release
-#if [ ! -e ${eddydir}/${EddyJacFlag} ]; then 
+#if [ ! -e ${eddydir}/${EddyJacFlag} ]; then
 #    echo "LSR resampling has been used. Eddy Output has already been combined."
 #    cp ${eddydir}/Pos.bval ${datadir}/bvals
 #    cp ${eddydir}/Pos.bvec ${datadir}/bvecs
 #    $FSLDIR/bin/imcp ${eddydir}/eddy_unwarped_images ${datadir}/data
 #else
+if [ ${SelectBestB0} -eq 1 ]; then
+	cut -d' ' -f2- ${eddydir}/Pos_Neg.bvals >${datadir}/bvals       # removes first value from bvals
+	cut -d' ' -f2- ${eddydir}/Pos_Neg.bvecs >${datadir}/bvecs_noRot # removes first value from bvecs
+fi
 if [ ${CombineDataFlag} -eq 2 ]; then
-	${FSLDIR}/bin/imcp  ${eddydir}/eddy_unwarped_images ${datadir}/data
-	cp ${eddydir}/Pos_Neg.bvals ${datadir}/bvals
-	cp ${datadir}/bvecs ${datadir}/bvecs_noRot
-	cp ${eddydir}/eddy_unwarped_images.eddy_rotated_bvecs ${datadir}/bvecs
+	# remove first volume as this is the reference b0, which was added to the dataset before running eddy
+	if [ ${SelectBestB0} -eq 1 ]; then
+		${FSLDIR}/bin/fslroi ${eddydir}/eddy_unwarped_images ${datadir}/data 1 -1
+		cut -d' ' -f2- ${eddydir}/eddy_unwarped_images.eddy_rotated_bvecs ${datadir}/bvecs # removes first value from bvecs
+	else
+		${FSLDIR}/bin/imcp ${eddydir}/eddy_unwarped_images ${datadir}/data
+		cp ${eddydir}/Pos_Neg.bvals ${datadir}/bvals
+		cp ${datadir}/bvecs ${datadir}/bvecs_noRot
+		cp ${eddydir}/eddy_unwarped_images.eddy_rotated_bvecs ${datadir}/bvecs
+	fi
 else
 	echo "JAC resampling has been used. Eddy Output is now combined."
-	PosVols=`wc ${eddydir}/Pos.bval | awk {'print $2'}`
-	NegVols=`wc ${eddydir}/Neg.bval | awk {'print $2'}`    #Split Pos and Neg Volumes
-	${FSLDIR}/bin/fslroi ${eddydir}/eddy_unwarped_images ${eddydir}/eddy_unwarped_Pos 0 ${PosVols}
-	${FSLDIR}/bin/fslroi ${eddydir}/eddy_unwarped_images ${eddydir}/eddy_unwarped_Neg ${PosVols} ${NegVols}
+	PosVols=$(wc ${eddydir}/Pos.bval | awk {'print $2'})
+	NegVols=$(wc ${eddydir}/Neg.bval | awk {'print $2'}) # Split Pos and Neg Volumes
+	${FSLDIR}/bin/fslroi ${eddydir}/eddy_unwarped_images ${eddydir}/eddy_unwarped_Pos ${SelectBestB0} ${PosVols} # ignore extra first volume if ${SelectBestB0} is 1
+	${FSLDIR}/bin/fslroi ${eddydir}/eddy_unwarped_images ${eddydir}/eddy_unwarped_Neg $((PosVols + ${SelectBestB0})) ${NegVols}
 	# Note: 'eddy_combine' is apparently hard-coded to use "data" as the output NIFTI file name
 	${FSLDIR}/bin/eddy_combine ${eddydir}/eddy_unwarped_Pos ${eddydir}/Pos.bval ${eddydir}/Pos.bvec ${eddydir}/Pos_SeriesVolNum.txt \
-             ${eddydir}/eddy_unwarped_Neg ${eddydir}/Neg.bval ${eddydir}/Neg.bvec ${eddydir}/Neg_SeriesVolNum.txt ${datadir} ${CombineDataFlag}
+		${eddydir}/eddy_unwarped_Neg ${eddydir}/Neg.bval ${eddydir}/Neg.bvec ${eddydir}/Neg_SeriesVolNum.txt ${datadir} ${CombineDataFlag}
 
 	${FSLDIR}/bin/imrm ${eddydir}/eddy_unwarped_Pos
 	${FSLDIR}/bin/imrm ${eddydir}/eddy_unwarped_Neg
-	cp ${datadir}/bvals ${datadir}/bvals_noRot
-	cp ${datadir}/bvecs ${datadir}/bvecs_noRot
-    
+	if [ ${SelectBestB0} -eq 0 ]; then
+		cp ${datadir}/bvecs ${datadir}/bvecs_noRot
+	fi
+
 	#rm ${eddydir}/Pos.bv*
 	#rm ${eddydir}/Neg.bv*
 
- 
 	# Divide Eddy-Rotated bvecs to Pos and Neg
-	line1=`awk 'NR==1 {print; exit}' ${eddydir}/eddy_unwarped_images.eddy_rotated_bvecs`
-	line2=`awk 'NR==2 {print; exit}' ${eddydir}/eddy_unwarped_images.eddy_rotated_bvecs`
-	line3=`awk 'NR==3 {print; exit}' ${eddydir}/eddy_unwarped_images.eddy_rotated_bvecs`   
+	line1=$(awk 'NR==1 {print; exit}' ${eddydir}/eddy_unwarped_images.eddy_rotated_bvecs)
+	line2=$(awk 'NR==2 {print; exit}' ${eddydir}/eddy_unwarped_images.eddy_rotated_bvecs)
+	line3=$(awk 'NR==3 {print; exit}' ${eddydir}/eddy_unwarped_images.eddy_rotated_bvecs)
 	Posline1=""
 	Posline2=""
 	Posline3=""
-	for ((i=1; i<=$PosVols; i++)); do
-	    Posline1="$Posline1 `echo $line1 | awk -v N=$i '{print $N}'`"
-	    Posline2="$Posline2 `echo $line2 | awk -v N=$i '{print $N}'`"
-	    Posline3="$Posline3 `echo $line3 | awk -v N=$i '{print $N}'`"
+	for ((i = $((SelectBestB0 + 1)); i <= $((PosVols + ${SelectBestB0})); i++)); do
+		Posline1="$Posline1 $(echo $line1 | awk -v N=$i '{print $N}')"
+		Posline2="$Posline2 $(echo $line2 | awk -v N=$i '{print $N}')"
+		Posline3="$Posline3 $(echo $line3 | awk -v N=$i '{print $N}')"
 	done
-	echo $Posline1 > ${eddydir}/Pos_rotated.bvec
-	echo $Posline2 >> ${eddydir}/Pos_rotated.bvec
-	echo $Posline3 >> ${eddydir}/Pos_rotated.bvec
-	
+	echo $Posline1 >${eddydir}/Pos_rotated.bvec
+	echo $Posline2 >>${eddydir}/Pos_rotated.bvec
+	echo $Posline3 >>${eddydir}/Pos_rotated.bvec
+
 	Negline1=""
 	Negline2=""
 	Negline3=""
-	Nstart=$((PosVols + 1 ))
-	Nend=$((PosVols + NegVols))
-	for  ((i=$Nstart; i<=$Nend; i++)); do
-	    Negline1="$Negline1 `echo $line1 | awk -v N=$i '{print $N}'`"
-	    Negline2="$Negline2 `echo $line2 | awk -v N=$i '{print $N}'`"
-	    Negline3="$Negline3 `echo $line3 | awk -v N=$i '{print $N}'`"
+	Nstart=$((PosVols + 1 + ${SelectBestB0}))
+	Nend=$((PosVols + NegVols + ${SelectBestB0}))
+	for ((i = $Nstart; i <= $Nend; i++)); do
+		Negline1="$Negline1 $(echo $line1 | awk -v N=$i '{print $N}')"
+		Negline2="$Negline2 $(echo $line2 | awk -v N=$i '{print $N}')"
+		Negline3="$Negline3 $(echo $line3 | awk -v N=$i '{print $N}')"
 	done
-	echo $Negline1 > ${eddydir}/Neg_rotated.bvec
-	echo $Negline2 >> ${eddydir}/Neg_rotated.bvec
-	echo $Negline3 >> ${eddydir}/Neg_rotated.bvec
-	
+	echo $Negline1 >${eddydir}/Neg_rotated.bvec
+	echo $Negline2 >>${eddydir}/Neg_rotated.bvec
+	echo $Negline3 >>${eddydir}/Neg_rotated.bvec
+
 	# Average Eddy-Rotated bvecs. Get for each direction the two b matrices, average those and then eigendecompose the average b-matrix to get the new bvec and bval.
 	# Also outputs an index file (1-based) with the indices of the input (Pos/Neg) volumes that have been retained in the output
 	${globalscriptsdir}/average_bvecs.py ${eddydir}/Pos.bval ${eddydir}/Pos_rotated.bvec ${eddydir}/Neg.bval ${eddydir}/Neg_rotated.bvec ${datadir}/avg_data ${eddydir}/Pos_SeriesVolNum.txt ${eddydir}/Neg_SeriesVolNum.txt
-	
+
 	mv ${datadir}/avg_data.bval ${datadir}/bvals
 	mv ${datadir}/avg_data.bvec ${datadir}/bvecs
 	rm -f ${datadir}/avg_data.bv??
@@ -107,44 +121,43 @@ imcp ${eddydir}/eddy_unwarped_images.eddy_cnr_maps ${datadir}/cnr_maps
 # 'eddy' can return negative values in some low signal locations, so use -abs for determining the fov mask
 ${FSLDIR}/bin/fslmaths ${datadir}/data -abs -Tmin -bin -fillh ${datadir}/fov_mask
 
-	 
-if [ ! $GdCoeffs = "NONE" ] ; then
-    echo "Correcting for gradient nonlinearities"
-	  # Note: data in the warped directory is eddy-current and suspectibility distortion corrected (via 'eddy'), but prior to gradient distortion correction
-	  # i.e., "data_posteddy_preGDC" would be another way to think of it
-	  warpedDir=${datadir}/warped
-    mkdir -p ${warpedDir}
-    ${FSLDIR}/bin/immv ${datadir}/data ${warpedDir}/data_warped
-    ${FSLDIR}/bin/immv ${datadir}/fov_mask ${warpedDir}/fov_mask_warped
-    ${FSLDIR}/bin/immv ${datadir}/cnr_maps ${warpedDir}/cnr_maps_warped
-
-    # Dilation outside of the field of view to minimise the effect of the hard field of view edge on the interpolation
-	  DiffRes=`${FSLDIR}/bin/fslval ${warpedDir}/data_warped pixdim1`
-    DilateDistance=`echo "$DiffRes * 4" | bc`  # Extrapolates the diffusion data up to 4 voxels outside of the FOV
-    ${CARET7DIR}/wb_command -volume-dilate ${warpedDir}/data_warped.nii.gz $DilateDistance NEAREST ${warpedDir}/data_dilated.nii.gz
-
-    # apply gradient distortion correction
-    ${globalscriptsdir}/GradientDistortionUnwarp.sh --workingdir="${datadir}" --coeffs="${GdCoeffs}" --in="${warpedDir}/data_dilated" --out="${datadir}/data" --owarp="${datadir}/fullWarp"
-    ${FSLDIR}/bin/immv ${datadir}/fullWarp ${warpedDir}
-    ${FSLDIR}/bin/immv ${datadir}/fullWarp_abs ${warpedDir}
-    ${FSLDIR}/bin/imrm ${warpedDir}/data_dilated
-
-    # Transform CNR maps
-    ${CARET7DIR}/wb_command -volume-dilate ${warpedDir}/cnr_maps_warped.nii.gz $DilateDistance NEAREST ${warpedDir}/cnr_maps_dilated.nii.gz
-    ${FSLDIR}/bin/applywarp --rel --interp=spline -i ${warpedDir}/cnr_maps_dilated -r ${warpedDir}/cnr_maps_dilated -w ${warpedDir}/fullWarp -o ${datadir}/cnr_maps
-    ${FSLDIR}/bin/imrm ${warpedDir}/cnr_maps_dilated
-
-    # Transform field of view mask (using conservative trilinear interpolation with high threshold)
-    ${FSLDIR}/bin/applywarp --rel --interp=trilinear -i ${warpedDir}/fov_mask_warped -r ${warpedDir}/fov_mask_warped -w ${warpedDir}/fullWarp -o ${datadir}/fov_mask
-    ${FSLDIR}/bin/fslmaths ${datadir}/fov_mask -thr 0.999 -bin ${datadir}/fov_mask
-
-    echo "Computing gradient coil tensor to correct for gradient nonlinearities"
-    ${FSLDIR}/bin/calc_grad_perc_dev --fullwarp=${warpedDir}/fullWarp -o ${datadir}/grad_dev
-    ${FSLDIR}/bin/fslmerge -t ${datadir}/grad_dev ${datadir}/grad_dev_x ${datadir}/grad_dev_y ${datadir}/grad_dev_z
-    ${FSLDIR}/bin/fslmaths ${datadir}/grad_dev -div 100 ${datadir}/grad_dev #Convert from % deviation to absolute
-    ${FSLDIR}/bin/imrm ${datadir}/grad_dev_?
-    ${FSLDIR}/bin/imrm ${datadir}/trilinear
-    ${FSLDIR}/bin/imrm ${warpedDir}/data_dilated_vol1
+if [ ! $GdCoeffs = "NONE" ]; then
+	echo "Correcting for gradient nonlinearities"
+	# Note: data in the warped directory is eddy-current and suspectibility distortion corrected (via 'eddy'), but prior to gradient distortion correction
+	# i.e., "data_posteddy_preGDC" would be another way to think of it
+	warpedDir=${datadir}/warped
+	mkdir -p ${warpedDir}
+	${FSLDIR}/bin/immv ${datadir}/data ${warpedDir}/data_warped
+	${FSLDIR}/bin/immv ${datadir}/fov_mask ${warpedDir}/fov_mask_warped
+	${FSLDIR}/bin/immv ${datadir}/cnr_maps ${warpedDir}/cnr_maps_warped
+
+	# Dilation outside of the field of view to minimise the effect of the hard field of view edge on the interpolation
+	DiffRes=$(${FSLDIR}/bin/fslval ${warpedDir}/data_warped pixdim1)
+	DilateDistance=$(echo "$DiffRes * 4" | bc) # Extrapolates the diffusion data up to 4 voxels outside of the FOV
+	${CARET7DIR}/wb_command -volume-dilate ${warpedDir}/data_warped.nii.gz $DilateDistance NEAREST ${warpedDir}/data_dilated.nii.gz
+
+	# apply gradient distortion correction
+	${globalscriptsdir}/GradientDistortionUnwarp.sh --workingdir="${datadir}" --coeffs="${GdCoeffs}" --in="${warpedDir}/data_dilated" --out="${datadir}/data" --owarp="${datadir}/fullWarp"
+	${FSLDIR}/bin/immv ${datadir}/fullWarp ${warpedDir}
+	${FSLDIR}/bin/immv ${datadir}/fullWarp_abs ${warpedDir}
+	${FSLDIR}/bin/imrm ${warpedDir}/data_dilated
+
+	# Transform CNR maps
+	${CARET7DIR}/wb_command -volume-dilate ${warpedDir}/cnr_maps_warped.nii.gz $DilateDistance NEAREST ${warpedDir}/cnr_maps_dilated.nii.gz
+	${FSLDIR}/bin/applywarp --rel --interp=spline -i ${warpedDir}/cnr_maps_dilated -r ${warpedDir}/cnr_maps_dilated -w ${warpedDir}/fullWarp -o ${datadir}/cnr_maps
+	${FSLDIR}/bin/imrm ${warpedDir}/cnr_maps_dilated
+
+	# Transform field of view mask (using conservative trilinear interpolation with high threshold)
+	${FSLDIR}/bin/applywarp --rel --interp=trilinear -i ${warpedDir}/fov_mask_warped -r ${warpedDir}/fov_mask_warped -w ${warpedDir}/fullWarp -o ${datadir}/fov_mask
+	${FSLDIR}/bin/fslmaths ${datadir}/fov_mask -thr 0.999 -bin ${datadir}/fov_mask
+
+	echo "Computing gradient coil tensor to correct for gradient nonlinearities"
+	${FSLDIR}/bin/calc_grad_perc_dev --fullwarp=${warpedDir}/fullWarp -o ${datadir}/grad_dev
+	${FSLDIR}/bin/fslmerge -t ${datadir}/grad_dev ${datadir}/grad_dev_x ${datadir}/grad_dev_y ${datadir}/grad_dev_z
+	${FSLDIR}/bin/fslmaths ${datadir}/grad_dev -div 100 ${datadir}/grad_dev #Convert from % deviation to absolute
+	${FSLDIR}/bin/imrm ${datadir}/grad_dev_?
+	${FSLDIR}/bin/imrm ${datadir}/trilinear
+	${FSLDIR}/bin/imrm ${warpedDir}/data_dilated_vol1
 fi
 
 # mask out any data outside the field of view
diff --git a/DiffusionPreprocessing/scripts/run_eddy.sh b/DiffusionPreprocessing/scripts/run_eddy.sh
index db3e94aa3..63e90cb2b 100755
--- a/DiffusionPreprocessing/scripts/run_eddy.sh
+++ b/DiffusionPreprocessing/scripts/run_eddy.sh
@@ -1,19 +1,19 @@
 #!/bin/bash
 #~ND~FORMAT~MARKDOWN~
 #~ND~START~
-# 
+#
 # # run_eddy.sh
-# 
+#
 # ## Copyright Notice
 #
 # Copyright (C) 2012-2019 The Human Connectome Project
-# 
+#
 # * Washington University in St. Louis
 # * University of Minnesota
 # * Oxford University
-# 
+#
 # ## Author(s)
-# 
+#
 # * Stamatios Sotiropoulos - Analysis Group, FMRIB Centre
 # * Saad Jbabdi - Analysis Group, FMRIB Center
 # * Jesper Andersson - Analysis Group, FMRIB Center
@@ -25,42 +25,39 @@
 # [Human Connectome Project][HCP] (HCP) Pipelines
 #
 # ## License
-# 
+#
 # See the [LICENSE](https://github.com/Washington-University/Pipelines/blob/master/LICENCE.md) file
-# 
+#
 # ## Description
-# 
+#
 # This script runs FSL's eddy command as part of the Human Connectome Project's
-# Diffusion Preprocessing 
-# 
+# Diffusion Preprocessing
+#
 # ## Prerequisite Installed Software
-# 
+#
 # * [FSL][FSL] - FMRIB's Software Library (version 5.0.7 or later)
-# 
+#
 #   FSL's environment setup script must also be sourced
-# 
+#
 # ## Prerequisite Environment Variables
-# 
+#
 # See output of usage function: e.g. <code>$ ./run_eddy.sh --help</code>
-# 
+#
 # <!-- References -->
-# 
+#
 # [HCP]: http://www.humanconnectome.org
 # [FSL]: http://fsl.fmrib.ox.ac.uk
-# 
+#
 #~ND~END~
 
-
 # Load Function Libraries
 source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib
 
-
 # --------------------------------------------------------------------------------
 #  Usage Description Function
 # --------------------------------------------------------------------------------
 
-show_usage()
-{
+show_usage() {
 	cat <<EOF
 
 Usage: ${g_script_name} PARAMETER...
@@ -143,7 +140,7 @@ EOF
 #  ${useGpuVersion}   - Set to "True" if user has requested an attempt to use
 #                       the GPU-enabled version of eddy
 #  ${workingdir}      - User specified working directory
-#  ${produceDetailedOutlierStats} 
+#  ${produceDetailedOutlierStats}
 #                     - Set to "True" if user has requested that the GPU-enabled version
 #                       of eddy produce detailed statistics about outliers after each iteration
 #  ${replaceOutliers} - Set to "True" if user has requested that the GPU-enabled version
@@ -161,10 +158,9 @@ EOF
 #  Support Functions
 # --------------------------------------------------------------------------------
 
-get_options()
-{
+get_options() {
 	local arguments=($@)
-	
+
 	# global output variables
 	useGpuVersion="False"
 	produceDetailedOutlierStats="False"
@@ -180,97 +176,96 @@ get_options()
 	unset ol_nstd_val
 	extra_eddy_args=""
 	g_cuda_version=""
-	
+
 	# parse arguments
 	local index=0
 	local numArgs=${#arguments[@]}
 	local argument
-	
-	while [ ${index} -lt ${numArgs} ]
-	do
+
+	while [ ${index} -lt ${numArgs} ]; do
 		argument=${arguments[index]}
-		
+
 		case ${argument} in
-			-h | --help)
-				show_usage
-				exit 0
-				;;
-			-g | --gpu)
-				useGpuVersion="True"
-				index=$(( index + 1 ))
-				;;
-			--wss)
-				produceDetailedOutlierStats="True"
-				index=$(( index + 1 ))
-				;;
-			--repol)
-				replaceOutliers="True"
-				index=$(( index + 1 ))
-				;;
-			-w | --workingdir)
-				workingdir=${arguments[$(( index + 1 ))]}
-				index=$(( index + 2 ))
-				;;
-			-w=* | --workingdir=*)
-				workingdir=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--nvoxhp=*)
-				nvoxhp=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--sep_offs_move)
-				sep_offs_move="True"
-				index=$(( index + 1 ))
-				;;
-			--rms)
-				rms="True"
-				index=$(( index + 1 ))
-				;;
-			--ff=*)
-				ff_val=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--dont_peas)
-				dont_peas="--dont_peas"
-				index=$(( index + 1 ))
-				;;
-			--fwhm=*)
-				fwhm_value=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--resamp=*)
-				resamp_value=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--ol_nstd=*)
-				ol_nstd_val=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			--extra-eddy-arg=*)
-				extra_eddy_arg=${argument#*=}
-				extra_eddy_args+=" ${extra_eddy_arg} "
-				index=$(( index + 1 ))
-				;;
-			--cuda-version=*)
-				g_cuda_version=${argument#*=}
-				index=$(( index + 1 ))
-				;;
-			*)
-				show_usage
-				echo "ERROR: Unrecognized Option: ${argument}"
-				exit 1
-				;;
+		-h | --help)
+			show_usage
+			exit 0
+			;;
+		-g | --gpu)
+			useGpuVersion="True"
+			index=$((index + 1))
+			;;
+		--wss)
+			produceDetailedOutlierStats="True"
+			index=$((index + 1))
+			;;
+		--repol)
+			replaceOutliers="True"
+			index=$((index + 1))
+			;;
+		-w | --workingdir)
+			workingdir=${arguments[$((index + 1))]}
+			index=$((index + 2))
+			;;
+		-w=* | --workingdir=*)
+			workingdir=${argument#*=}
+			index=$((index + 1))
+			;;
+		--nvoxhp=*)
+			nvoxhp=${argument#*=}
+			index=$((index + 1))
+			;;
+		--sep_offs_move)
+			sep_offs_move="True"
+			index=$((index + 1))
+			;;
+		--rms)
+			rms="True"
+			index=$((index + 1))
+			;;
+		--ff=*)
+			ff_val=${argument#*=}
+			index=$((index + 1))
+			;;
+		--dont_peas)
+			dont_peas="--dont_peas"
+			index=$((index + 1))
+			;;
+		--fwhm=*)
+			fwhm_value=${argument#*=}
+			index=$((index + 1))
+			;;
+		--resamp=*)
+			resamp_value=${argument#*=}
+			index=$((index + 1))
+			;;
+		--ol_nstd=*)
+			ol_nstd_val=${argument#*=}
+			index=$((index + 1))
+			;;
+		--extra-eddy-arg=*)
+			extra_eddy_arg=${argument#*=}
+			extra_eddy_args+=" ${extra_eddy_arg} "
+			index=$((index + 1))
+			;;
+		--cuda-version=*)
+			g_cuda_version=${argument#*=}
+			index=$((index + 1))
+			;;
+		*)
+			show_usage
+			echo "ERROR: Unrecognized Option: ${argument}"
+			exit 1
+			;;
 		esac
 	done
-	
+
 	# check required parameters
 	if [ -z ${workingdir} ]; then
 		show_usage
 		echo "  Error: <working-dir> not specified - Exiting without running eddy"
 		exit 1
 	fi
-	
+
 	# report options
 	echo "-- ${g_script_name}: Specified Command-Line Options - Start --"
 	echo "   workingdir: ${workingdir}"
@@ -302,8 +297,7 @@ get_options()
 #   g_stdEddy
 #   g_gpuEnabledEddy
 #
-determine_eddy_tools_for_supported_six_series()
-{
+determine_eddy_tools_for_supported_six_series() {
 	g_stdEddy="${FSLDIR}/bin/eddy_openmp"
 
 	if [ "${useGpuVersion}" = "True" ]; then
@@ -325,7 +319,7 @@ determine_eddy_tools_for_supported_six_series()
 				# They have an ${FSLDIR}/bin/eddy. So use it.
 				g_gpuEnabledEddy="${FSLDIR}/bin/eddy"
 			elif [ -e ${FSLDIR}/bin/eddy_cuda ]; then
-				# They have an ${FSLDIR}/bin/eddy_cuda. So use it. 
+				# They have an ${FSLDIR}/bin/eddy_cuda. So use it.
 				g_gpuEnabledEddy="${FSLDIR}/bin/eddy_cuda"
 			else
 				# If they have neither an FSLDIR/bin/eddy or FSLDIR/bin/eddy_cuda,
@@ -357,15 +351,14 @@ determine_eddy_tools_for_supported_six_series()
 #   g_stdEddy - path to the standard (non-GPU) version of eddy
 #   g_gpuEnabledEddy - path to GPU-enabled version of eddy
 #
-determine_eddy_tools_to_use()
-{
+determine_eddy_tools_to_use() {
 	local fsl_version_file
 	local fsl_version
 	local fsl_version_array
 	local fsl_primary_version
 	local fsl_secondary_version
 	local fsl_tertiary_version
-	
+
 	# get the current version of FSL in use
 	fsl_version_file="${FSLDIR}/etc/fslversion"
 
@@ -381,7 +374,7 @@ determine_eddy_tools_to_use()
 	# FSL X.Y.Z would have X as primary, Y as secondary, and Z as tertiary versions
 
 	fsl_version_array=(${fsl_version//./ })
-	
+
 	fsl_primary_version="${fsl_version_array[0]}"
 	fsl_primary_version=${fsl_primary_version//[!0-9]/}
 
@@ -390,14 +383,14 @@ determine_eddy_tools_to_use()
 
 	fsl_tertiary_version="${fsl_version_array[2]}"
 	fsl_tertiary_version=${fsl_tertiary_version//[!0-9]/}
-	
-	if [[ $(( ${fsl_primary_version} )) -lt 5 ]]; then
+
+	if [[ $((${fsl_primary_version})) -lt 5 ]]; then
 		# e.g. 4.x.x
 		log_Err_Abort "FSL 5.0.7 or greater is required."
 
-	elif [[ $(( ${fsl_primary_version} )) -eq 5 ]]; then
+	elif [[ $((${fsl_primary_version})) -eq 5 ]]; then
 		# e.g. 5.x.x
-		if [[ $(( ${fsl_secondary_version} )) -gt 0 ]]; then
+		if [[ $((${fsl_secondary_version})) -gt 0 ]]; then
 			# e.g. 5.1.x, 5.2.x, 5.3.x, etc.
 			# There aren't any 5.1.x, 5.2.x, 5.3.x, etc. versions that we know
 			# at the time this code was written. We don't expect any such
@@ -406,7 +399,7 @@ determine_eddy_tools_to_use()
 			log_Err_Abort "FSL version ${fsl_version} is currently unsupported"
 		else
 			# e.g. 5.0.x
-			if [[ $(( ${fsl_tertiary_version} )) -le 8 ]]; then
+			if [[ $((${fsl_tertiary_version})) -le 8 ]]; then
 				# 5.0.7 or 5.0.8
 				g_stdEddy="${FSLDIR}/bin/eddy"
 				g_gpuEnabledEddy="${FSLDIR}/bin/eddy.gpu"
@@ -423,11 +416,11 @@ determine_eddy_tools_to_use()
 			fi
 		fi
 
-	elif [[ $(( ${fsl_primary_version} )) -eq 6 ]]; then
+	elif [[ $((${fsl_primary_version})) -eq 6 ]]; then
 		# e.g. 6.x.x
-		if [[ $(( ${fsl_secondary_version} )) -eq 0 ]]; then
+		if [[ $((${fsl_secondary_version})) -eq 0 ]]; then
 			# e.g. 6.0.x
-			if [[ $(( ${fsl_tertiary_version} )) -eq 0 ]]; then
+			if [[ $((${fsl_tertiary_version})) -eq 0 ]]; then
 				# 6.0.0
 				log_Err_Abort "FSL version ${fsl_version} is currently unsupported"
 			else
@@ -446,11 +439,11 @@ determine_eddy_tools_to_use()
 			determine_eddy_tools_for_supported_six_series
 		fi
 
-	elif [[ $(( ${fsl_primary_version} )) -gt 6 ]]; then
+	elif [[ $((${fsl_primary_version})) -gt 6 ]]; then
 		# e.g. 7.x.x
 		# These versions do not exist that we know of at this writing.
 		# For now, we'll assume that they will work like the 6.0.1 version.
-		determine_eddy_tools_for_supported_six_series		
+		determine_eddy_tools_for_supported_six_series
 
 	else
 		# If we reach here, the primary version is:
@@ -466,27 +459,26 @@ determine_eddy_tools_to_use()
 	fi
 }
 
-# 
+#
 # Function Description
 #  Main processing of script
 #
-#  Gets user specified command line options, runs appropriate eddy 
+#  Gets user specified command line options, runs appropriate eddy
 #
-main()
-{
+main() {
 	# Get Command Line Options
 	#
 	# Global Variables Set:
 	#  See documentation for get_options function
 	get_options "$@"
-	
+
 	# Determine the eddy tools to use
 	determine_eddy_tools_to_use
 
 	local stdEddy="${g_stdEddy}"
 	local gpuEnabledEddy="${g_gpuEnabledEddy}"
-	
-	# Determine which eddy executable to use based upon whether 
+
+	# Determine which eddy executable to use based upon whether
 	# the user requested use of the GPU-enabled version of eddy
 	# and whether the requested version of eddy can be found.
 
@@ -507,9 +499,9 @@ main()
 			log_Err_Abort "Non-GPU-enabled version of eddy NOT found: ${stdEddy}"
 		fi
 	fi
-	
+
 	log_Msg "eddy executable command to use: ${eddyExec}"
-	
+
 	# Add option to eddy command for producing detailed outlier stats after each
 	# iteration if user has requested that option _and_ the GPU-enabled version
 	# of eddy is to be used.  Also add option to eddy command for replacing
@@ -522,28 +514,28 @@ main()
 	rmsOption=""
 	ff_valOption=""
 	ol_nstd_option=""
-	
+
 	if [ "${eddyExec}" = "${gpuEnabledEddy}" ]; then
 		if [ "${produceDetailedOutlierStats}" = "True" ]; then
 			outlierStatsOption="--wss"
 		fi
-		
+
 		if [ "${replaceOutliers}" = "True" ]; then
 			replaceOutliersOption="--repol"
 		fi
-		
+
 		if [ "${nvoxhp}" != "" ]; then
 			nvoxhpOption="--nvoxhp=${nvoxhp}"
 		fi
-		
+
 		if [ "${sep_offs_move}" = "True" ]; then
 			sep_offs_moveOption="--sep_offs_move"
 		fi
-		
+
 		if [ "${rms}" = "True" ]; then
 			rmsOption="--rms"
 		fi
-		
+
 		if [ "${ff_val}" != "" ]; then
 			ff_valOption="--ff=${ff_val}"
 		fi
@@ -554,7 +546,7 @@ main()
 			ol_nstd_option="--ol_nstd=${ol_nstd_val}"
 		fi
 	fi
-	
+
 	log_Msg "outlier statistics option: ${outlierStatsOption}"
 	log_Msg "replace outliers option: ${replaceOutliersOption}"
 	log_Msg "nvoxhp option: ${nvoxhpOption}"
@@ -564,11 +556,11 @@ main()
 	log_Msg "ol_nstd_option: ${ol_nstd_option}"
 
 	# Main processing - Run eddy
-	
-	topupdir=`dirname ${workingdir}`/topup
-	
+
+	topupdir=$(dirname ${workingdir})/topup
+
 	${FSLDIR}/bin/imcp ${topupdir}/nodif_brain_mask ${workingdir}/
-	
+
 	eddy_command="${eddyExec} "
 	eddy_command+="${outlierStatsOption} "
 	eddy_command+="${replaceOutliersOption} "
@@ -588,20 +580,20 @@ main()
 	eddy_command+="--out=${workingdir}/eddy_unwarped_images "
 	eddy_command+="--flm=quadratic "
 
-	if [ ! -z "${dont_peas}" ] ; then
+	if [ ! -z "${dont_peas}" ]; then
 		eddy_command+="--dont_peas "
 	fi
 
-	if [ ! -z "${resamp_value}" ] ; then
+	if [ ! -z "${resamp_value}" ]; then
 		eddy_command+="--resamp=${resamp_value} "
 	fi
 
-	if [ ! -z "${ol_nstd_option}" ] ; then
+	if [ ! -z "${ol_nstd_option}" ]; then
 		eddy_command+="${ol_nstd_option} "
 	fi
-	
-	if [ ! -z "${extra_eddy_args}" ] ; then
-		for extra_eddy_arg in ${extra_eddy_args} ; do
+
+	if [ ! -z "${extra_eddy_args}" ]; then
+		for extra_eddy_arg in ${extra_eddy_args}; do
 			eddy_command+=" ${extra_eddy_arg} "
 		done
 	fi
@@ -610,7 +602,7 @@ main()
 	log_Msg "${eddy_command}"
 	${eddy_command}
 	eddyReturnValue=$?
-	
+
 	log_Msg "Completed with return value: ${eddyReturnValue}"
 	exit ${eddyReturnValue}
 }
@@ -626,8 +618,8 @@ g_script_name=$(basename "${0}")
 
 # Allow script to return a Usage statement, before any other output
 if [ "$#" = "0" ]; then
-    show_usage
-    exit 1
+	show_usage
+	exit 1
 fi
 
 # Verify that HCPPIPEDIR Environment variable is set
@@ -637,8 +629,8 @@ if [ -z "${HCPPIPEDIR}" ]; then
 fi
 
 # Load function libraries
-source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@"         # Debugging functions; also sources log.shlib
-source ${HCPPIPEDIR}/global/scripts/opts.shlib                 # Command line option functions
+source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib
+source ${HCPPIPEDIR}/global/scripts/opts.shlib         # Command line option functions
 
 opts_ShowVersionIfRequested $@
 
@@ -656,4 +648,4 @@ log_Check_Env_Var FSLDIR
 #
 # Invoke the 'main' function to get things started
 #
-main $@
+main "$@"
diff --git a/DiffusionPreprocessing/scripts/run_topup.sh b/DiffusionPreprocessing/scripts/run_topup.sh
index 028852684..4fb0cb2d1 100755
--- a/DiffusionPreprocessing/scripts/run_topup.sh
+++ b/DiffusionPreprocessing/scripts/run_topup.sh
@@ -11,7 +11,7 @@ topup_config_file=${configdir}/b02b0.cnf
 
 ${FSLDIR}/bin/topup --imain=${workingdir}/Pos_Neg_b0 --datain=${workingdir}/acqparams.txt --config=${topup_config_file} --out=${workingdir}/topup_Pos_Neg_b0 -v --fout=${workingdir}/topup_Pos_Neg_b0_field.nii.gz
 
-dimt=`${FSLDIR}/bin/fslval ${workingdir}/Pos_b0 dim4`
+dimt=$(${FSLDIR}/bin/fslval ${workingdir}/Pos_b0 dim4)
 dimt=$((${dimt} + 1))
 
 echo "Applying topup to get a hifi b0"
@@ -20,9 +20,9 @@ ${FSLDIR}/bin/fslroi ${workingdir}/Neg_b0 ${workingdir}/Neg_b01 0 1
 ${FSLDIR}/bin/applytopup --imain=${workingdir}/Pos_b01,${workingdir}/Neg_b01 --topup=${workingdir}/topup_Pos_Neg_b0 --datain=${workingdir}/acqparams.txt --inindex=1,${dimt} --out=${workingdir}/hifib0
 
 if [ ! -f ${workingdir}/hifib0.nii.gz ]; then
-    echo "run_topup.sh -- ERROR -- ${FSLDIR}/bin/applytopup failed to generate ${workingdir}/hifib0.nii.gz"
-    # Need to add mechanism whereby scripts that invoke this script (run_topup.sh)
-    # check for a return code to determine success or failure
+	echo "run_topup.sh -- ERROR -- ${FSLDIR}/bin/applytopup failed to generate ${workingdir}/hifib0.nii.gz"
+	# Need to add mechanism whereby scripts that invoke this script (run_topup.sh)
+	# check for a return code to determine success or failure
 fi
 
 ${FSLDIR}/bin/imrm ${workingdir}/Pos_b0*
@@ -32,10 +32,9 @@ echo "Running BET on the hifi b0"
 ${FSLDIR}/bin/bet ${workingdir}/hifib0 ${workingdir}/nodif_brain -m -f 0.2
 
 if [ ! -f ${workingdir}/nodif_brain.nii.gz ]; then
-    echo "run_topup.sh -- ERROR -- ${FSLDIR}/bin/bet failed to generate ${workingdir}/nodif_brain.nii.gz"
-    # Need to add mechanism whereby scripts that invoke this script (run_topup.sh)
-    # check for a return code to determine success or failure
+	echo "run_topup.sh -- ERROR -- ${FSLDIR}/bin/bet failed to generate ${workingdir}/nodif_brain.nii.gz"
+	# Need to add mechanism whereby scripts that invoke this script (run_topup.sh)
+	# check for a return code to determine success or failure
 fi
 
 echo -e "\n END: run_topup"
-
diff --git a/global/config/best_b0.cnf b/global/config/best_b0.cnf
new file mode 100755
index 000000000..ff6560a25
--- /dev/null
+++ b/global/config/best_b0.cnf
@@ -0,0 +1,33 @@
+# This configuration file was originally written by Matteo Bastiani and ? as part of the developing HCP pipeline
+# The original configuration file was released under the Apache license 2.0 (https://git.fmrib.ox.ac.uk/matteob/dHCP_neo_dMRI_pipeline_release/blob/master/LICENSE)
+
+# Resolution (knot-spacing) of warps in mm
+--warpres=20,14
+# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel)
+--subsamp=2,2
+# FWHM of gaussian smoothing
+--fwhm=8,4
+# Maximum number of iterations
+--miter=5,5
+# Relative weight of regularisation
+--lambda=0.0005,0.00001
+# If set to 1 lambda is multiplied by the current average squared difference
+--ssqlambda=1
+# Regularisation model
+--regmod=bending_energy
+# If set to 1 movements are estimated along with the field
+--estmov=1
+# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient
+--minmet=0
+# Quadratic or cubic splines
+--splineorder=3
+# Precision for calculation and storage of Hessian
+--numprec=double
+# Linear or spline interpolation
+--interp=spline
+# If set to 1 the images are individually scaled to a common mean intensity 
+--scale=1
+# To ensure that field is written
+--fout=field_test
+# To ensure that "unwarped" images are written
+--iout=iout_test
diff --git a/global/scripts/select_dwi_vols b/global/scripts/select_dwi_vols
new file mode 100755
index 000000000..aa3d933fb
--- /dev/null
+++ b/global/scripts/select_dwi_vols
@@ -0,0 +1,82 @@
+#!/bin/bash
+
+if [ "$4" == "" ];then
+    echo ""
+    echo "Usage: select_dwi_vols <data> <bvals> <output> <approx_bval> [other options]"
+    echo ""
+    echo " Optional arguments:"
+    echo ""
+    echo " -b <value>        additional bvalues to be extracted (can be used multiple times)"
+    echo " -m                output mean instead of concat"
+    echo " -v                output variance instead of concat"
+    echo " -obv <bvecs>      produce <output>.bval and <output>.bvec"
+    echo " -db <delta_bvals> only include volumes with offset in b-value less than"
+    echo "                   this value in either direction (default: 100 s/mm2)"
+    echo ""
+    exit 1
+fi
+
+# parse compulsory arguments
+d=$1
+bvals=`cat $2`
+o=$3
+b=$4
+shift
+shift
+shift
+# parse optional arguments
+outmean=""
+outvar=""
+outbv=""
+delta_bvals="100"
+while [ ! -z "$1" ]
+do
+  case "$1" in
+      -m) outmean="-m";;
+      -v) outvar="-v";;
+      -obv) outbv=$2;shift;;
+      -b) b="$b $2";shift;;
+      -db) delta_bvals="$2";shift;;
+  esac
+  shift
+done
+
+if [ "$outmean" == "-m" ] && [ "$outvar" == "-v" ];then
+    echo "Error: you can either use -m or -v but not both"
+    exit 1
+fi
+
+list=""
+cmd="awk '{"
+echo "Extracting bvals=$b"
+for bval in $b;do
+    cnt=0
+    for i in $bvals;do
+	j=`echo $i | awk -F"E" 'BEGIN{OFMT="%10.10f"} {print $1 * (10 ^ $2)}' `
+	j=${j/.*}
+	j=`echo "$j - $bval" |bc | awk ' { if($1>=0) { print $1} else {print $1*-1 }}'`
+	if [ $j -lt ${delta_bvals} ];then
+	    if [ "${list}" == "" ];then
+		list="${cnt}"
+        cmd="$cmd print \$$(($cnt+1)) \" \""
+	    else
+		list="${list},${cnt}"
+		cmd="$cmd \" \" \$$(($cnt+1))"
+	    fi
+	fi
+	cnt=$(($cnt + 1))
+    done
+done
+echo volume list = $list
+cmd="$cmd }'"
+
+
+if [ "$outbv" != "" ];then
+    bcmd="echo $bvals | `echo $cmd` > ${o}.bval"
+    eval $bcmd
+    vcmd="cat $outbv | `echo $cmd` > ${o}.bvec"
+    eval $vcmd
+fi
+
+echo $FSLDIR/bin/fslselectvols -i $d -o $o --vols=$list $outmean $outvar
+$FSLDIR/bin/fslselectvols -i $d -o $o --vols=$list $outmean $outvar