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

set VERSION = 'longmc dev';
set scriptname = `basename $0`

set tpNid = (); # cross-sectional subject name
set longbaseid = (); # lnog base  subject name
set subjid = (); # eventually, subjid = $tpNid.long.$longbaseid
set DoConf2Hires = 0;

set tmpdir = ();
set cleanup = 1;
set LF = ();

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`

#if($#tmpdir == 0) then
#  if(-dw /scratch)   set tmpdir = /scratch/tmpdir.longmc.$$
#  if(! -dw /scratch) set tmpdir = $outdir/tmpdir.longmc.$$
#endif
#mkdir -p $tmpdir

set subjdir = $SUBJECTS_DIR/$subjid
set mdir = $subjdir/mri
mkdir -p $mdir/orig $mdir/transforms $subjdir/scripts

# Set up log file
if($#LF == 0) set LF = $subjdir/scripts/longmc.log
if($LF != /dev/null) rm -f $LF
echo "Log file for longmc" >> $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

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

# Create the orig.mgz and rawavg.mgz from the cross by concatenating
# LTAs to the base and resampling. Note that T2 and FLAIR do not
# need to be handled here in either the normal stream or conf2hires

# When there are multiple T1 inputs (eg, 001.mgz, 002.mgz), 
# each will have an LTA (and -iscale.txt). In cross, MC is done by
# running robust_template using the ???.mgz and ???.lta. Here,
# we create new LTAs by concatenating the existing LTAs with
# the lta to base and then just run robust_template with the
# original ???.mgz and the new LTAs.

# These are the outputs
set origvol = $mdir/orig.mgz;
set rawvol = $mdir/rawavg.mgz; 

if($#CrossList == 1) then
  # Single input in cross, directly resample to base space
  set srcvol = $CrossList[1];
  # Create orig.mgz. Note: in cross, this would be created
  # directly from the rawavg. But with one input, srcvol
  # and rawavg are the same in the cross.
  # Always use cubic here?
  set cmd = (mri_convert -at $tpNtobase_regfile \
    -odt uchar -rt cubic $srcvol $origvol)
  echo "\n $cmd \n" |& tee -a $LF 
  $cmd |& tee -a $LF
  if($status) goto error_exit;
  # Create rawavg image. 
  if(! $DoConf2Hires) then
    # Use cubic here? Might not matter
    set cmd = (mri_convert -at $tpNtobase_regfile -rt cubic $srcvol $rawvol)
    echo "\n $cmd \n" |& tee -a $LF 
    $cmd |& tee -a $LF
    if($status) goto error_exit;
  else
    # Normally, rawavg not used for anything, but in conf2hires it is
    # Create rawavg image (not used by the long stream but used by conf2hires)
    # Need to compute the reg that takes rawavg into template space
    # and then change the header, not the pix, to make it in reg with base
    set raw2confreg = $subjdir/mri/transforms/rawavgcross2conf.lta
    set tmpfile = `fs_temp_file --suffix .dat`
    set cmd = (tkregister2_cmdl --noedit --mov ${SUBJECTS_DIR}/${tpNid}/mri/rawavg.mgz \
      --targ ${SUBJECTS_DIR}/${tpNid}/mri/orig.mgz --regheader --reg $tmpfile \
      --ltaout $raw2confreg)
    echo "\n $cmd \n" |& tee -a $LF 
    $cmd |& tee -a $LF
    if($status) goto error_exit;
    rm -f $tmpfile
    set raw2basereg = $subjdir/mri/transforms/rawavgcross2base.lta
    set cmd = (mri_concatenate_lta $raw2confreg $tpNtobase_regfile $raw2basereg)
    echo "\n $cmd \n" |& tee -a $LF 
    $cmd |& tee -a $LF
    if($status) goto error_exit;
    set rawvol  = $subjdir/mri/rawavg.mgz
    set cmd = (mri_vol2vol --mov ${SUBJECTS_DIR}/${tpNid}/mri/rawavg.mgz \
     --targ $SUBJECTS_DIR/$longbaseid/mri/orig.mgz --reg $raw2basereg --no-resample --o $rawvol)
    echo "\n $cmd \n" |& tee -a $LF 
    $cmd |& tee -a $LF
    if($status) goto error_exit;
  endif
else
  # Multiple inputs in cross subject
  set LongLtas = ()
  set LongIscales = ()
  set ConcatLtas = ()
  foreach cross ($CrossList)
    set stem = `fname2stem $cross`
    set stembase = `basename $stem`
    set lta = $stem.lta
    set iscale = $stem-iscale.txt
    # copy over lta and iscale
    set cmd = (cp $lta $iscale $mdir)
    echo "\n $cmd \n" |& tee -a $LF 
    $cmd |& tee -a $LF
    if($status) goto error_exit;
    # Keep in a (somewhat redundant) list
    #set LongLtas    = ($LongLtas    $LongLtas $SUBJECTS_DIR/$tpNid/mri/orig/$stembase.lta) # Not needed
    set LongIscales = ($LongIscales $SUBJECTS_DIR/$tpNid/mri/orig/$stembase-iscale.txt)
    # create a new lta to the base template
    set concatlta = $stem-to-base.lta
    set cmd = (mri_concatenate_lta $lta $tpNtobase_regfile $concatlta)
    echo "\n $cmd \n" |& tee -a $LF 
    $cmd |& tee -a $LF
    if($status) goto error_exit;
    set ConcatLtas = ($ConcatLtas $concatlta)
  end

  # Use mri_robust_template to create the rawavg (as with cross stream)
  # Do this even for conf2hires because it has to be resampled anyway
  set cmd = (mri_robust_template --mov $CrossList --ixforms $ConcatLtas --average 1 \
    --noit --template $rawvol --iscalein $LongIscales)
  echo "\n $cmd \n" |& tee -a $LF 
  $cmd |& tee -a $LF
  if($status) goto error_exit;

  # Convert the rawavg to orig (as with cross stream)
  set cmd = (mri_convert -odt uchar $rawvol $origvol)
  echo "\n $cmd \n" |& tee -a $LF 
  $cmd |& tee -a $LF
  if($status) goto error_exit;

endif

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

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

# Done
echo " " |& 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 "Longmc-Run-Time-Sec $tSecRun" |& tee -a $LF
echo "Longmc-Run-Time-Min $tRunMin" |& tee -a $LF
echo "Longmc-Run-Time-Hours $tRunHours" |& tee -a $LF
echo " " |& tee -a $LF
echo "longmc 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 "-long":
    case "-longitudinal":
      if ( $#argv < 2) goto arg1err;
      # get the cross subject name to use for timepoint
      set tpNid = $argv[1]; shift;
      set tpNid = `basename $tpNid`; # remove trailing /
      # get the subject to use for the base subject
      set longbaseid = $argv[1]; shift;
      set longbaseid = `basename $longbaseid`; # remove trailing /
      # and create subjid to reflect its longitudinal relation to longbaseid
      set subjid = ${tpNid}.long.${longbaseid}
      breaksw

    case "-s":
    case "--s":
      # Override, good for testing, must be put after -long
      if($#argv < 1) goto arg1err;
      set subjid = $argv[1]; shift;
      breaksw

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

    case "--conf2hires":
    case "-conf2hires":
      set DoConf2Hires = 1;
      breaksw
    case "--no-conf2hires":
    case "-no-conf2hires":
      set DoConf2Hires = 0;
      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":
    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($#subjid == 0) then
  echo "ERROR: must spec subject with -long"
  exit 1;
endif
foreach s ($tpNid $longbaseid)
  if(! -e $SUBJECTS_DIR/$s) then
    echo "ERROR: cannot find $s"
    exit 1;
  endif
end

# in order to create orig.mgz in LONG we need at least 001.mgz in CROSS:
#set CrossList = (`ls ${SUBJECTS_DIR}/${tpNid}/mri/orig/[0-9][0-9][0-9].mgz`);
set CrossList = (`find ${SUBJECTS_DIR}/${tpNid}/mri/orig -iname "[0-9][0-9][0-9].mgz" -print`)
if ( $#CrossList == 0) then
  echo "ERROR: no CROSS run data found in ${SUBJECTS_DIR}/${tpNid}/mri/orig. Make sure to" |& tee -a $LF
  echo "have a volume called 001.mgz there." |& tee -a $LF
  echo "If you have a second run of data call it 002.mgz, etc." |& tee -a $LF
  echo "See also: http://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/Conversion" |& tee -a $LF
  exit 1
endif

if($#CrossList > 1) then
  foreach cross ($CrossList)
    set stem = `fname2stem $cross`
    set lta = $stem.lta
    if(! -e $lta) then
      echo "ERROR: cannot find $lta"
      exit 1;
    endif
    # Not sure the iscale is really needed
    set iscale = $stem-iscale.txt
    if(! -e $iscale) then
      echo "ERROR: cannot find $iscale"
      exit 1;
    endif
  end
endif

set tpNtobase_regfile = $SUBJECTS_DIR/$longbaseid/mri/transforms/${tpNid}_to_${longbaseid}.lta
if(! -e $tpNtobase_regfile) then
  echo "ERROR: cannot find $tpNtobase_regfile"
  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 ""
  echo "longmc "
  echo " -long CrossTPname Basename"
  echo " -conf2hires, -no-conf2hires"
  echo " -sd SUBJECTS_DIR"
  echo " -s subjectname : override, must be after -long"
  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

Performes the motion correction step for the longitudinal recon-all stream
when creating the longitudinal timepoint of a subject.