#!/bin/bash

#
# Script name: cmore_diff.sh
#
# Description: Script to preprocess dMRI images and calculate MD
#  (adapted pipeline from dMRI sequence with 3 diffusion directions)
#
# Authors: Fidel Alfaro-Almagro, Stephen M. Smith & Mark Jenkinson
#
# Copyright 2020 University of Oxford
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#

echo "Start: `date`";
set -x ;

scriptDir="`dirname $0`" ;
subject="$1" ;

# Clean previous runnings of this pipeline
rm   -rf $subject/dMRI/dMRI/ ;
mkdir -p $subject/dMRI/dMRI/dtifit_output ;
cd       $subject/dMRI/dMRI/ ;

# Create symbolic links to the important files and copy bvals and bvecs if absent
ln -s ../raw/AP.nii.gz ./AP.nii.gz ;
ln -s ../raw/PA.nii.gz ./B0_PA.nii.gz ;
for ED in "AP" "PA" ; do
    if [ ! -f ../raw/$ED.bval ] ; then
        cp $scriptDir/$ED.bval ../raw/$ED.bval
        cp $scriptDir/$ED.bvec ../raw/$ED.bvec
    fi
done
ln -s ../raw/AP.bval ./AP.bval ;
ln -s ../raw/AP.bvec ./AP.bvec ;
ln -s ../raw/PA.bval ./PA.bval ;
ln -s ../raw/PA.bvec ./PA.bvec ;

# Generate aqparams file for topup
if [[ -e ../raw/AP.json && -e ../raw/PA.json ]]; then
    echo "Found .json files for dMRI - using to extract total readout time for topup..."
    APECHOSPACING=`grep EffectiveEchoSpacing ../raw/AP.json | awk '{ print $2 }' | sed 's/,//'`
    APROLINES=`fslinfo ../raw/AP.nii.gz | grep dim1 | head -1 | awk '{print $2}'`
    APTOTROTIME=`echo "$APECHOSPACING * ($APROLINES - 1)" | bc | awk '{printf "%f", $0}'`
    echo "AP total readout time: $APTOTROTIME s"

    PAECHOSPACING=`grep EffectiveEchoSpacing ../raw/PA.json | awk '{ print $2 }' | sed 's/,//'`
    PAROLINES=`fslinfo ../raw/PA.nii.gz | grep dim1 | head -1 | awk '{print $2}'`
    PATOTROTIME=`echo "$PAECHOSPACING * ($PAROLINES - 1)" | bc | awk '{printf "%f", $0}'`
    echo "PA total readout time: $PATOTROTIME s"
else
    echo "Cannot find dMRI .json files"
    APTOTROTIME=0.0335
    PATOTROTIME=0.0335
    echo "Assuming AP and PA total readout times are $APTOTROTIME and $PATOTROTIME s (appropriate for GRAPPA = 2)"fi
fi

echo "0 -1 0 $APTOTROTIME" > ./acqparams.txt
echo "0 1 0 $PATOTROTIME" >> ./acqparams.txt
echo 'Diffusion parameters:'
cat ./acqparams.txt


$FSLDIR/bin/eddy_correct ./AP.nii.gz ./Data.nii.gz 0 spline ;
$FSLDIR/bin/fslroi       ./AP.nii.gz ./B0_AP.nii.gz 0 1 ;
$FSLDIR/bin/fslmerge -t  ./B0_4D_AP_PA.nii.gz ./B0_AP.nii.gz ./B0_PA.nii.gz ;

$FSLDIR/bin/topup \
    --imain=./B0_4D_AP_PA.nii.gz \
    --datain=./acqparams.txt \
    --config=$scriptDir/topup_b02b0.cnf \
    --out=TOPUP \
    --fout=TOPUP_fout \
    --jacout=TOPUP_jout \
    --iout=TOPUP_iout ;

# Average output
$FSLDIR/bin/fslmaths TOPUP_iout -Tmean TOPUP_iout_Tmean

# Convert fieldmap to rad/s
$FSLDIR/bin/fslmaths TOPUP_fout -mul 2 -mul 3.1415926536 TOPUP_fout_rad_per_s

# Run bet on the corrected b0 image for passing into oxford_asl
$FSLDIR/bin/bet TOPUP_iout_Tmean TOPUP_iout_Tmean_bet

# Apply topup to the different B0s to get a combined B0
$FSLDIR/bin/applytopup \
    --imain=B0_AP,B0_PA \
    --topup=TOPUP \
    --datain=./acqparams.txt \
    --inindex=1,2 \
    --out=./hifib0.nii.gz ;

# Run a loose (-f 0.2) brain extraction on the B0
$FSLDIR/bin/bet ./hifib0.nii.gz ./hifib0_brain.nii.gz -R -m -f 0.2;

# Fill holes, and slightly erode the brain mask
$FSLDIR/bin/fslmaths \
    ./hifib0_brain_mask.nii.gz \
    -fillh \
    -kernel sphere 2.5 \
    -ero hifib0_brain_mask_ero.nii.gz

# Run applytopup on the data
$FSLDIR/bin/applytopup \
    --imain=./Data.nii.gz \
    --datain=./acqparams.txt \
    --inindex=1 \
    --topup=TOPUP \
    --out=./Data_corr.nii.gz \
    --method=jac \
    --interp=spline ;

# Perform GDC on the data
$BB_BIN_DIR/bb_pipeline_tools/bb_GDC_half_voxel \
    --workingdir=./data_GDC \
    --in=./hifib0.nii.gz \
    --out=./hifib0_ud.nii.gz \
    --owarp=./hifib0_ud_warp.nii.gz ;

# Estimate a registration from B0 to T1 space
$FSLDIR/bin/flirt \
    -in ./hifib0_ud.nii.gz \
    -ref ../../T1/T1.nii.gz \
    -out ./hifib0_ud_to_T1.nii.gz \
    -omat ./hifib0_ud_to_T1.mat \
    -cost mutualinfo \
    -interp spline ;

# Invert the previous transformation
${FSLDIR}/bin/convert_xfm \
    -omat ./T1_to_hifib0_ud.mat \
    -inverse ./hifib0_ud_to_T1.mat ;

# Apply GDC to the data after eddy_correct
${FSLDIR}/bin/applywarp --rel \
    -i ./Data_corr.nii.gz \
    -r ./hifib0_ud.nii.gz \
    -o ./Data_corr_ud.nii.gz \
    -w ./hifib0_ud_warp.nii.gz \
    --interp=spline ;

# Run dtifit on the data
$FSLDIR/bin/dtifit \
    -m ./hifib0_brain_mask_ero.nii.gz \
    -k ./Data_corr_ud.nii.gz \
    -r ./AP.bvec \
    -b ./AP.bval \
    -o ./dtifit_output/dti ;

set +x ;
echo "End: `date`" ;