#!/bin/bash

module load fsl
module load mrtrix3
module load matlab/r2014b

timestamp() {
  date +"%T"
}

SCRIPTDIR=/home/soostermeijer/Monash076/Sanne/Copy_ADS_Simon
PROJDIR=/home/soostermeijer/Monash076/Sanne/Copy_ADS_Simon
DATADIR=$PROJDIR/T3/dicom_denamed_mri2

SUBJIDS=$(<$DATADIR/c.txt) # SubjectIDs scratch

# Reference image for EDDY correction
refImg=0 

#MODULES
	MODULE0=1 #mrconvert
	MODULE1=1 #EDDY
	MODULE2=1 #GEN DTI&FA
	MODULE3=1 #Registration

for SUBJ in $SUBJIDS; do
	timestamp
	# Dynamic directories
	DICOMDIR=$DATADIR/$SUBJ
	DIFFDIR=$DATADIR/$SUBJ/DiffusionBETf.01
	
	# clean outputs and reinitialise
	if [ ! -d $DIFFDIR ]; then mkdir $DIFFDIR; else echo "Cleaning outputs"; rm -rf $DIFFDIR; mkdir $DIFFDIR; fi

	################################ MODULE 0: MRtrix mrconvert #######################################
	if [ $MODULE0 = "1" ]; then
		echo "Running MODULE 0: MRtrix mrconvert: $SUBJ"
		# scan dicom folder and write selection options to selectseries.txt
		echo q | mrinfo $DICOMDIR &> $DIFFDIR/selectseries.txt
		# select specified series and write diffusion-weighted images to dwi.mif
		SERIES=DTI_MGH_60Dir
		echo `grep ${SERIES} $DIFFDIR/selectseries.txt | cut -c 3-4` | mrinfo -export_grad_mrtrix $DIFFDIR/gradients.b $DICOMDIR
		# echo `grep ${SERIES} $DIFFDIR/selectseries.txt | cut -c 3-4` | mrinfo -export_grad_fsl $DIFFDIR/bvecs $DIFFDIR/bvals $DICOMDIR
		echo `grep ${SERIES} $DIFFDIR/selectseries.txt | cut -c 3-4` | mrconvert $DICOMDIR $DIFFDIR/dwi.nii
		echo "Finished MODULE 0: MRtrix mrconvert: $SUBJ"
	else
		echo "Skipping MODULE 0: MRtrix mrconvert: $SUBJ"
	fi
	###################################################################################################

	####################################### MODULE 1: EDDY ############################################
	if [ $MODULE1 = "1" ]; then
		echo "Running MODULE 1: FSL EDDY: $SUBJ"
		cd $DIFFDIR
		#motion correction
		eddy_correct $DIFFDIR/dwi.nii $DIFFDIR/dwi_ec $refImg
		# rotate bvecs
		matlab -nodesktop -nosplash -r "SCRIPTDIR='${SCRIPTDIR}'; addpath(genpath(SCRIPTDIR)); SUBJ='${SUBJ}'; rot_gradients; exit"
		# gen mask (mrtrix method)
		# dwi2mask dwi_ec.nii.gz nodif_brain.nii.gz -grad gradients_rot.b

		# gen mask (FSL method)
		fslroi dwi_ec.nii.gz nodif.nii.gz $refImg 1		
		bet nodif.nii.gz nodif_brain.nii.gz -m -f 0.1

		# erode mask
		# fslmaths nodif_brain_mask.nii.gz -ero $DIFFDIR/eroded_mask.nii.gz
		echo "Finished MODULE 1: FSL EDDY: $SUBJ"
	else
		echo "Skipping MODULE 1: FSL EDDY: $SUBJ"
	fi
	###################################################################################################

	####################################### MODULE 2: DTI ############################################
	if [ $MODULE2 = "1" ]; then
		echo "Running MODULE 2: MRtrix DTI & FA: $SUBJ"
		# convert to .mif format and combine in gradients
		mrconvert -grad $DIFFDIR/gradients_rot.b $DIFFDIR/dwi_ec.nii.gz $DIFFDIR/dwi_ec.mif 
		# gen diffusion tensor
		dwi2tensor $DIFFDIR/dwi_ec.mif $DIFFDIR/dti.mif
		# gen fa and ev maps
		tensor2metric -fa $DIFFDIR/fa.mif -vector $DIFFDIR/ev.mif -mask $DIFFDIR/nodif_brain.nii.gz $DIFFDIR/dti.mif
		# convert fa to .nii in order to run registration step below
		mrconvert $DIFFDIR/fa.mif $DIFFDIR/fa.nii
		echo "Finished MODULE 2: MRtrix DTI & FA: $SUBJ"
	else
		echo "Skipping MODULE 2: MRtrix DTI & FA: $SUBJ"
	fi
	###################################################################################################

	#################################### MODULE 3: Registration #########################################
	if [ $MODULE3 = "1" ]; then
		echo "Running MODULE 3: Registration: $SUBJ"
		# remove outputs
		# rm $DIFFDIR/fa_flirt.nii.gz
		# rm $DIFFDIR/flirt.mat
		# rm $DIFFDIR/fnirt_warp.nii.gz
		# rm $DIFFDIR/fa_fnirt.nii.gz
		# rm $DIFFDIR/fa_reor.nii
		
		cd $DIFFDIR
		fslreorient2std fa.nii fa_reor.nii
		gunzip fa_reor.nii.gz
		
		flirt -ref $FSLDIR/data/standard/FMRIB58_FA_1mm -in $DIFFDIR/fa_reor.nii -cost normmi -searchcost normmi -dof 12 -interp trilinear -out $DIFFDIR/fa_flirt -omat $DIFFDIR/flirt.mat
		fnirt --in=$DIFFDIR/fa_reor.nii --aff=$DIFFDIR/flirt.mat --cout=$DIFFDIR/fnirt_warp --config=FA_2_FMRIB58_1mm
		applywarp --ref=$FSLDIR/data/standard/FMRIB58_FA_1mm --in=$DIFFDIR/fa_reor.nii --warp=$DIFFDIR/fnirt_warp --out=$DIFFDIR/fa_fnirt
		
		echo "Finished MODULE 3: Registration: $SUBJ"
	else
		echo "Skipping MODULE 3: Registration: $SUBJ"
	fi
	timestamp
	#####################################################################################################
done