#!/bin/tcsh -f
# samseg-long - sources
if(-e $FREESURFER_HOME/sources.csh) then
  source $FREESURFER_HOME/sources.csh
endif

set VERSION = 'samseg-long dev';
set scriptname = `basename $0`

set outdir = ();
set inputlist = ()
set DoMC = ()
setenv OMP_NUM_THREADS 1
set ForceUpdate = 0;
set tmpdir = ();
set cleanup = 1;
set LF = ();
set SaveWarp = 1;
set SavePosteriors = 0
set DoMRF = 0;
set mrfgca = $FREESURFER_HOME/average/samseg/samseg.talairach.m3z.mrf.gca

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
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`

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

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

# Set up log file
if($#LF == 0) set LF = $outdir/log/samseg-long.log
if($LF != /dev/null) rm -f $LF
echo "Log file for samseg-long" >> $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
echo "pid $$" | tee -a $LF
if($?PBS_JOBID) then
  echo "pbsjob $PBS_JOBID"  >> $LF
endif

#========================================================
if($DoMC) then
  mkdir -p $outdir/mc

  # set names for ltas and output volumes
  @ nth = 0
  set RunLtas = ()
  set RunIscales = ()
  set MCVolumes = ()
  foreach input ($inputlist)
    @ nth = $nth + 1
    set nthstr = `printf %03d $nth`
    set nthlta = $outdir/mc/input$nthstr.lta
    set RunLtas=($RunLtas $nthlta)
    set nthiscale = $outdir/mc/input$nthstr-iscale.txt
    set RunIscales=($RunIscales $nthiscale)
    set mcvol = $outdir/mc/mc.$nthstr.mgz
    set MCVolumes = ($MCVolumes $mcvol)
  end

  set template = $outdir/mc/mctemplate.mgz
  set cmd = (mri_robust_template --mov $inputlist --average 1 --template $template \
    --satit  --inittp 1 --fixtp  --noit   --iscale  --iscaleout $RunIscales \
    --subsample 200 --lta $RunLtas --mapmov $MCVolumes)
  echo $cmd | tee -a $LF
  set ud = `UpdateNeeded $template $inputlist`
  if($ud || $ForceUpdate) then
    fs_time $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    echo "\nMotion correction does not need to be updated\n" |& tee -a $LF
  endif
else
  echo "Not performing motion correction" |& tee -a $LF
  set MCVolumes = ($inputlist)
endif

mkdir -p $outdir/inputs
pushd $outdir/inputs
@ nth = 0;
set inputlist2 = ()
foreach invol ($MCVolumes)
  @ nth = $nth + 1
  set nthstr = `printf %03d $nth`
  # filename cannot have a dot in it
  set volname = input_$nthstr.mgz 
  set ext = `fname2ext $invol`
  if($ext != mgz) then
    set ud = `UpdateNeeded $volname $invol`
    if($ud || $ForceUpdate) then
      set cmd = (mri_convert $invol $volname)
      echo $cmd | tee -a $LF
      $cmd |& tee -a $LF
      if($status) goto error_exit;
    else
      echo "Conversion to mgz does not need to be updated"
    endif
  else
    if($DoMC) then
      set mcvol = ../mc/mc.$nthstr.mgz
      set cmd = (ln -sf $mcvol $volname)  
    else
      set cmd = (ln -sf $invol $volname)  
    endif
    echo $cmd | tee -a $LF    
    $cmd |& tee -a $LF    
    if($status) goto error_exit;
  endif
  set inputlist2 = ($inputlist2 -t $outdir/inputs/$volname)
end
popd

echo "\n\n" | tee -a $LF

set cmd = (run_samseg_long $inputlist2 -o $outdir)
if($SaveWarp) set cmd = ($cmd --save-warp)
if($SavePosteriors) set cmd = ($cmd --save-posteriors)

echo $cmd | tee -a $LF

fs_time $cmd |& tee -a $LF
if($status) goto error_exit;

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

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

# Done
echo "\n\n" | tee -a $LF
set tSecEnd = `date '+%s'`;
@ tSecRun = $tSecEnd - $tSecStart;
set tRunMin = `echo $tSecRun/50|bc -l`
set tRunMin = `printf %5.2f $tRunMin`
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 "Samseg-Long-Run-Time-Sec $tSecRun" |& tee -a $LF
echo "Samseg-Long-Run-Time-Min $tRunMin" |& tee -a $LF
echo "Samseg-Long-Run-Time-Hours $tRunHours" |& tee -a $LF
echo " " |& tee -a $LF
echo "samseg-long 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 "--o":
      if($#argv < 1) goto arg1err;
      set outdir = $argv[1]; shift;
      breaksw

    case "--i":
      if($#argv < 1) goto arg1err;
      if(! -e $argv[1]) then
        echo "ERROR: cannot find $argv[1]"
        exit 1;
      endif
      set inputlist = ($inputlist $argv[1]); shift;
      breaksw

    case "--mc":
      set DoMC = 1;
      breaksw

    case "--no-mc":
    case "--nomc":
      set DoMC = 0;
      breaksw

    case "--threads":
      if($#argv < 1) goto arg1err;
      setenv OMP_NUM_THREADS $argv[1]; shift
      breaksw

    case "--save-posteriors":
      set SavePosteriors = 1;
      breaksw

    case "--force-update":
      set ForceUpdate = 1
      breaksw

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

    case "--nolog":
    case "--no-log":
      set LF = /dev/null
      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:
      echo ERROR: Flag $flag unrecognized. 
      echo $cmdline
      exit 1
      breaksw
  endsw

end

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

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

if($#outdir == 0) then
  echo "ERROR: must spec outdir"
  exit 1;
endif
if($#inputlist < 2) then
  echo "ERROR: must spec at least two inputs"
  exit 1;
endif
if($#DoMC == 0) then
  echo "ERROR: must spec either --mc or --no-mc"
  exit 1;
endif

set inputlisttmp = ()
foreach f ($inputlist)
  set ftmp = `getfullpath $f`
  set inputlisttmp = ($inputlisttmp $ftmp)
end
set inputlist = ($inputlisttmp)

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 ""
  echo "samseg-long -o outputdir "
  echo " --i input1 --i input2 ..."
  echo " --mc, --no-mc : align all inputs using robust register (choose one)"
  echo " --threads nthreads"
  echo "  --save-posteriors : save posterior probs"
  echo " --force-update"
  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

For longitudinal samseg analysis. All inputs must be a single modality. If
you are going to run surface analysis, then they all must be T1w.

samseg-long --i timepoint1.mgz --i timepoint2.mgz --mc --o samseglongdir

This will create the samseglongdir folder. This folder will have a base
folder and time point folders called tp001, tp002, etc. Each of these is a 
samseg output folder.

Create a base subject by transfering all the base samseg analysis to the
recon-all base subject. The name of the base subject can be anything:
   samseg2recon --base --s subject_base  --samseg samseglongdir
Now run recon-all on the base
   recon-all -s subject_base -autorecon2-samseg -autorecon3

Create a subject for each time point by transfering all the samseg analysis 
from that TP to the subject for that TP. The subject name can be anything
   samseg2recon --long 2 --s long.tp002  --samseg samseglongdir
Then run
   recon-all -long-samseg subject_base long.tp002 -autorecon2-samseg -autorecon3