lnifmri_prep03_sdc

	#!/bin/bash
#
# LNiF MRI Lab Pipeline for preprocessing of
# structural and functional MRI dataset
# stefano.tambalo@unitn.it - 20200310
# rev. 20210513
#
# Correct Susceptibility Induced Distortion:
# -   Fieldmap (if present)
# -   AP/PA rev. ph. encoding (if present)
# -   Nonlin. reg to T1
#
########################################################################
export LC_ALL=C

shopt -s extglob
Usage() {
    printf "\n***************************************************\n"
    printf "LNiF MRI Lab Preprocessing Pipeline\n"
    printf "CIMeC - Center for Mind/Brain Science\n\n"
    printf "\tFSLDIR not set!\n\n"
    printf "\tRequires:\n\t\tFSL (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSL)\n"
    printf "\n\n(c) stefano.tambalo@unitn.it, 2020.\n"
    printf "***************************************************\n\n"
    exit 1
}
[ "$FSLDIR" = "" ] && Usage

Start() {
	stamp=`date +%Y%m%d_%H%M%S`
    printf "\n***************************************************\n"
    printf "LNiF MRI Lab Preprocessing Pipeline\n"
    printf "CIMeC - Center for Mind/Brain Science\n\n"
    printf "Process functional data (suscept. dist. correction)\n"    

}

NextStep() {
    printf "\n***************************************************\n"
    printf "LNiF MRI Lab Preprocessing Pipeline\n"
    printf "CIMeC - Center for Mind/Brain Science\n\n"
    printf "Complete. Exec. Time: ${exectime}\nNow run lnifmri_prep04_filter"
    printf "\n\n***************************************************\n\n"
    exit 1
}

# 20210513 - This is disabled. See heredoc in the topup SDC section
# #Formatted configuration file with default params to run TOPUP
# arg="# ##############################################################\n\
# # Configuration file for AP/PA unwarping of ep2d_bold timeseries\n\
# # LNIF MRI LAB @ UniTN - stefano.tambalo@unitn.it\n\
# # ##############################################################\n\
# # Resolution (knot-spacing) of warps in mm\n\
# --warpres=6,4,1\n\
# # Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel)\n\
# --subsamp=1,1,1\n\
# # FWHM of gaussian smoothing\n\
# --fwhm=4,2,0\n\
# # Maximum number of iterations\n\
# --miter=5,10,20\n\
# # Relative weight of regularisation\n\
# --lambda=0.000015,0.0000000005,0.00000000001\n\
# # If set to 1 lambda is multiplied by the crrent average squared difference\n\
# --ssqlambda=1\n\
# # Regularisation model\n\
# --regmod=bending_energy\n\
# # If set to 1 movements are estimated along with the field\n\
# --estmov=1,1,0\n\
# # 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient\n\
# --minmet=0,0,1\n\
# # Quadratic or cubic splines\n\
# --splineorder=3\n\
# # Precision for calculation and storage of Hessian\n\
# --numprec=double\n\
# # Linear or spline interpolation\n\
# --interp=spline\n\
# # If set to 1 the images are individually scaled to a common mean intensity\n\
# --scale=1\n";

exectime=0
SECONDS=0
QCDIR=`pwd`/qc
LOGDIR=`pwd`/logs

Start

#Wait for 01_anat and 02_preproc to finish
while true; do
	if [[ -f ${LOGDIR}/lnifmri_prep01_greenlight.txt  && -f ${LOGDIR}/lnifmri_prep02_greenlight.txt ]]; then
		stamp=`date +%Y%m%d_%H%M%S`
		printf "\t%s: Greenlight from 01_anat and 02_preproc: starting SDC process\n\n" "${stamp}"
		break
  	fi

	stamp=`date +%Y%m%d_%H%M%S`
	printf "\t%s: Waiting greenlight from 01_anat and 02_preproc\n" "${stamp}"
	sleep 180

done

for epiimg in `cat ${LOGDIR}/lnifmri_prep_inputlist.txt` ; do
    epiname=$($FSLDIR/bin/remove_ext ${epiimg})

    epifilename=$(basename ${epiname})

    readarray -d / -t fileparts <<< "${epiimg}"
    nid=$((${#fileparts[@]}-1))
    sid=$((${#fileparts[@]}-3))
    nameID=$(echo ${fileparts[${nid}]})
    subID=$(echo ${fileparts[${sid}]})
    nameIDnoext=$($FSLDIR/bin/remove_ext ${nameID})

	#define different images/masks to load
	printf -v STR "${subID}_ses-01_T1w"
	t1struc=${subID}/$STR.anat/T1_biascorr.nii.gz
	bstruc=${subID}/$STR.anat/T1_biascorr_brain.nii.gz
	wmstruc=${subID}/$STR.anat/T1_wmedge.nii.gz
	csfstruc=${subID}/$STR.anat/T1_csfedge.nii.gz
	brainmask=${subID}/$STR.anat/T1_biascorr_brain_mask.nii.gz
	lesionmask=${subID}/$STR.anat/lesionmask.nii.gz
	lesionmaskinv=${subID}/$STR.anat/lesionmaskinv.nii.gz
	std=$FSLDIR/data/standard/MNI152_T1_2mm.nii.gz
	
	printf "Scan Name: ${subID} - ${nameID}\n"

	#Extract core parameters from JSON
	manuf=$(awk '$1 ~ /"Manufacturer"/ {print $2}' ${epiname}_parameters.json | tr -d ',')
	case "${manuf}" in
		"\"Siemens\"")
			phenctag="PhaseEncodingDirection"
			echosp=$(awk '$1 ~ /"DerivedVendorReportedEchoSpacing"/ {print $2}' ${epiname}_parameters.json | tr -d ',')
		;;
		"\"Philips\"")
			phenctag="PhaseEncodingAxis"
			# % Calculate Echo Spacing parameter for Philips-derived dataset
			# %
			# % based on: % calc_echo_spacing_philips_mod.m
			# % Formula from http://www.spinozacentre.nl/wiki/index.php/NeuroWiki:Current_developments#B0_correction
			# %
			# % stefano.tambalo@unitn.it - 20210604
			# % *********************************************************************
			# % variables, please adapt for each dataset (can be found in *.PAR file)
			# % matrix size in y direction
			matrixsize_phase_enc_dir=$(awk '$1 ~ /"AcquisitionMatrixPE"/ {print $2}' ${epiname}_parameters.json | tr -d ',')
			# % EPI factor
			epifactor=$(awk '$1 ~ /"EchoTrainLength"/ {print $2}' ${epiname}_parameters.json | tr -d ',')
			# % water-fat shift per pixel
			water_fat_shift_pixel=10.604
			# % magnetic field strength (T)
			fieldstrength_tesla=$(awk '$1 ~ /"MagneticFieldStrength"/ {print $2}' ${epiname}_parameters.json | tr -d ',')
			# % *********************************************************************
			# % Do not edit below this line
			water_fat_diff_ppm=3.35
			# % gyromagnetic ratio for proton (1H)
			resonance_freq_mhz_tesla=42.576
			echo_train_length=$(echo "${epifactor}+1" | bc -l)
			# % water_fat_shift_hz 3T = 434.215 Hz
			water_fat_shift_hz=$(echo "scale=4;${fieldstrength_tesla}*${water_fat_diff_ppm}*${resonance_freq_mhz_tesla}" | bc -l)
			BW_hz_pixel=$(echo "scale=4;${water_fat_shift_hz}/${water_fat_shift_pixel}" | bc -l)
			totBW=$(echo "scale=4;${BW_hz_pixel}*${echo_train_length}" | bc -l)
			echo_spacing_sec=$(echo "scale=4;1/${totBW}" | bc -l)
			echo_spacing_msec=$(echo "scale=4;${echo_spacing_sec}*1000" | bc -l)
			echosp=${echo_spacing_sec}
		;;
	esac

	acqmtxpe=$(awk '$1 ~ /"AcquisitionMatrixPE"/ {print $2}' ${epiname}_parameters.json | tr -d ',')
	totreadtime=`echo "scale=4;${echosp}*(${acqmtxpe}-1)" | bc -l`
    phenc=$(awk '$1 ~ /"'"${phenctag}"'"/ {print $2}' ${epiname}_parameters.json | tr -d ',')

	printf "Manufacturer: %s\n" "${manuf}"
	printf "Echo Spacing: %s\n" "${echosp}"

	if [ "${phenc}" = "\"j-\"" ];
		then
			pedir="-y"
		else
			pedir="y"
	fi

	# Run SDC preparation only once per subID
	if [ ! -e ${subID}/sdc/lnifmri_prep03_greenlight.txt ]; then
		# #################################################################################################
		# PREPARE FIELDMAP
		if 	[ -e ${subID}/sdc/*gre_field_mapping*e1.nii.gz -a \
			  -e ${subID}/sdc/*gre_field_mapping*e2.nii.gz -a \
			  -e ${subID}/sdc/*gre_field_mapping*ph.nii.gz ]; then
			printf "\tPreparing fielmap\n"
			runfmap=1
			#Set array to store Echo Times
			declare ediff

			#Process Magnitude GRE dataset
			for n in `ls ${subID}/sdc/*gre_field_mapping*_e[12].nii.gz` ; do
			    g=`$FSLDIR/bin/remove_ext ${n}`

				# printf "\tReformat JSON file\n"
	    	    cat ${g}.json | tr -d '\n' > ${g}_temp1.json
	       		sed 's/,\t\"/,\n\"/g' ${g}_temp1.json > ${g}_temp2.json
	       	 	cat ${g}_temp2.json | tr -d '\t' > ${g}_parameters.json

			    # printf "\tRunning BET on ${n}\n"
			    $FSLDIR/bin/bet $n ${g}_brain \
			    				-f 0.4 \
			    				-g 0.1 \
			    				-R;
			    
			    #Erode brain-extracted image to exclude noisy voxels
			    $FSLDIR/bin/fslmaths ${g}_brain -ero ${g}_brain_eroded

    			# Add box kernel?
				# see: https://doi.org/10.1371/journal.pone.0185647, page 5
			    # $FSLDIR/bin/fslmaths ${g}_brain -kernel box 5 -ero ${g}_brain_eroded
				
				#Parse EchoTime from JSON
				etime=$(awk '$1 ~ /"EchoTime"/ {print $2}' ${g}_parameters.json | tr -d ',')
			    ediff+=(${etime})
			done

			#Set arguments for fsl_prepare_fieldmap
			phaimg=`ls ${subID}/sdc/*gre_field_mapping*ph.nii.gz`

			#Consider 1st echo of magnitude image (higher SNR)
			magbrain=`ls ${subID}/sdc/*gre_field_mapping*e1_brain_eroded.nii.gz`
			magimg=`ls ${subID}/sdc/*gre_field_mapping*e1.nii.gz`
			fmap="${subID}/sdc/ep2d_fmap_fieldmap"

			#Compute echo time difference
			echdif=`echo "(${ediff[1]} - ${ediff[0]})*1000" | bc -l`

			#Derive Fieldmap
			printf "\tRunning fieldmap estimation with EchoTime diff: ${echdif}\n\n"
			$FSLDIR/bin/fsl_prepare_fieldmap SIEMENS $phaimg $magbrain $fmap $echdif
			
			# Add gaussian smoothing?
			# see: https://doi.org/10.1371/journal.pone.0185647, page 5
			# $FSLDIR/bin/fslmaths ${fmap} -s 5 ${fmap}

		else
			runfmap=0
			printf "\n\tFieldmap data incomplete or not available.\n"
		fi

		#################################################################################################
		# PREPARE AP/PA FOR TOPUP
		if 	[ -e ${subID}/sdc/*ep2d_se_distcorr_AP_*.nii.gz -a -e ${subID}/sdc/*ep2d_se_distcorr_PA_*.nii.gz ]; then
			printf "\n\tPreparing config file for topup\n"
			runtopup=1
			#Default parameters config-file for topup
			cnf="${subID}/sdc/ep2d_topup_cfgpars.cnf"
			cnfname=`basename ${cnf} .cnf`
			#Print config-file - Updated 20210513
			cat << EOF > ${cnf}
# ##############################################################
# Configuration file for AP/PA unwarping of ep2d_bold timeseries
# LNIF MRI LAB @ UniTN - stefano.tambalo@unitn.it
# ##############################################################
# Resolution (knot-spacing) of warps in mm
--warpres=6,4,1
# Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel)
--subsamp=1,1,1
# FWHM of gaussian smoothing
--fwhm=4,2,0
# Maximum number of iterations
--miter=5,10,20
# Relative weight of regularisation
--lambda=0.000015,0.0000000005,0.00000000001
# 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,1,0
# 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient
--minmet=0,0,1
# 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
EOF
			#Prepare config file
			TUPDIR=${subID}/func/sdc_tup
			mkdir -m 777 ${TUPDIR}
			#Write config-file for topup
			printf -v par "%s/%s_blipudpars.par" "${TUPDIR}" "${nameIDnoext}"
			>${TUPDIR}/${nameIDnoext}_blipudpars.par
			#>${par}
						
			printf "\tMerging opposite phase encoding SE-EPI\n\n"
			ls -d1 ${subID}/sdc/*ep2d_se_distcorr_AP*.nii.gz >> ${subID}/sdc/lnifmri_prep03_phenctomerge.txt
			ls -d1 ${subID}/sdc/*ep2d_se_distcorr_PA*.nii.gz >> ${subID}/sdc/lnifmri_prep03_phenctomerge.txt

			phencdir=0
			for p in `cat ${subID}/sdc/lnifmri_prep03_phenctomerge.txt` ; do
				#Collapse volumes in case of 4D acquisition
				printf -v phencimg "ep2d_phenc_${phencdir}.nii.gz" "$(( ++phencdir ))"
				$FSLDIR/bin/fslmaths ${p} -Tmean ${subID}/sdc/${phencimg}

				#Parsing of AP/PA JSON files to properly compose config file for TOPUP
				parsphenc=${p%.nii.gz}
				#Alternate phase encoding
			    altphenc=$(awk '$1 ~ /"'"${phenctag}"'"/ {print $2}' ${parsphenc}.json | tr -d ',')
				case "${manuf}" in
					"\"Siemens\"")
				    case "${altphenc}" in
						"\"j-\"")
							# AP - append to config file
							printf "%s %s\n" "0 -1 0" "${totreadtime}" >> ${par}
							;;
						"\"j\"")
							# PA - append to config file
							printf "%s %s\n" "0 1 0" "${totreadtime}" >> ${par}
							;;
						"\"i-\"")
							# LR - append to config file
							printf "%s %s\n" "-1 0 0" "${totreadtime}" >> ${par}
							;;
						"\"i\"")
							# RL - append to config file
							printf "%s %s\n" "1 0 0" "${totreadtime}" >> ${par}
							;;
					esac
					;;
					"\"Philips\"")
					# Here we don't have info on the PE dir, we assume that alphabetical order of
					# SDC scans is correct.
					# stefano.tambalo@unitn.it - 20210604
					#
					# See: https://github.com/nipreps/sdcflows/issues/5#issuecomment-586417652
					# Note that until Philips specifies the phase encoding polarity,
					# dcm2niix will only be able to populate PhaseEncodingAxis (e.g. i or j)
					# rather than PhaseEncodingDirection (disambiguating j from j-).
				    case "${altphenc}" in
						"\"j\"")
							# AP + PA - overwrite config file
							printf "%s %s\n%s %s" "0 -1 0" "${totreadtime}" "0 1 0" "${totreadtime}" > ${par}
							;;
						"\"i\"")
							# LR + RL - overwrite config file
							printf "%s %s\n%s %s" "-1 0 0" "${totreadtime}" "1 0 0" "${totreadtime}" > ${par}
							;;
					esac
					;;
				esac					
			done

			$FSLDIR/bin/fslmerge -t ${subID}/sdc/ep2d_topup_fieldmap `ls -d1 ${subID}/sdc/*ep2d_phenc*.nii.gz`
		
		else
			runtopup=0
			printf "\n\tAP/PA SE data incomplete or not available.\n\n"
		fi

		# Mark SDC preparation as complete
		stamp=`date +%Y%m%d_%H%M%S`
		printf "LNiF MRI Lab Preprocessing Pipeline\nDo not remove this file!\nCreated: ${stamp}" > ${subID}/sdc/lnifmri_prep03_greenlight.txt
	fi

	#################################################################################################
	# UNWARPING OF EPI
	#Load bold EPI scan
	epivol="${epiname}_onevol"
	$FSLDIR/bin/fslroi ${epiimg} ${epivol} 3 1
	
	# printf "\tScan Name: ${subID} - ${nameID}\n"
		
	##########################################################################################
	#
	#	Fieldmap method
	#
	##########################################################################################

	if [[ "${runfmap}" -eq "0" ]]; then
		printf "\n\tSkipping Fieldmap unwarping\n";
	else
		[ ! -d ${subID}/func/sdc_fmp ] && { mkdir -m 777 ${subID}/func/sdc_fmp; }
		FMAPDIR=${subID}/func/sdc_fmp
		printf "\n\tSaving fieldmap unwarped data to: $FMAPDIR\n"

		#Call epi_reg w/fieldmap est. (Single step)
		printf "\tRunning epi_reg\n\n"
		$LNIFMRIDIR/lnifmri_util_epireg	--echospacing=${echosp} \
							--wmseg=${wmstruc} \
							--fmap=${fmap} \
							--fmapmag=${magimg} \
							--fmapmagbrain=${magbrain} \
							--pedir=${pedir} \
							--epi=${epivol} \
							--t1=${t1struc} \
							--t1brain=${bstruc} \
							--out=$FMAPDIR/${nameIDnoext} \
							-v --noclean \
							--weight=${lesionmaskinv};

		printf "\n\tWarping EPI to T1\n"
		$FSLDIR/bin/applywarp	-i ${epiimg} \
								-r ${t1struc} \
								-w $FMAPDIR/${nameIDnoext}_warp.nii.gz \
								--interp=spline \
								-o $FMAPDIR/${nameIDnoext};
		$FSLDIR/bin/fslroi $FMAPDIR/${nameIDnoext} ${subID}/func/${epifilename}_fmap_checkwarpedtostruc 3 1
		#$FSLDIR/bin/fslroi $FMAPDIR/${nameIDnoext} ${epiname}_fmap_checkwarpedtostruc 3 1
		#printf "${epiname}_fmap_checkwarpedtostruc.nii.gz ${t1struc}\n" >> $QCDIR/qc_filelist_check2struct.txt
		printf "${subID}/func/${epifilename}_fmap_checkwarpedtostruc.nii.gz ${t1struc}\n" >> $QCDIR/qc_filelist_check2struct.txt

		##########################################################################################
		#
		#	EPI to MNI
		#
		##########################################################################################	
		#Warping 4D EPI in standard (2mm MNI) space
		printf "\tCreating EPI to MNI composite warps\n"
		$FSLDIR/bin/convertwarp	-r ${std} \
								--warp1=$FMAPDIR/${nameIDnoext}_warp \
								--warp2=${subID}/${STR}.anat/T1_to_MNI_nonlin_field.nii.gz \
								-j $FMAPDIR/${nameIDnoext}_jac2MNI.nii.gz \
								-o $FMAPDIR/${nameIDnoext}_warp2MNI.nii.gz;

		printf "\tWarping EPI to MNI\n"
		$FSLDIR/bin/applywarp	-i ${epivol} \
								-r ${std} \
								-w $FMAPDIR/${nameIDnoext}_warp2MNI.nii.gz \
								--interp=spline \
								-o ${epiname}_fmap_checkwarpedtostd;
		#$FSLDIR/bin/fslroi $FMAPDIR/${nameIDnoext}_wrp2std ${epiname}_fmap_checkwarpedtostd 30 1

		printf "\tFieldmap SDC DONE.\n"
	fi

	##########################################################################################
	#
	#	AP/PA method
	#
	##########################################################################################

	if [[ "${runtopup}" -eq "0" ]]; then
		printf "\n\tSkipping Topup unwarping\n";
	else
		[ ! -d ${subID}/func/sdc_tup ] && { mkdir -m 777 ${subID}/func/sdc_tup; }
		# TUPDIR=${subID}/func/sdc_tup
		printf "\n\tSaving TOPUP unwarped data to: $TUPDIR\n"

		printf "\tRunning topup\n"
		$FSLDIR/bin/topup	--imain=${subID}/sdc/ep2d_topup_fieldmap.nii.gz \
							--datain=${par} \
							--config=${cnf} \
							--out=$TUPDIR/${nameIDnoext} \
							--iout=$TUPDIR/${nameIDnoext}_unw \
							--fout=$TUPDIR/${nameIDnoext}_field;
		#Hz to rad/s
		printf "\tConverting FMAP values: Hz -> rad/s\n"
		$FSLDIR/bin/fslmaths $TUPDIR/${nameIDnoext}_field -mul 6.2832 $TUPDIR/${nameIDnoext}_field_rads
	    $FSLDIR/bin/fslroi $TUPDIR/${nameIDnoext}_unw $TUPDIR/${nameIDnoext}_unw_onevol 0 1
	    $FSLDIR/bin/bet	$TUPDIR/${nameIDnoext}_unw_onevol $TUPDIR/${nameIDnoext}_unw_onevol_brain \
	    				-f 0.4 \
	    				-g 0.1 \
	    				-R;
	    #Erode brain-extracted image to exclude noisy voxels
	    $FSLDIR/bin/fslmaths $TUPDIR/${nameIDnoext}_unw_onevol_brain -ero $TUPDIR/${nameIDnoext}_unw_onevol_brain_eroded

		printf "\tRunning epi_reg with topup fieldmap\n\n"
		$LNIFMRIDIR/lnifmri_util_epireg	--epi=${epivol} \
							--t1=${t1struc} \
							--t1brain=${bstruc} \
							--wmseg=${wmstruc} \
							--fmap=$TUPDIR/${nameIDnoext}_field_rads \
							--fmapmag=$TUPDIR/${nameIDnoext}_unw_onevol \
							--fmapmagbrain=$TUPDIR/${nameIDnoext}_unw_onevol_brain_eroded \
							--echospacing=${echosp} \
							--pedir=${pedir} \
							--out=$TUPDIR/${nameIDnoext} \
							-v --noclean \
							--weight=${lesionmaskinv};

		printf "\n\tWarping EPI to T1\n"
		$FSLDIR/bin/applywarp	-i ${epiimg} \
								-r ${t1struc} \
								-w $TUPDIR/${nameIDnoext}_warp.nii.gz \
								--interp=spline \
								-o $TUPDIR/${nameIDnoext};
		$FSLDIR/bin/fslroi $TUPDIR/${nameIDnoext} ${subID}/func/${epifilename}_tup_checkwarpedtostruc 3 1
		#$FSLDIR/bin/fslroi $TUPDIR/${nameIDnoext} ${epiname}_tup_checkwarpedtostruc 3 1
		#printf "${epiname}_tup_checkwarpedtostruc.nii.gz ${t1struc}\n" >> $QCDIR/qc_filelist_check2struct.txt
		printf "${subID}/func/${epifilename}_tup_checkwarpedtostruc.nii.gz ${t1struc}\n" >> $QCDIR/qc_filelist_check2struct.txt
		
		##########################################################################################
		#
		#	EPI to MNI
		#
		##########################################################################################

		printf "\tCreating EPI to MNI composite warps\n"
		$FSLDIR/bin/convertwarp	-r ${std} \
								--warp1=$TUPDIR/${nameIDnoext}_warp \
								--warp2=${subID}/${STR}.anat/T1_to_MNI_nonlin_field.nii.gz \
								-j $TUPDIR/${nameIDnoext}_jac2MNI.nii.gz \
								-o $TUPDIR/${nameIDnoext}_warp2MNI.nii.gz;

		printf "\tWarping EPI to MNI\n"
		$FSLDIR/bin/applywarp	-i ${epivol} \
								-r ${std} \
								-w $TUPDIR/${nameIDnoext}_warp2MNI.nii.gz \
								--interp=spline \
								-o ${epiname}_tup_checkwarpedtostd;

		printf "\tTOPUP SDC DONE.\n"
	fi
	
	##########################################################################################
	#
	#	NO UNWARPING - Only T1 Registration
	#
	##########################################################################################

	#RAW (distorted) EPI to T1w
	[ ! -d ${subID}/func/sdc_unc ] && { mkdir -m 777 ${subID}/func/sdc_unc; }
	UNCDIR=${subID}/func/sdc_unc
	printf "\n\tSaving SD UNCORRECTED data to: $UNCDIR\n\n"

	$LNIFMRIDIR/lnifmri_util_epireg	--epi=${epivol} \
						--t1=${t1struc} \
						--t1brain=${bstruc} \
						--wmseg=${wmstruc} \
						--out=$UNCDIR/${nameIDnoext} \
						-v --noclean;

	printf "\n\tRegistration of EPI to T1\n"
	$FSLDIR/bin/flirt	-in ${epiimg} \
						-ref ${t1struc} \
						-applyxfm \
						-init $UNCDIR/${nameIDnoext}.mat \
						-interp spline \
						-out $UNCDIR/${nameIDnoext} \
						-refweight ${lesionmaskinv};
						
	$FSLDIR/bin/fslroi $UNCDIR/${nameIDnoext} ${subID}/func/${epifilename}_unc_checkwarpedtostruc 3 1
	#$FSLDIR/bin/fslroi $UNCDIR/${nameIDnoext} ${epiname}_unc_checkwarpedtostruc 3 1
	#printf "${epiname}_unc_checkwarpedtostruc.nii.gz ${t1struc}\n" >> $QCDIR/qc_filelist_check2struct.txt
	printf "${subID}/func/${epifilename}_unc_checkwarpedtostruc.nii.gz ${t1struc}\n" >> $QCDIR/qc_filelist_check2struct.txt

	##########################################################################################
	#
	#	EPI to MNI
	#
	##########################################################################################

	printf "\tCreating EPI to MNI composite warps\n"
	$FSLDIR/bin/convertwarp	-r ${std} \
							--premat=$UNCDIR/${nameIDnoext}.mat \
							--warp2=${subID}/${STR}.anat/T1_to_MNI_nonlin_field.nii.gz \
							-j $UNCDIR/${nameIDnoext}_jac2MNI.nii.gz \
							-o $UNCDIR/${nameIDnoext}_warp2MNI.nii.gz;

	printf "\tWarping EPI to MNI\n"
	$FSLDIR/bin/applywarp	-i ${epivol} \
							-r ${std} \
							-w $UNCDIR/${nameIDnoext}_warp2MNI.nii.gz \
							--interp=spline \
							-o ${epiname}_unc_checkwarpedtostd;

	printf "\tnonlinreg SDC DONE.\n\n"

	ls -d1 `pwd`/sub-*/func/sdc*/${nameID} >> ${LOGDIR}/lnifmri_prep_sdclist.txt

done #for epiimg

#Prepare HTML report for struc and std registration of unwarped images
printf "\tCreating SDC report\n\n"
$FSLDIR/bin/slicesdir -o `cat $QCDIR/qc_filelist_check2struct.txt`
mv slicesdir $QCDIR/lnifmri_prep_func2struc_QC
printf "\n\n"
ls -d1 sub-*/func/*checkwarpedtostd.nii.gz >> $QCDIR/qc_filelist_check2std.txt
$FSLDIR/bin/slicesdir -p ${std} `cat $QCDIR/qc_filelist_check2std.txt`
mv slicesdir $QCDIR/lnifmri_prep_func2std_QC
printf "\n"

#ls -d1 ./sub-*/func/sdc*/*ep2d_bold*[0-9]_+([0-9]).nii.gz > ${LOGDIR}/lnifmri_prep_sdclist.txt

duration=$SECONDS
printf -v exectime "$(($duration / 60))m:$(($duration % 60))s."

NextStep