#!/bin/bash

if [ $# -ne 1 ];
then

echo "Please provide path to working directory. E.g., $0 /Your/data/folder/ds001226-download"

else

dir=$1

for f in sub-CON01 sub-CON02
do

#TOPUP correction for susceptibility-induced distortions

#extract volumes without diffusion weighting
fslroi $dir/$f/ses-preop/dwi/${f}_ses-preop_acq-AP_dwi.nii.gz $dir/$f/ses-preop/dwi/nodif_AP 0 1
fslroi $dir/$f/ses-preop/dwi/${f}_ses-preop_acq-PA_dwi.nii.gz $dir/$f/ses-preop/dwi/nodif_PA 0 1

#merge AP and PA nodif images
fslmerge -t $dir/$f/ses-preop/dwi/AP_PA_b0 $dir/$f/ses-preop/dwi/nodif_AP $dir/$f/ses-preop/dwi/nodif_PA 
 
#create acqparams.txt (get total readout time from json file)
rm -rf $dir/$f/ses-preop/dwi/acqparams.txt

echo "0 -1 0 0.0266003" >> $dir/$f/ses-preop/dwi/acqparams.txt
echo "0 1 0 0.0266003" >> $dir/$f/ses-preop/dwi/acqparams.txt

#run topup

topup --imain=$dir/$f/ses-preop/dwi/AP_PA_b0 \
      --datain=$dir/$f/ses-preop/dwi/acqparams.txt \
      --config=b02b0.cnf \
      --out=$dir/$f/ses-preop/dwi/topup_AP_PA_b0 \
      --iout=$dir/$f/ses-preop/dwi/topup_AP_PA_b0_iout \
      --fout=$dir/$f/ses-preop/dwi/topup_AP_PA_b0_fout


#check the topup results
fsleyes $dir/$f/ses-preop/dwi/topup_AP_PA_b0_iout.nii.gz $dir/$f/ses-preop/dwi/AP_PA_b0.nii.gz

#run EDDY for motion and eddy current correction

#generate a brain mask using the corrected b0 images

#average the two corrected b0 volumes
fslmaths $dir/$f/ses-preop/dwi/topup_AP_PA_b0_iout -Tmean $dir/$f/ses-preop/dwi/hifi_nodif

#apply BET on the averaged b0
bet $dir/$f/ses-preop/dwi/hifi_nodif $dir/$f/ses-preop/dwi/hifi_nodif_brain -m -f 0.2

#check BET has generated an accurate brain mask
fsleyes $dir/$f/ses-preop/dwi/hifi_nodif_brain

#create index.txt file (for a dataset with 102 volumes, such as the one being used for this demo - to find out number of volumes use 'fslinfo data.nii')

rm -rf $dir/$f/ses-preop/dwi/index.txt

indx=""
for ((i=1; i<=102; i+=1)); 
do indx="$indx 1"; 
done
echo $indx > $dir/$f/ses-preop/dwi/index.txt


#run EDDY
eddy_openmp --imain=$dir/$f/ses-preop/dwi/${f}_ses-preop_acq-AP_dwi.nii.gz \
            --mask=$dir/$f/ses-preop/dwi/hifi_nodif_brain_mask \
            --index=$dir/$f/ses-preop/dwi/index.txt \
            --acqp=$dir/$f/ses-preop/dwi/acqparams.txt \
            --bvecs=$dir/$f/ses-preop/dwi/${f}_ses-preop_acq-AP_dwi.bvec \
            --bvals=$dir/$f/ses-preop/dwi/${f}_ses-preop_acq-AP_dwi.bval \
            --fwhm=0 \
            --topup=$dir/$f/ses-preop/dwi/topup_AP_PA_b0 \
            --flm=quadratic \
            --out=$dir/$f/ses-preop/dwi/eddy_unwarped_images \
            --cnr_maps
#           --data_is_shelled


#run eddy quality crontrol 
eddy_quad $dir/$f/ses-preop/dwi/eddy_unwarped_images \
          -idx $dir/$f/ses-preop/dwi/index.txt \
          -par $dir/$f/ses-preop/dwi/acqparams.txt \
          -m $dir/$f/ses-preop/dwi/hifi_nodif_brain_mask \
          -b $dir/$f/ses-preop/dwi/${f}_ses-preop_acq-AP_dwi.bval 

#DTIFIT - fitting the diffusion tensor

#select only bvals 0, 700 and 1200 for tensor fitting
select_dwi_vols $dir/$f/ses-preop/dwi/eddy_unwarped_images.nii.gz \
                $dir/$f/ses-preop/dwi/${f}_ses-preop_acq-AP_dwi.bval \
                $dir/$f/ses-preop/dwi/data_for_dtifit \
                0 -b 700 -b 1200 \
                -obv $dir/$f/ses-preop/dwi/eddy_unwarped_images.eddy_rotated_bvecs

#fit the diffusion tensor model to this subset of data
dtifit --data=$dir/$f/ses-preop/dwi/data_for_dtifit.nii.gz \
       --mask=$dir/$f/ses-preop/dwi/hifi_nodif_brain_mask \
       --bvecs=$dir/$f/ses-preop/dwi/data_for_dtifit.bvec \
       --bvals=$dir/$f/ses-preop/dwi/data_for_dtifit.bval \
       --out=$dir/$f/ses-preop/dwi/dti \
       --wls --sse  

#check dtifit output
fsleyes $dir/$f/ses-preop/dwi/dti_FA.nii.gz $dir/$f/ses-preop/dwi/dti_MD.nii.gz $dir/$f/ses-preop/dwi/dti_V1.nii.gz $dir/$f/ses-preop/dwi/dti_sse.nii.gz


done

fi