#!/bin/bash


#TOPUP correction for susceptibility-induced distortions

#extract volumes without diffusion weighting
fslroi sub-CON01/ses-preop/dwi/sub-CON01_ses-preop_acq-AP_dwi.nii.gz sub-CON01/ses-preop/dwi/nodif_AP 0 1
fslroi sub-CON01/ses-preop/dwi/sub-CON01_ses-preop_acq-PA_dwi.nii.gz sub-CON01/ses-preop/dwi/nodif_PA 0 1

#merge AP and PA nodif images
fslmerge -t sub-CON01/ses-preop/dwi/AP_PA_b0 sub-CON01/ses-preop/dwi/nodif_AP sub-CON01/ses-preop/dwi/nodif_PA 
 
#have a look at the merged nodif_PA file
fsleyes sub-CON01/ses-preop/dwi/AP_PA_b0 

#create acqparams.txt (get total readout time from json file)
rm -rf sub-CON01/ses-preop/dwi/acqparams.txt

echo "0 -1 0 0.0266003" >> sub-CON01/ses-preop/dwi/acqparams.txt
echo "0 1 0 0.0266003" >> sub-CON01/ses-preop/dwi/acqparams.txt

#run topup

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


#check the topup results
fsleyes sub-CON01/ses-preop/dwi/topup_AP_PA_b0_fout

fsleyes sub-CON01/ses-preop/dwi/topup_AP_PA_b0_iout.nii.gz sub-CON01/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 sub-CON01/ses-preop/dwi/topup_AP_PA_b0_iout -Tmean sub-CON01/ses-preop/dwi/hifi_nodif

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

#check BET has generated an accurate brain mask
fsleyes sub-CON01/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 sub-CON01/ses-preop/dwi/index.txt

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


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

#check eddy_output
fsleyes sub-CON01/ses-preop/dwi/eddy_unwarped_images

#run eddy quality crontrol 
eddy_quad sub-CON01/ses-preop/dwi/eddy_unwarped_images \
          -idx sub-CON01/ses-preop/dwi/index.txt \
          -par sub-CON01/ses-preop/dwi/acqparams.txt \
          -m sub-CON01/ses-preop/dwi/hifi_nodif_brain_mask \
          -b sub-CON01/ses-preop/dwi/sub-CON01_ses-preop_acq-AP_dwi.bval 

#DTIFIT - fitting the diffusion tensor

#select only bvals 0, 700 and 1200 for tensor fitting
select_dwi_vols sub-CON01/ses-preop/dwi/eddy_unwarped_images.nii.gz \
                sub-CON01/ses-preop/dwi/sub-CON01_ses-preop_acq-AP_dwi.bval \
                sub-CON01/ses-preop/dwi/data_for_dtifit \
                0 -b 700 -b 1200 \
                -obv sub-CON01/ses-preop/dwi/eddy_unwarped_images.eddy_rotated_bvecs

#check bvals file
gedit sub-CON01/ses-preop/dwi/data_for_dtifit.bval

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

#check dtifit output
fsleyes sub-CON01/ses-preop/dwi/dti_FA.nii.gz sub-CON01/ses-preop/dwi/dti_MD.nii.gz sub-CON01/ses-preop/dwi/dti_V1.nii.gz sub-CON01/ses-preop/dwi/dti_sse.nii.gz