#!/bin/tcsh -f

#########################################################################
#
# Auto grad-flip tester (v2.1, Sep 2015)
#   by PA Taylor (UCT, AIMS)
#
# Simple script to test what 'flip', if any, should likely be
# performed for a data set when using '1dDW_Grad_o_Mat.
#
#
# + USING
# The script takes two arguments to run:
# 
# $ @GradFlipTest  'STRING_CMD'  DWI_FILE
#
# where:
#    STRING_CMD :is the string command you would use in 1dDW_Grad_o_Mat
#                to input the gradient vector file, for example:
#                '-in_grad_cols GRAD_FILE'  (and make sure to put quotes
#                around the entry).  If you need to include other 
#                options, such as 'in_bvals FILE' or something, you can
#                do that in this single string command, as well.
#
#   DWI_FILE    :the name of the DWI set to be processed and possibly
#                averaged together in parts at the same time as the grad
#                file.
#
# + OUTPUT
#   Mostly, this code will recommend using either 'No flip', or one of
#   the {-flip_x|-flip_y|-flip_z}.  It will produce a directory called 
#   '_tmp_TESTFLIP/' to store temporary files, of which there are many.
#   It will also prompt you the user to visually check the tract results
#   as well with some simple example scripts.
#
########################################################################

# v2.1, Sep,2015:  linear DTI fits, for speed
# v2.2, Sep,2015:  echo edu to print commands

# trap for help or bad input   9 Feb 2015 [rickr]
if ( $#argv != 2 ) then
   echo "usage `basename $0` 'STRING_CMD' DWI_FILE"
   exit
endif

set in_CMD = "$1"
set in_DSET = "$2"

set FLIPS = ( "NO flip" '-flip_x' "-flip_y" "-flip_z" )
set FNOMS = ( "F0" "FX" "FY" "FZ" )

set dir = "_tmp_TESTFLIP"
set ALL_ARR =  ( )
set SUMA_CALL = ""

set min_tr_len = 30

# make a temporary output directory
if( -d "$dir" ) then
    echo "Directory '$dir' exists already."
else
    printf "\n\nMaking the temporary directory '$dir' for the files."
    mkdir $dir
endif

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
#                 Start the looping
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

printf "\n\n\n++ Starting to test...\n"

set i = '1'

# Setup names
set out_GRAD  = "$dir/$FNOMS[$i]_GRAD.dat"
set out_DSET  = "$dir/$FNOMS[$i]_DSET.nii.gz"
set out_DT    = "$dir/$FNOMS[$i]_DT"
set out_TRK   = "$dir/$FNOMS[$i]_WB"
set out_GRID  = "$dir/$FNOMS[$i]_WB_000.grid"
set out_TRACT = "$dir/$FNOMS[$i]_WB_000.niml.tract"
set out_mask  = "$dir/mask.nii.gz"

printf "\n\n++ Grad_o_Matting (just needed once).\n\n\n"

# Only do Grad_o_Mat with no flipping-- only fit tensors once, and
# then flip eigenvectors individually
1dDW_Grad_o_Mat                          \
    -echo_edu                            \
    $in_CMD                              \
    -proc_dset      $in_DSET             \
    -pref_dset      $out_DSET            \
    -out_grad_cols  $out_GRAD            \
    -overwrite

# mask based on 0th brick of 1st set
3dAutomask                               \
    -echo_edu                            \
    -overwrite                           \
    -prefix $out_mask                    \
    $out_DSET'[0]'

printf "\n\n++ Tensor fitting (just needed once).\n\n\n"

# Tensor fitting of non-flipped
3dDWItoDT                                \
    -echo_edu                            \
    -eigs -linear -sep_dsets             \
    -mask $out_mask                      \
    -prefix $out_DT                      \
    $out_GRAD                            \
    $out_DSET                            \
    -overwrite


set ext = `3dinfo -av_space ${out_DT}_FA*HEAD`


printf "\n\n++ Make the new eigenvectors for each flip.\n\n\n"


# Make other flips
foreach j ( `seq 1 1 3` )

    # flip_x (i.e., 0th brick in V? files)
    set out_DTx = "$dir/$FNOMS[2]_DT_V${j}"
    3dcalc                                  \
        -a ${out_DT}_V${j}$ext              \
        -expr 'a*(-equals(l,0)+equals(l,1)+equals(l,2))'  \
        -prefix $out_DTx                    \
        -overwrite
        
    set out_DTy = "$dir/$FNOMS[3]_DT_V${j}"
    3dcalc                                  \
        -a ${out_DT}_V${j}$ext              \
        -expr 'a*(equals(l,0)-equals(l,1)+equals(l,2))'  \
        -prefix $out_DTy                    \
        -overwrite

    set out_DTz = "$dir/$FNOMS[4]_DT_V${j}"
    3dcalc                                  \
        -a ${out_DT}_V${j}$ext              \
        -expr 'a*(equals(l,0)+equals(l,1)-equals(l,2))'  \
        -prefix $out_DTz                    \
        -overwrite

end

printf "\n\n++ Now for some tracking with each set of eigenvectors.\n\n"

foreach j ( `seq 1 1 4` )

    set out_DTx_niml = "$dir/$FNOMS[$j].niml.opts"
    # Clear file first, then fill
    printf '' > $out_DTx_niml
    printf "<DTIFILE_opts\n" >> $out_DTx_niml
    printf "dti_V1='$dir/$FNOMS[$j]_DT_V1$ext'\n" >> $out_DTx_niml
    printf "dti_V2='$dir/$FNOMS[$j]_DT_V2$ext'\n" >> $out_DTx_niml
    printf "dti_V3='$dir/$FNOMS[$j]_DT_V3$ext'\n" >> $out_DTx_niml
    printf "dti_FA='${out_DT}_FA$ext'\n" >> $out_DTx_niml
    printf "dti_MD='${out_DT}_MD$ext'\n" >> $out_DTx_niml
    printf "dti_L1='${out_DT}_L1$ext' />\n" >> $out_DTx_niml

    # long tracts
    set out_TRK   = "$dir/$FNOMS[$j]_WB"
    3dTrackID -mode DET                      \
        -echo_edu                            \
        -alg_Nseed_X 1                       \
        -alg_Nseed_Y 1                       \
        -alg_Nseed_Z 1                       \
        -alg_Thresh_Len $min_tr_len          \
        -no_indipair_out                     \
        -logic OR                            \
        -dti_list $out_DTx_niml              \
        -netrois $out_mask                   \
        -prefix $out_TRK                     \
        -overwrite

end


########################################################################
#                           Summary stuff
########################################################################

foreach i ( `seq 1 1 4` )

    # always underlay single FA, which wouldn't change
    set out_DT    = "$dir/$FNOMS[1]_DT"    

    set out_TRK   = "$dir/$FNOMS[$i]_WB"
    set out_GRID  = "$dir/$FNOMS[$i]_WB_000.grid"
    set out_TRACT = "$dir/$FNOMS[$i]_WB_000.niml.tract"

    # grep for the numbers of tracts in each case
    set HasNT = `grep -A 1 "\# NT" $out_GRID`
    printf "Checking GRID file: $out_GRID\n"

    set ALL_ARR = ( $ALL_ARR "$HasNT[3]" )

    # for building output command
    set SUMA_CALL = "${SUMA_CALL}# For '$FLIPS[$i]':\n"
    set SUMA_CALL = "${SUMA_CALL}suma -vol ./${out_DT}_FA*HEAD -tract ./$out_TRACT &\n"

end

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
#                      Find the max of the tracts
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

set L = `printf $#ALL_ARR`
set MaxI = "1"
set MaxV = "$ALL_ARR[$MaxI]"

printf "\n"
foreach i ( `seq 1 1 $L` )

    if ( "$ALL_ARR[$i]" > "$MaxV" ) then
        set MaxI = "$i"
        set MaxV = "$ALL_ARR[$MaxI]"
    endif

    printf "\n The number of tracks for '$FLIPS[$i]' was\t $ALL_ARR[$i]"

end

printf "\n\n --> THEREFORE, I guess that the best flip is: '$FLIPS[$MaxI]'\n\n\n"
printf "All results are in directory './$dir/'\n"
printf "You may check the tracking results to verify, e.g.:\n"
printf "$SUMA_CALL\n"
printf "DONE.\n\n"