#!/bin/tcsh -f
# sratio - signed ratio
#   +A/B if A>B
#   -B/A if B>A

#
# sratio
#
# Script to computed a voxel-wise signed ratio. 
# Used to use fslmaths, but does not anymore
#
# Original Author: Doug Greve
#
# Copyright © 2021 The General Hospital Corporation (Boston, MA) "MGH"
#
# Terms and conditions for use, reproduction, distribution and contribution
# are found in the 'FreeSurfer Software License Agreement' contained
# in the file 'LICENSE' found in the FreeSurfer distribution, and here:
#
# https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferSoftwareLicense
#
# Reporting: freesurfer@nmr.mgh.harvard.edu
#
#

setenv FSLOUTPUTTYPE NIFTI
set tmpdir = ();
set VERSION = 'sratio dev';
set cleanup = 1;
set LF = ();
set A = ();
set B = ();
set sAdivB = ();
set DoAbs = 0;
set MaskThresh = ()

set inputargs = ($argv);
set PrintHelp = 0;
if($#argv == 0) goto usage_exit;
set n = `echo $argv | grep -e -help | wc -l` 
if($n != 0) then
  set PrintHelp = 1;
  goto usage_exit;
endif
set n = `echo $argv | grep -e -version | wc -l` 
if($n != 0) then
  echo $VERSION
  exit 0;
endif

source $FREESURFER_HOME/sources.csh

goto parse_args;
parse_args_return:
goto check_params;
check_params_return:

set StartTime = `date`;
set tSecStart = `date '+%s'`;
set year  = `date +%Y`
set month = `date +%m`
set day   = `date +%d`
set hour   = `date +%H`
set min    = `date +%M`

set outdir = `dirname $sAdivB`

mkdir -p $outdir
pushd $outdir > /dev/null
set outdir = `pwd`;
popd > /dev/null

if($#tmpdir == 0) then
  set tmpdir = `fs_temp_dir --scratch`
endif
mkdir -p $tmpdir

# Set up log file
set LF = /dev/null
echo "Log file for sratio" >> $LF
date  | tee -a $LF
echo "" | tee -a $LF
echo "setenv SUBJECTS_DIR $SUBJECTS_DIR" | tee -a $LF
echo "cd `pwd`"  | tee -a $LF
echo $0 $inputargs | tee -a $LF
echo "" | tee -a $LF
cat $FREESURFER_HOME/build-stamp.txt | tee -a $LF
echo $VERSION | tee -a $LF
uname -a  | tee -a $LF
if($?PBS_JOBID) then
  echo "pbsjob $PBS_JOBID"  >> $LF
endif

#========================================================
mkdir -p $tmpdir

if($DoAbs) then
  set A2 = $tmpdir/A.nii
  set cmd = (mri_concat $A --abs --o $A2)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) goto error_exit;
  set A = $A2
  set B2 = $tmpdir/B.nii
  set cmd = (mri_concat $B --abs --o $B2)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) goto error_exit;
  set B = $B2
endif

# A > B
set AgtB = $tmpdir/AgtB.nii
set cmd = (fscalc $A -sub $B -o $AgtB)
echo $cmd
$cmd
if($status) exit 1;

# B > A (A < B)
set BgtA = $tmpdir/BgtA.nii
set cmd = (fscalc $B -sub $A -o $BgtA)
echo $cmd
$cmd
if($status) exit 1;

# A/B
set AdivB = $tmpdir/AdivB.nii
set cmd = (fscalc $A -div $B -o $AdivB)
echo $cmd
$cmd
if($status) exit 1;

# -B/A
set BdivA = $tmpdir/BdivA.nii
set cmd = (fscalc $B -div $A -mul -1 -o $BdivA) # signed by -1
echo $cmd
$cmd
if($status) exit 1;

# A > B mask
set AgtBmask = $tmpdir/AgtB.mask.nii
set cmd = (mri_binarize --i $AgtB --min 0 --o $AgtBmask)
echo $cmd
$cmd
if($status) exit 1;

# Anum = A/B where A > B
set Anum = $tmpdir/Anum.nii
set cmd = (fscalc $AdivB -mul $AgtBmask -o $Anum)
echo $cmd
$cmd
if($status) exit 1;

# B > A mask
set BgtAmask = $tmpdir/BgtA.mask.nii
set cmd = (mri_binarize --i $BgtA --min 0 --o $BgtAmask)
echo $cmd
$cmd
if($status) exit 1;

# Bnum = -B/A where  B > A
set Bnum = $tmpdir/Bnum.nii
set cmd = (fscalc $BdivA -mul $BgtAmask -o $Bnum)
echo $cmd
$cmd
if($status) exit 1;

# Now just add Anum and Bnum
set sAdivBtmp = $tmpdir/sAdivB.nii
set cmd = (fscalc $Anum -add $Bnum -o $sAdivBtmp)
echo $cmd
$cmd
if($status) exit 1;

if($#MaskThresh) then
  # Mask and convert to output
  set mask = $tmpdir/mask.nii
  set cmd = (mri_concat --abs $A $B --max --o $mask)
  echo $cmd
  $cmd
  if($status) exit 1;
  set cmd = (mri_binarize --i $mask --min $MaskThresh --o $mask)
  echo $cmd
  $cmd
  if($status) exit 1;
  set cmd = (mri_mask $sAdivBtmp $mask $sAdivB)
  echo $cmd
  $cmd
  if($status) exit 1;
else
  # Just convert to output
  set cmd = (mri_convert $sAdivBtmp $sAdivB)
  echo $cmd
  $cmd
  if($status) exit 1;
endif


#========================================================

# Cleanup
if($cleanup) rm -rf $tmpdir

# Done
echo " " |& tee -a $LF
set tSecEnd = `date '+%s'`;
@ tSecRun = $tSecEnd - $tSecStart;
set tRunHours = `echo $tSecRun/3600|bc -l`
set tRunHours = `printf %5.2f $tRunHours`
echo "Started at $StartTime " |& tee -a $LF
echo "Ended   at `date`" |& tee -a $LF
echo "Sratio-Run-Time-Sec $tSecRun" |& tee -a $LF
echo "Sratio-Run-Time-Hours $tRunHours" |& tee -a $LF
echo " " |& tee -a $LF
echo "tkmedit -f $A -aux $B -fminmax 1.01 1.2 -ov $sAdivB" |& tee -a $LF
echo "" |& tee -a $LF
echo "sratio Done" |& tee -a $LF
exit 0

###############################################

############--------------##################
error_exit:
echo "ERROR:"

exit 1;
###############################################

############--------------##################
parse_args:
set cmdline = ($argv);
while( $#argv != 0 )

  set flag = $argv[1]; shift;
  
  switch($flag)

    case "--abs":
      set DoAbs = 1;
      breaksw

    case "--mask-thresh":
      if($#argv < 1) goto arg1err;
      set MaskThresh = $argv[1]; shift;
      breaksw

    case "--tmp":
    case "--tmpdir":
      if($#argv < 1) goto arg1err;
      set tmpdir = $argv[1]; shift;
      set cleanup = 0;
      breaksw

    case "--nocleanup":
      set cleanup = 0;
      breaksw

    case "--cleanup":
      set cleanup = 1;
      breaksw

    case "--debug":
      set verbose = 1;
      set echo = 1;
      breaksw

    default:
      if($#A == 0) then
        set A = $flag;
        if(! -e $A) then
          echo "ERROR: cannot find $A"
          exit 1;
        endif
      else if( $#B == 0 ) then
        set B = $flag;
        if(! -e $B) then
          echo "ERROR: cannot find $B"
          exit 1;
        endif
      else if($#sAdivB == 0) then 
        set sAdivB = $flag;
      else
        echo ERROR: Flag $flag unrecognized. 
        echo $cmdline
        exit 1
      endif
      breaksw
  endsw

end

goto parse_args_return;
############--------------##################

############--------------##################
check_params:

if($#A == 0) then
  echo "ERROR: must spec A"
  exit 1;
endif
if($#B == 0) then
  echo "ERROR: must spec B"
  exit 1;
endif
if($#sAdivB == 0) then
  echo "ERROR: must spec output"
  exit 1;
endif

goto check_params_return;
############--------------##################

############--------------##################
arg1err:
  echo "ERROR: flag $flag requires one argument"
  exit 1
############--------------##################
arg2err:
  echo "ERROR: flag $flag requires two arguments"
  exit 1
############--------------##################

############--------------##################
usage_exit:
  echo "sratio options A B sAdivB"
  echo "  A/B if A>B, -B/A if B>A"
  echo " Options:"
  echo "   --abs : compute absolute value of both A and B before sratio"
  echo "   --mask-thresh thresh : theshold based on max(abs(A),abs(B))>thresh"
  echo ""

  if(! $PrintHelp) exit 1;
  echo $VERSION
  cat $0 | awk 'BEGIN{prt=0}{if(prt) print $0; if($1 == "BEGINHELP") prt = 1 }'
exit 1;

#---- Everything below here is printed out as part of help -----#
BEGINHELP