#!/bin/tcsh
set nonomatch

# output usage if needed
if ($#argv < 4) then
    echo "usage: topupSingleFolder.script <firstDirDWI.nrrd> <OppositDirDWI.nrrd> <outputBase> <numCores>"
    exit
endif

setenv OMP_NUM_THREADS $argv[4]

echo using $OMP_NUM_THREADS cores

# Data Directories
set BaseDir       = /proj/NIRAL/studies/IBIS/VSA
set DATADIR       = $BaseDir/Data_dMRI
set ProcessDIR    = $BaseDir/processing/dMRI
set ConfigDIR     = $ProcessDIR/config

# Tools
set DTIprepCMD    = /proj/NIRAL/tools/DTIPrep1.2.11b
set unuCMD        = /proj/NIRAL/tools/unu
set DWIconvertCMD = /proj/NIRAL/tools/DWIConvert_4.6

#FSL
set FSLHOME       = /proj/NIRAL/tools/fsl_6.0.3

setenv FSLDIR $FSLHOME
setenv PATH $FSLHOME/bin:${PATH}
source $FSLHOME/etc/fslconf/fsl.csh
echo FSL set to $FSLDIR

module remove python
setenv PATH /proj/NIRAL/tools/Anaconda/anaconda3/bin:$PATH
#use NIRAL anaconda for eddyquant

#topup params
set DTIPrepProto  = $ConfigDIR/protocol_DTIPrep1.2.11_IBIS_VSA.xml
set topupConfig   = $ConfigDIR/fsl_regb02b0.cnf

# Report Files
set todayDate = `date +%y_%m_%d`
set reportFile = $ProcessDIR/reports/Topup_report_${todayDate}.txt
#if (-e $reportFile) rm $reportFile
touch $reportFile

# processing pipeline version
set procVersion = v1.0

#minimum bval to count as b=0 shot for topup computation
set minB0 = 20
#pattern for topup parameter file
set topupParamPattern = "0 PE 0 0.0924"

#setup input data
# current script assumes we have input data in NRRD format
set P1file = $argv[1]
set P2file = $argv[2]
set P1P2base = $argv[3]
set dMRIfolder = $P1file:h

#setup output data
set DTIPrepFolder = $dMRIfolder/DTIPrep1.2.11

set P1dtiprepQC   = $DTIPrepFolder/$P1file:t:r_QCed.nrrd
set P2dtiprepQC   = $DTIPrepFolder/$P2file:t:r_QCed.nrrd
set P1NIFTI       = $dMRIfolder/$P1file:t:r_QCed.nii.gz
set P2NIFTI       = $dMRIfolder/$P2file:t:r_QCed.nii.gz
set B0indexP1File = $dMRIfolder/$P1file:t:r_QCed.b0Index
set B0indexP2File = $dMRIfolder/$P2file:t:r_QCed.b0Index
set B0P1File      = $dMRIfolder/$P1file:t:r_QCed_B0.nii.gz
set B0P2File      = $dMRIfolder/$P2file:t:r_QCed_B0.nii.gz

set P1P2NIFTI     = ${P1P2base}.nii.gz
set B0indexP1P2File = ${P1P2base}.b0Index
set B0P1P2File    = ${P1P2base}_B0.nii.gz

set B0P1P2Corr    = ${P1P2base}_B0_Topup.nii.gz
set P1P2Topup     = ${P1P2base}_topup_fieldcoef.nii.gz
set P1P2Field     = ${P1P2base}_field.nii.gz
set topupParam    = ${P1P2base}_topupParams.txt
set B0P1P2Avg     = ${P1P2base}_AvgB0_Topup.nii.gz
set B0P1P2Mask    = ${P1P2base}_brainmask.nii.gz
set P1P2Corr      = ${P1P2base}_topupEddy_corr.nii.gz
set P1P2CorrNoNeg = ${P1P2base}_topupEddy_corr_noNeg.nii.gz
set P1P2CorrNRRD  = ${P1P2base}_topupEddy_corr_noNeg.nrrd
set indexFile     = ${P1P2base}_indexAPPA.txt
set quadReport    = $P1P2Corr:r:r.qc/qc.pdf

# This is the pipeline
# 1. Run DTIPrep on each file & convert to NIFTI
# 2. Create B0 files and B0 index files
# 3. Combine the B0 and DWI files
# 4. Run topup
# 5. Extract mask from average topup corr B0
# 6. Run Eddy
# 7. Make Corr DWI positive only
# 8. Run QUAD
# 9. Convert to NRRD

#################################################
#################################################
# 1. Run DTIPrep on each file & convert to NIFTI
#################################################
#################################################

echo Running DTIPrep and conversion
echo And ensure that number of slices is even

if (! -e $DTIPrepFolder) mkdir $DTIPrepFolder

if (! -e $P1dtiprepQC) then
#    SKIPPING AS FOR THIS CASE P1 is already run through DTIPrep
    cp $P1file $P1dtiprepQC
#    echo running DTIPrep P1 $P1dtiprepQC >> $reportFile
#    echo $DTIprepCMD -c -w $P1file --xmlProtocol $DTIPrepProto -f $DTIPrepFolder
endif
if (! -e $P2dtiprepQC) then
    echo running DTIPrep P1 $P2dtiprepQC >> $reportFile
    $DTIprepCMD -c -w $P2file --xmlProtocol $DTIPrepProto -f $DTIPrepFolder
endif

if (! -e $P1NIFTI && -e $P1dtiprepQC) then
    echo converting $P1dtiprepQC  >> $reportFile
    echo padding  1 slice >> $reportFile
    unu pad -v 0 -i $P1dtiprepQC -o $P1dtiprepQC:r_pad.nrrd -min 0 0 0 0 -max M M M M+1 -b pad 
    $DWIconvertCMD --conversionMode NrrdToFSL --outputBVectors $P1NIFTI:r:r.bvecs --outputBValues $P1NIFTI:r:r.bvals --outputVolume $P1NIFTI --inputVolume $P1dtiprepQC:r_pad.nrrd 
endif
if (! -e $P2NIFTI && -e $P2dtiprepQC) then
    echo converting $P2dtiprepQC  >> $reportFile
    echo padding  1 slice >> $reportFile
    unu pad -v 0 -i $P2dtiprepQC -o $P2dtiprepQC:r_pad.nrrd -min 0 0 0 0 -max M M M M+1 -b pad 
    $DWIconvertCMD --conversionMode NrrdToFSL --outputBVectors $P2NIFTI:r:r.bvecs --outputBValues $P2NIFTI:r:r.bvals --outputVolume $P2NIFTI --inputVolume $P2dtiprepQC:r_pad.nrrd
endif



#################################################
#################################################
# 2. Create B0 files and B0 index files
#################################################
#################################################

echo creating B0 files

if (! -e $B0indexP1File && -e $P1NIFTI:r:r.bvals) then
    echo creating b0 index file $B0indexP1File from $P1NIFTI:r:r.bvals >> $reportFile
    set allbval = `cat $P1NIFTI:r:r.bvals`
    set numbval = $#allbval
    set curIndex = 1
    set indexSet = ()
    while ($curIndex <= $numbval)
	@ dwiIndex = $curIndex - 1
	if ($allbval[$curIndex] < $minB0) then
	    set indexSet = ($indexSet $dwiIndex)
	endif
	    @ curIndex = $curIndex + 1
    end
    echo $indexSet >! $B0indexP1File
endif
if (! -e $B0indexP2File && -e $P2NIFTI:r:r.bvals) then
    echo creating b0 index file $B0indexP2File from $P2NIFTI:r:r.bvals >> $reportFile
    set allbval = `cat $P2NIFTI:r:r.bvals`
    set numbval = $#allbval
    set curIndex = 1
    set indexSet = ()
    while ($curIndex <= $numbval)
	@ dwiIndex = $curIndex - 1
	if ($allbval[$curIndex] < $minB0) then
	    set indexSet = ($indexSet $dwiIndex)
	endif
	    @ curIndex = $curIndex + 1
    end
    echo $indexSet >! $B0indexP2File
endif

if (! -e $B0P1File && -e $B0indexP1File) then
    echo creating B0 $B0P1File >> $reportFile
    set B0IndexSet = `cat $B0indexP1File`	
    set numbval = $#B0IndexSet	

    set tmpDir = $B0P1File:h/tmpDir
    if (-e $tmpDir) rm -r $tmpDir
    mkdir $tmpDir
    set curIndex = 1
    set fileSet = ()
    while ($curIndex <= $numbval)
	set indexString = "$curIndex"
	if ($curIndex < 10) then
	    set indexString = "00$curIndex"
	endif
	if ($curIndex < 100) then
	    set indexString = "0$curIndex"
	endif
	set fileSet = ($fileSet $tmpDir/B0_$indexString.nii.gz)
	set index = $B0IndexSet[$curIndex]
	fslroi $P1NIFTI $tmpDir/B0_$indexString.nii.gz $index 1
	@ curIndex = $curIndex + 1
    end
    fslmerge  -t $B0P1File $fileSet
    rm -r $tmpDir
endif
if (! -e $B0P2File && -e $B0indexP2File) then
    echo creating B0 $B0P2File >> $reportFile
    set B0IndexSet = `cat $B0indexP2File`	
    set numbval = $#B0IndexSet	

    set tmpDir = $B0P1File:h/tmpDir
    if (-e $tmpDir) rm -r $tmpDir
    mkdir $tmpDir
    set curIndex = 1
    set fileSet = ()
    while ($curIndex <= $numbval)
	set indexString = "$curIndex"
	if ($curIndex < 10) then
	    set indexString = "00$curIndex"
	endif
	if ($curIndex < 100) then
	    set indexString = "0$curIndex"
	endif
	set fileSet = ($fileSet $tmpDir/B0_$indexString.nii.gz)
	set index = $B0IndexSet[$curIndex]
	fslroi $P2NIFTI $tmpDir/B0_$indexString.nii.gz $index 1
	@ curIndex = $curIndex + 1
    end
    fslmerge  -t $B0P2File $fileSet
    rm -r $tmpDir
endif

#################################################
#################################################
# 3. Combine the B0 and DWI files
#################################################
#################################################
echo merge B0s AP and PA

if (! -e $B0P1P2File && -e $B0P1File && -e $B0P2File) then
    echo merging both PE B0  into $B0P1P2File:t >> $reportFile
    fslmerge -t $B0P1P2File $B0P1File  $B0P2File
endif

if (! -e $P1P2NIFTI && -e $P1NIFTI && -e $P2NIFTI) then
    echo merging both PE DWI into $P1P2NIFTI:t >> $reportFile
    fslmerge -t $P1P2NIFTI $P1NIFTI $P2NIFTI
    cat  $P1NIFTI:r:r.bvals $P2NIFTI:r:r.bvals >! $P1P2NIFTI:r:r.bvals
    cat  $P1NIFTI:r:r.bvecs $P2NIFTI:r:r.bvecs >! $P1P2NIFTI:r:r.bvecs
endif

if (! -e $B0indexP1P2File && -e $P1P2NIFTI:r:r.bvals) then
    echo creating b0 index file $B0indexP1P2File from $P1P2NIFTI:r:r.bvals >> $reportFile
    set allbval = `cat $P1P2NIFTI:r:r.bvals`
    set numbval = $#allbval
    set curIndex = 1
    set indexSet = ()
    while ($curIndex <= $numbval)
	@ dwiIndex = $curIndex - 1
	if ($allbval[$curIndex] < $minB0) then
	    set indexSet = ($indexSet $dwiIndex)
	endif
	    @ curIndex = $curIndex + 1
    end
    echo $indexSet >! $B0indexP1P2File
endif

#################################################
#################################################
# 4. Run topup
#################################################
#################################################

echo Run topup

if (! -e $topupParam) then
    echo generating topupParam $topupParam >> $reportFile
    set P1num = `wc -w $B0indexP1File`
    set P2num = `wc -w $B0indexP2File`
    set counter = 1
    while ($counter <= $P1num[1])
	echo $topupParamPattern:s/PE/1/ >>! $topupParam
	@ counter = $counter + 1	
    end	
    set counter = 1
    while ($counter <= $P2num[1])
	echo $topupParamPattern:s/PE/-1/ >>! $topupParam
	@ counter = $counter + 1
    end
endif

if ( -e $B0P1P2File && ! -e $B0P1P2Corr) then
    echo running topup $B0P1P2Corr >> $reportFile
    topup --imain=$B0P1P2File --datain=$topupParam --out=$P1P2Topup:r:r:s/_fieldcoef//  --fout=$P1P2Field --iout=$B0P1P2Corr --config=$topupConfig
endif

#################################################
#################################################
# 5. Extract mask from average topup corr B0
#################################################
#################################################

if ( -e $B0P1P2Corr && ! -e $B0P1P2Mask) then
    echo generating brain masks $B0P1P2Mask >> $reportFile
    fslmaths $B0P1P2Corr -Tmean $B0P1P2Avg
    bet $B0P1P2Avg $B0P1P2Mask -m -v
endif

#################################################
#################################################
# 6. Run Eddy
#################################################
#################################################

echo running eddy

if (! -e $indexFile) then
    echo generating $indexFile >> $reportFile
    set spaces = `fslinfo  $P1P2NIFTI | grep dim4`
    set numDWIs = $spaces[2]
    set b0Info = (`cat $B0indexP1P2File` $numDWIs)
    set numBaseline = $#b0Info
    set counter = 0
    set row = 1
    set indexSet = ()
    while ($counter < $numDWIs)
	@ nextrow = $row + 1
	if ($counter ==  $b0Info[$nextrow] ) then
	    @ row = $nextrow
	endif
	set indexSet = ($indexSet $row)
	@ counter = $counter + 1
    end
    echo $indexSet >! $indexFile
endif

if (-e $P1P2NIFTI && -e $B0P1P2Corr && -e $B0P1P2Mask && ! -e $P1P2Corr) then
    echo applying eddy and topup $P1P2Corr:t >> $reportFile
    eddy_openmp --imain=$P1P2NIFTI --mask=$B0P1P2Mask:r:r_mask.nii.gz --acqp=$topupParam --index=$indexFile --bvecs=$P1P2NIFTI:r:r.bvecs --bvals=$P1P2NIFTI:r:r.bvals --topup=$P1P2Topup:r:r:s/_fieldcoef//  --out=$P1P2Corr:r:r --data_is_shelled --verbose --estimate_move_by_susceptibility
endif

#################################################
#################################################
# 7. Make Corr DWI positive only
#################################################
#################################################

echo making noNeg
if (-e $P1P2Corr && ! -e $P1P2CorrNoNeg) then
    echo removing negative values $P1P2CorrNoNeg >> $reportFile
    fslmaths $P1P2Corr -thr 0 $P1P2CorrNoNeg
endif 

#################################################
#################################################
# 8. Run QUAD
#################################################
#################################################

echo running Quad
if (! -e $quadReport && -e $P1P2Corr) then
    echo running quad for $P1P2Corr >> $reportFile
    /nas/longleaf/apps/fsl/6.0.3/fsl/bin/eddy_quad $P1P2Corr:r:r -idx $indexFile  -par $topupParam  -m $B0P1P2Mask:r:r_mask.nii.gz -b $P1P2NIFTI:r:r.bvals
endif

#################################################
#################################################
# 9. Convert to NRRD
#################################################
#################################################

echo conversion back to NRRD
if (! -e $P1P2CorrNRRD && -e $P1P2CorrNoNeg) then
    echo converting back to NRRD $P1P2CorrNRRD >> $reportFile

    $DWIconvertCMD --conversionMode FSLToNrrd --inputBVectors $P1P2NIFTI:r:r.bvecs --inputBValues $P1P2NIFTI:r:r.bvals --fslNIFTIFile $P1P2CorrNoNeg --outputVolume $P1P2CorrNRRD
endif


### 10. Cleanup some huge intermediate files to reduce footprint
#rm $P1P2Corr $P1P2CorrNoNeg $P1dtiprepQC:r_pad.nrrd

###################
###################
exit
###################
###################