run

#!/bin/bash

#show commands running
set -x

#terminate if any command fails
set -e


###############################################################################
# Adapted from a run file from Flywheel
CONTAINER="[garikoitz/rtppreproc]"
echo -e "$CONTAINER  Initiated"


###############################################################################
# Built to flywheel-v0 spec.
FLYWHEEL_BASE=/flywheel/v0
OUTPUT_DIR=$FLYWHEEL_BASE/output
INPUT_DIR=$FLYWHEEL_BASE/input
MANIFEST=$FLYWHEEL_BASE/manifest.json
CONFIG_FILE=$FLYWHEEL_BASE/config.json




###############################################################################
# Configure the ENV

## add N4BiasFieldCorrection to path
# chmod +x /etc/fsl/5.0/fsl.sh
source /etc/fsl/5.0/fsl.sh
export ANTSPATH=/usr/lib/ants
PATH=$PATH:/usr/lib/ants
PATH=$PATH:/mrtrix3/bin
FSLDIR=/usr/share/fsl/5.0
PATH=$PATH:$FSLDIR/bin
LD_LIBRARY_PATH=/usr/lib/fsl/5.0:/usr/share/fsl/5.0/bin
FSLBROWSER=/etc/alternatives/x-www-browser
FSLCLUSTER_MAILOPTS=n
FSLLOCKDIR=
FSLMACHINELIST=
FSLMULTIFILEQUIT=TRUE
FSLOUTPUTTYPE=NIFTI_GZ
FSLREMOTECALL=
FSLTCLSH=/usr/bin/tclsh
FSLWISH=/usr/bin/wish
POSSUMDIR=/usr/share/fsl/5.0
# add base fw so that it can find the awkflipbvec
PATH=$PATH:/$FLYWHEEL_BASE
# Same for python3 in miniconda
export PATH=$PATH:/conda/miniconda3/bin
###############################################################################
# Initialize config parameters

# rpe='  '
acqd='  '
denoise='  '
degibbs='  '
eddy='  '
bias='  '
ricn='  '
norm='  '
nval='  '
anatalign='  '
doreslice='  '
reslice='  '
saveoutput='  '
bias_method='  '
antsb='  '
antsc='  '
antss='  '
eddy_data_is_shelled='  '
eddy_slm='  '
eddy_niter='  '
eddy_repol='  '
eddy_mporder='  '
eddy_slspec='  '

###############################################################################
# Parse config options from CONFIG file or MANIFEST

function parse_config {

  CONFIG_FILE=$FLYWHEEL_BASE/config.json
  MANIFEST_FILE=$FLYWHEEL_BASE/manifest.json

  if [[ -f $CONFIG_FILE ]]; then
    echo -e "$(cat $CONFIG_FILE | jq -r '.config.'$1)"
  else
    CONFIG_FILE=$MANIFEST_FILE
    echo -e "$(cat $MANIFEST_FILE | jq -r '.config.'$1'.default')"
  fi

}

# String parsing
# config_rpe="$(parse_config 'rpe')"
config_acqd="$(parse_config 'acqd')"
config_denoise="$(parse_config 'denoise')"
config_degibbs="$(parse_config 'degibbs')"
config_eddy="$(parse_config 'eddy')"
config_bias="$(parse_config 'bias')"
config_ricn="$(parse_config 'ricn')"
config_norm="$(parse_config 'norm')"
config_nval="$(parse_config 'nval')"
config_anatalign="$(parse_config 'anatalign')"
config_doreslice="$(parse_config 'doreslice')"
config_reslice="$(parse_config 'reslice')"
config_saveoutput="$(parse_config 'saveoutput')"
config_bias_method="$(parse_config 'bias_method')"
config_antsb="$(parse_config 'antsb')"
config_antsc="$(parse_config 'antsc')"
config_antss="$(parse_config 'antss')"
config_eddy_data_is_shelled="$(parse_config 'eddy_data_is_shelled')"
config_eddy_slm="$(parse_config 'eddy_slm')"
config_eddy_niter="$(parse_config 'eddy_niter')"
config_eddy_repol="$(parse_config 'eddy_repol')"
config_eddy_mporder="$(parse_config 'eddy_mporder')"
config_eddy_slspec="$(parse_config 'eddy_slspec')"
config_topup_lambda="$(parse_config 'topup_lambda')"


###############################################################################
# Generate paths to input data


# We do not use their method of having it in the config
DIFF=$(find $INPUT_DIR/DIFF/* -name "*.nii*" | head -1)
BVAL=$(find $INPUT_DIR/BVAL/* -name "*bval*" | head -1)
BVEC=$(find $INPUT_DIR/BVEC/* -name "*bvec*" | head -1)


## optional t1w, only required if alignment with anatomical required 
if [[ -d  $INPUT_DIR/ANAT ]]; then
    ANAT=$(find $INPUT_DIR/ANAT/* -name "*.nii*" | head -1)  ## optional
fi
## optional FS, alignment to FS mask if FS dir exists

if [[ -d $INPUT_DIR/FSMASK ]]; then
    FS=true
	FSMASK=$(find $INPUT_DIR/FSMASK/* -name "*.nii*" | head -1) 
fi
## optional reverse phase encoded (rpe) inputs
if [[ -d  $INPUT_DIR/RDIF ]]; then
    RDIF=$(find $INPUT_DIR/RDIF/* -name "*.nii*" | head -1)  ## optional
    RBVL="NO"
	RBVC="NO"
fi

if [[ -d  $INPUT_DIR/RBVL ]]; then
    RBVL=$(find $INPUT_DIR/RBVL/*  -name "*bval*" | head -1) ## optional
fi
if [[ -d  $INPUT_DIR/RBVC ]]; then
    RBVC=$(find $INPUT_DIR/RBVC/* -name "*bvec*" | head -1) ## optional
fi




###############################################################################
# From now on is their code. Just copy the variables to their variable names
# A we want
cd $OUTPUT_DIR





## acquisition direction: RL, PA, IS
ACQD=$config_acqd

## switches for potentially optional steps
DO_DENOISE=$config_denoise
DO_DEGIBBS=$config_degibbs
DO_EDDY=$config_eddy
DO_BIAS=$config_bias
DO_RICN=$config_ricn
DO_NORM=$config_norm
DO_ANATALIGN=$config_anatalign
NEW_RES=$config_reslice
NORM=$config_nval
DO_RESLICE=$config_doreslice
DO_SAVEOUTPUT=$config_saveoutput

## assign output space of final data if anat align not called
out=proc

## diffusion file that changes name based on steps performed
difm=dwi
mask=b0_dwi_brain_mask



## create / remove old tmp folders / previous run files explicitly
rm -rf $OUTPUT_DIR/tmp $OUTPUT_DIR/eddyqc cor1.b cor2.b corr.b
## create temp folders explicitly
mkdir -p $OUTPUT_DIR/tmp

common="-quiet -force"


## read in eddy options
# RPE=$config_rpe  ## optional

## switch and advanced options for bias correction
BIAS_METHOD=$config_bias_method
ANTSB=$config_antsb
ANTSC=$config_antsc
ANTSS=$config_antss

## fill in arguments common to all dwifslpreproc calls
common_fslpreproc="-eddy_mask ${mask}.mif -eddyqc_all $OUTPUT_DIR/eddyqc -scratch $OUTPUT_DIR/tmp"

## add advanced options to eddy call
eddy_data_is_shelled=$config_eddy_data_is_shelled
eddy_slm=$config_eddy_slm
eddy_niter=$config_eddy_niter
eddy_repol=$config_eddy_repol
eddy_mporder=$config_eddy_mporder

eddy_options=" " ## must contain at least 1 space according to mrtrix doc
[ "$eddy_repol" == "true" ] && eddy_options="$eddy_options --repol"
[ "$eddy_data_is_shelled" == "true" ] && eddy_options="$eddy_options --data_is_shelled"
eddy_options="$eddy_options --slm=$eddy_slm"
eddy_options="$eddy_options --niter=$eddy_niter"
[ "$eddy_mporder" != "0" ] && eddy_options="$eddy_options --mporder=$eddy_mporder"



## add add advanced options for topup call
topup_lambda=$config_topup_lambda
topup_options=" "
[ "$topup_lambda" != "0.005,0.001,0.0001,0.000015,0.000005,0.0000005,0.00000005,0.0000000005,0.00000000001" ] && topup_options="$topup_options --lambda=$topup_lambda"




echo "Converting input files to mrtrix format..."

## convert input diffusion data into mrtrix format
mrconvert -fslgrad $BVEC $BVAL $DIFF raw1.mif --export_grad_mrtrix raw1.b $common

## if the second input exists
if [[ -e $RDIF ]]; then

    ## convert it to mrtrix format
    mrconvert -fslgrad $RBVC $RBVL $RDIF raw2.mif --export_grad_mrtrix raw2.b $common

fi


## determine the type of acquisition for dwipreproc eddy options
if [[ ! $RDIF ]];
then 

    ## if no second sequence, override to the only option
    RPE="none"

else
    
    ## grab the size of each sequence
    nb0F=`mrinfo -size raw1.mif | grep -oE '[^[:space:]]+$'`
    nb0R=`mrinfo -size raw2.mif | grep -oE '[^[:space:]]+$'`

    echo "Forward phase encoded dwi volume has $nb0F volumes."
    echo "Reverse phase encoded dwi volume has $nb0R volumes."
    
    ## check the size of the inputs
    if [ $nb0F -eq $nb0R ];
    then
    ## if they match, it's "all"
    RPE="all"
    ## just because the # of volumes match doesn't mean they're valid
    else
    ## if they don't, it's "pairs"
    RPE="pairs"

    ## if the last dim is even
    if [ $(($nb0R%2)) == 0 ];
    then
        ## pass the file - no assurance it's valid volumes, just a valid number of them
        echo "The RPE file has an even number of volumes. No change was made."
    else
        ## drop any volumes w/ a sufficiently high bval to be a direction - often makes an odd sequence even
        echo "The RPE file has an odd number of volumes. Only the b0 volumes were extracted."
        dwiextract -bzero raw2.mif raw2.mif $common
        ob0=`mrinfo -size raw2.mif | grep -oE '[^[:space:]]+$'`
        echo "This should be an even number: $ob0"
        ## this doesn't stop or exit if it's still odd...
    fi

    fi
    
fi

echo "RPE assigned as: $RPE"









echo "Identifying correct gradient orientation..."

if [[ $RPE == "all" ]]; then

    ## merge them
    mrcat raw1.mif raw2.mif raw.mif $common
    cat raw1.b raw2.b > raw.b

    echo "Creating processing mask..."

    ## create mask from merged data
    dwi2mask raw.mif ${mask}.mif $common

    ## check and correct gradient orientation and create corrected image
    dwigradcheck raw.mif -grad raw.b -mask ${mask}.mif -export_grad_mrtrix corr.b -scratch $OUTPUT_DIR/tmp $common
    mrconvert raw.mif -grad corr.b ${difm}.mif $common
else

    echo "Creating processing mask..."

    ## create mask
    dwi2mask raw1.mif ${mask}.mif $common

    ## check and correct gradient orientation and create corrected image
    dwigradcheck raw1.mif -grad raw1.b -mask ${mask}.mif -export_grad_mrtrix cor1.b -scratch $OUTPUT_DIR/tmp $common
    mrconvert raw1.mif -grad cor1.b ${difm}.mif $common

    if [[ -e raw2.mif ]]; then
	dwi2mask raw2.mif rpe_${mask}.mif $common
	cp raw2.b cor2.b
	mrconvert raw2.mif -grad cor2.b rpe_${difm}.mif $common
    fi

fi

## perform PCA denoising
if [ $DO_DENOISE == "true" ] || [ $DO_RICN == "true" ]; then

    echo "Performing PCA denoising..."
    dwidenoise -extent 5,5,5 -noise fpe_noise.mif -estimator Exp2 ${difm}.mif ${difm}_denoise.mif $common

    ## if the second volume exists, denoise as well and average the noise volumes together    
    if [ -e rpe_${difm}.mif ]; then
        dwidenoise -extent 5,5,5 -noise rpe_noise.mif -estimator Exp2 rpe_${difm}.mif rpe_${difm}_denoise.mif $common
        mrcalc fpe_noise.mif rpe_noise.mif -add 2 -divide noise.mif $common
    else
        cp fpe_noise.mif noise.mif
    fi

    ## if denoise is true, use the denoised volume on the next steps
    ## otherwise it just makes the noise for the rician
    if [ $DO_DENOISE == "true" ]; then
        difm=${difm}_denoise
    fi

fi



## if scanner artifact is found
if [[ $DO_DEGIBBS == "true" ]]; then

    echo "Performing Gibbs ringing correction..."
    mrdegibbs -nshifts 20 -minW 1 -maxW 3 ${difm}.mif ${difm}_degibbs.mif $common

    if [[ -e rpe_${difm}.mif ]]; then
	mrdegibbs -nshifts 20 -minW 1 -maxW 3 rpe_${difm}.mif rpe_${difm}_degibbs.mif $common
    fi

    difm=${difm}_degibbs

fi
 
## perform eddy correction with FSL
if [[ $DO_EDDY == "true" ]]; then

    if [[ $RPE == "none" ]]; then

		echo "Performing FSL eddy correction..."
		dwifslpreproc ${difm}.mif ${difm}_eddy.mif -rpe_none -pe_dir ${ACQD} -eddy_options "$eddy_options" $common_fslpreproc $common
		difm=${difm}_eddy
    fi

    if [[ $RPE == "pairs" ]]; then

		echo "Performing FSL topup and eddy correction ..."
    	## pull and merge the b0s
    	dwiextract -bzero ${difm}.mif fpe_b0.mif $common
    	dwiextract -bzero rpe_${difm}.mif rpe_b0.mif $common ## maybe redundant?


		nb0F=`mrinfo -size fpe_b0.mif | grep -oE '[^[:space:]]+$'`
		nb0R=`mrinfo -size rpe_b0.mif | grep -oE '[^[:space:]]+$'`
		if [ $nb0F -gt $nb0R ];
		then
		mrconvert fpe_b0.mif tmp.mif -coord 3 0:$(($nb0R-1)) $common
		mrcat tmp.mif rpe_b0.mif b0_pairs.mif -axis 3 $common 
		elif [ $nb0F -eq $nb0R ];
		then
		mrcat fpe_b0.mif rpe_b0.mif b0_pairs.mif -axis 3 $common
		else
		mrconvert rpe_b0.mif tmp.mif -coord 3 0:$(($nb0F-1)) $common
		mrcat fpe_b0.mif tmp.mif b0_pairs.mif -axis 3 $common
		fi




    	## call to dwifslpreproc w/ new options
    	dwifslpreproc ${difm}.mif ${difm}_eddy.mif -rpe_pair -se_epi b0_pairs.mif -pe_dir ${ACQD} -align_seepi -topup_options "$topup_options" -eddy_options "$eddy_options" $common_fslpreproc $common
		difm=${difm}_eddy

    fi

    if [[ $RPE == "all" ]]; then

		echo "Performing FSL eddy correction for merged input DWI sequences..."
    	dwifslpreproc ${difm}.mif ${difm}_eddy.mif -rpe_all -pe_dir ${ACQD} -topup_options "$topup_options" -eddy_options "$eddy_options" $common_fslpreproc $common
		difm=${difm}_eddy

    fi

fi

## rebuild mask after eddy motion - necessary?
dwi2mask ${difm}.mif ${mask}.mif $common




## compute bias correction with ANTs on dwi data
if [[ $DO_BIAS == "true" ]]; then

    if [ $BIAS_METHOD == "ants" ]; then
        echo "Performing bias correction with ANTs..."
        dwibiascorrect ants -ants.b $ANTSB -ants.c $ANTSC -ants.s $ANTSS -mask ${mask}.mif ${difm}.mif ${difm}_bias.mif -scratch $OUTPUT_DIR/tmp $common
    fi

    if [ $BIAS_METHOD == "fsl" ]; then
        echo "Performing bias correction with FSL..."
        dwibiascorrect fsl -mask ${mask}.mif ${difm}.mif ${difm}_bias.mif -scratch $OUTPUT_DIR/tmp $common   
    fi

	difm=${difm}_bias


fi

## perform Rician background noise removal
if [[ $DO_RICN == "true" ]]; then

    echo "Performing Rician background noise removal..."
    mrinfo ${difm}.mif -export_grad_mrtrix tmp.b  $common
    mrcalc noise.mif -finite noise.mif 0 -if lowbnoisemap.mif  $common
    mrcalc ${difm}.mif 2 -pow lowbnoisemap.mif 2 -pow -sub -abs -sqrt - $common | mrcalc - -finite - 0 -if tmp.mif $common
    difm=${difm}_ricn
    mrconvert tmp.mif -grad tmp.b ${difm}.mif $common
    rm -f tmp.mif tmp.b

fi

## perform intensity normalization of dwi data
if [[ $DO_NORM == "true" ]]; then

    echo "Performing intensity normalization..."

    ## create fa wm mask of input subject
    dwi2tensor -mask ${mask}.mif  ${difm}.mif - $common | tensor2metric  - -fa - $common | mrthreshold  -abs 0.5 - wm.mif $common


    ## normalize intensity of generous FA white matter mask to 1000
    dwinormalise individual -intensity $NORM ${difm}.mif wm.mif ${difm}_norm.mif $commont
    difm=${difm}_norm

fi






echo "Creating dwi space b0 reference images..."

## create b0 and mask image in dwi space on forward direction only
dwiextract ${difm}.mif - -bzero $common | mrmath - mean b0_dwi.mif -axis 3 $common

## compute dwi mask for processing
dwi2mask ${difm}.mif ${mask}.mif $common

## convert to nifti for alignment to anatomy later on
mrconvert b0_dwi.mif b0_dwi.nii.gz $common
mrconvert ${mask}.mif ${mask}.nii.gz $common

## apply mask to image
fslmaths b0_dwi.nii.gz -mas ${mask}.nii.gz b0_dwi_brain.nii.gz

## align diffusion data to T1 acpc anatomy
if [[ $DO_ANATALIGN == "true" ]]; then
        
	## create local copy of anat
    ## instead of a single copy, make sure the T1 is RAS, as we do the same with the DWI
    # cp $ANAT ./t1.nii.gz
    mrconvert $common -stride 1,2,3 $ANAT ./t1.nii.gz
	ANAT=t1
    
	echo "Running brain extraction on anatomy..."
	if [[ $FS == "true" ]]; then
		## use FS mask if FS dir exist 
		mrconvert $common  -stride 1,2,3 $FSMASK ./t1_brain.nii.gz
	else
       	## create t1 mask
       	bet ${ANAT}.nii.gz ${ANAT}_brain -R -B -m
	fi

    echo "Aligning dwi data with  anatomy..."

    # epi_reg from FSL was doing a poor job, so here I maintain the old code and 
	# add the new version with ANTs
	## compute BBR registration corrected diffusion data to  anatomy
    # epi_reg --epi=b0_dwi_brain.nii.gz --t1=${ANAT}.nii.gz --t1brain=${ANAT}_brain.nii.gz --out=dwi2anatalign
    ## apply the transform w/in mrtrix, correcting gradients
    # transformconvert dwi2anatalign.mat b0_dwi_brain.nii.gz ${ANAT}_brain.nii.gz flirt_import dwi2anatalign_mrtrix.mat $common
    # mrtransform -linear dwi2anatalign_mrtrix.mat ${difm}.mif ${difm}_anatalign.mif $common

	# This is the new code with ANTs
	export ANTSPATH=/usr/bin/
	## compute transform mat using ANTs
	antsRegistrationSyN.sh -d 3 -f ${ANAT}.nii.gz -m b0_dwi_brain.nii.gz -t r -o ants
	## Apply the transform matrix
	ConvertTransformFile 3 ants0GenericAffine.mat ants0GenericAffine.txt
	transformconvert ants0GenericAffine.txt itk_import dwi2anatalign_mrtrix.txt
	mrtransform -linear dwi2anatalign_mrtrix.txt ${difm}.mif ${difm}_anatalign.mif $common

	difm=${difm}_anatalign

    ## assign output space label
    out=anatalign

fi

if [[ $DO_RESLICE == "true" ]]; then

    echo "Reslicing diffusion data to requested isotropic voxel size..."

    ## sed to turn possible decimal into p
    VAL=`echo $NEW_RES | sed s/\\\./p/g`
	echo "Reslicing diffusion data to the requested isotropic voxel size of $VAL mm^3..."

    mrgrid ${difm}.mif regrid -voxel $NEW_RES ${difm}_${VAL}mm.mif $common
    difm=${difm}_${VAL}mm

else

    ## append voxel size in mm to the end of file, rename
    VAL=`mrinfo -spacing dwi.mif | awk {'print $1'} | sed s/\\\./p/g`
    echo VAL: $VAL
    mv ${difm}.mif ${difm}_${VAL}mm.mif
    difm=${difm}_${VAL}mm

fi

echo "Creating $out space b0 reference images..."

## create final b0 / mask
dwiextract ${difm}.mif - -bzero $common | mrmath - mean b0_${out}.mif -axis 3 $common
dwi2mask ${difm}.mif b0_${out}_brain_mask.mif $common

## create output space b0s
mrconvert b0_${out}.mif -stride 1,2,3,4 b0_${out}.nii.gz $common
mrconvert b0_${out}_brain_mask.mif -stride 1,2,3,4 b0_${out}_brain_mask.nii.gz $common
fslmaths b0_${out}.nii.gz -mas b0_${out}_brain_mask.nii.gz b0_${out}_brain.nii.gz

echo "Creating preprocessed dwi files in $out space..."

## convert to nifti / fsl output for storage
mrconvert ${difm}.mif -stride 1,2,3,4 dwi.nii.gz -export_grad_fsl dwi.bvecs dwi.bvals -export_grad_mrtrix ${difm}.b -json_export ${difm}.json $common



# Notes agter 1.0.9: this is no longer necessary, it works well with the correct bvec
# [Here we can change the strides to -1,2,3,4 or flip the bvec. Later one we would prefer to have RAS instead of LAS, then here we are going to change the x only. 
# This requires more work on the future and more thorough testing for when the data is coming incorrectly, as with the newest FSL versions we need to have it right from the beginning.
# awkflipbvec.sh dwi.bvecs dwi_noXflip.bvecs
# python3 $FLYWHEEL_BASE/flipX.py dwi.bvecs dwi_noXflip.bvecs]


## export a lightly structured text file (json?) of shell count / lmax
echo "Writing text file of basic sequence information..."

## parse the number of shells / determine if a b0 is found
if [ ! -f b0_dwi.mif ]; then
    echo "No b-zero volumes present"
    nshell=`mrinfo -shell_bvalues -bvalue_scaling false ${difm}.mif | wc -w`
    shell=$nshell
    b0s=0
    lmaxs=`dirstat ${difm}.b | grep lmax | awk '{print $8}' | sed "s|:||g"`
else
    nshell=`mrinfo -shell_bvalues -bvalue_scaling false ${difm}.mif | wc -w`
    shell=$(($nshell-1)) ## at least 1 b0 found
    b0s=`mrinfo -shell_sizes ${difm}.mif | awk '{print $1}'`
    lmaxs='0 '`dirstat ${difm}.b | grep lmax | awk '{print $8}' | sed "s|:||g"`
fi


## add file name to summary.txt
echo ${difm} > summary.txt

if [[ $shell -gt 1 ]]; then
    echo multi-shell: $shell total shells >> summary.txt
else
    echo single-shell: $shell total shell >> summary.txt
fi

## print the number of b0s
echo Number of b0s: $b0s >> summary.txt

echo >> summary.txt
echo shell / count / lmax >> summary.txt

## echo basic shell count summaries
mrinfo -shell_bvalues ${difm}.mif >> summary.txt
mrinfo -shell_sizes ${difm}.mif >> summary.txt

## echo max lmax per shell
lmaxs=`dirstat ${difm}.b | grep lmax | awk '{print $8}' | sed "s|:||g"`
echo $lmaxs >> summary.txt

## print into log
cat summary.txt

echo "Cleaning up working directory..."

## cleanup
if [[ $DO_SAVEOUTPUT == "false" ]]; then
    find . -maxdepth 1 -mindepth 1 -type f -name "*.mif" ! -name "${difm}.mif" -delete
    find . -maxdepth 1 -mindepth 1 -type f -name "*.b" ! -name "${difm}.b" -delete
    rm -f *.mif
    rm -f *.b
    rm -f *fast*.nii.gz
    rm -f *init.mat
    rm -f dwi2anatalign.nii.gz
    rm -rf ./tmp
fi
###############################################################################
# Check status and exit accordingly
# Ask Michael how he wants to do this, convert all the previous calls to functions?


## check for the output for success
if [[ -f ${OUTPUT_DIR}/dwi.nii.gz ]]; then
    echo "Changing permissions"
    chmod -R 755 ${OUTPUT_DIR} 
    echo "$CONTAINER  Success!"
    exit 0
else
    echo "$CONTAINER  Failure detected. Final dwi data files are missing!"
    exit 1
fi