#!/bin/csh -f
# xhemireg

# Create a new subject from an existing subject by left-right reversal
# of the volumes. Should be enough to be able to include in 
# make_average_subject. 
#
# surfaces: 
#  Apply xfm: orig, white, pial, smoothwm, inflated
#  Copy to contra hemi: area, area.pial, curv, curv.pial, avg_curv, 
#     inflated.{H,K}, sulc, thickness, jacobian_white, aparc, aparc.a2005s
#  Regenerate: qsphere, sphere, sphere.reg
#  Not sure: before, after, defect, uncorrected, *.nofix
#
# volumes:
#  Dont convert volumes in mri/orig
#

set VERSION = 'xhemireg dev';

set subject = ();
set srchemilist = ();
set surfreg = xhemi.reg
set UseNoNeg = 0
set NoRandomness = 1;
set RngSeed = 1234
set DoReg = 0;
set PrintHelp = 0;
set outdir = ();
set DoAllVol = 0;
set DoComputeTal = 1;
set DoEstimateTal = 0;
set gcaprep = 0
set CoRegThreads = 1

# Not all of the lists below are really needed (orig T1 brain needed by make_ave_vol)
set vollist = (orig T1 brain brainmask norm aseg aparc+aseg)
set vollist = ($vollist nu) 
# Note: T1 is needed by mris_make_average_surface
set ovlist = (curv sulc area volume thickness inflated.H )
set ovlist = ($ovlist area.pial area.mid curv.pial inflated.K)
set surflist = (orig white pial smoothwm sphere)
set surflist = ($surflist inflated)
set annotlist = ()
set annotlist = ($annotlist aparc BA aparc.a2009s aparc.DKTatlas)
set labellist = (cortex)

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

source $FREESURFER_HOME/sources.csh

goto parse_args;
parse_args_return:
goto check_params;
check_params_return:

if($#outdir == 0) then
  set xsubject = $subject/xhemi
  set outdir = $SUBJECTS_DIR/$xsubject
endif
mkdir -p $outdir/mri/transforms $outdir/surf $outdir/label $outdir/scripts 
set outdir = `getfullpath $outdir`
echo "outdir is $outdir" 

cd $SUBJECTS_DIR/$subject

if($DoAllVol) then
  # Get a list of all mgh and mgz volumes in mri
  cd mri
  set vollist = (`ls *.mgz *.mgh | sed 's/.mgz//g' | sed 's/.mgh//g'`);
  cd ..
endif

@ nthhemi = 0;
foreach srchemi ($srchemilist)
  @ nthhemi = $nthhemi + 1;
  if($srchemi == lh) set trghemi = rh;
  if($srchemi == rh) set trghemi = lh;

  set LF = $outdir/xhemireg.$srchemi.log
  if(-e $LF) mv $LF $LF.bak
  echo "Logfile is $LF"

  echo $0 | tee -a $LF
  echo $cmdargs | tee -a $LF
  echo $VERSION | tee -a $LF
  date | tee -a $LF
  uname -a | tee -a $LF
  pwd | tee -a $LF
  echo "setenv SUBJECTS_DIR $SUBJECTS_DIR "| tee -a $LF

  # Convert volumes, create LR registration
  if($nthhemi == 1) then
    # Create LR registration
    date | tee -a $LF

    # Determine which dimension gets reversed
    set xsign = "+1"
    set ysign = "+1"
    set zsign = "+1"
    set orientation = `mri_info --orientation mri/orig.mgz`
    set o1 = `echo $orientation | cut -c 1`
    set o2 = `echo $orientation | cut -c 2`
    set o3 = `echo $orientation | cut -c 3`
    if($o1 == L || $o1 == R) set xsign = "-1"
    if($o2 == L || $o2 == R) set zsign = "-1" # 2=z in tkreg
    if($o3 == L || $o3 == R) set ysign = "-1" # 3=y in tkreg

    echo orientation: $orientation $o1  $o2  $o3 | tee -a $LF

    # This is the same as --left-right-reverse-pix
    set reg = $outdir/lrrev.register.dat
    rm -f $reg
    echo "subject-unknown" >> $reg
    echo "1" >> $reg
    echo "1" >> $reg
    echo "1" >> $reg
    echo "$xsign 0 0 1" >> $reg
    echo " 0 $ysign 0 0" >> $reg
    echo " 0 0 $zsign 0" >> $reg
    echo " 0 0 0 1" >> $reg
    echo "round" >> $reg

    set regpure = $outdir/lrrev.pure.register.dat
    rm -f $regpure
    echo "subject-unknown" >> $regpure
    echo "1" >> $regpure
    echo "1" >> $regpure
    echo "1" >> $regpure
    echo "$xsign 0 0 0" >> $regpure
    echo " 0 $ysign 0 0" >> $regpure
    echo " 0 0 $zsign 0" >> $regpure
    echo " 0 0 0 1" >> $regpure
    echo "round" >> $regpure

    # Convert volumes. aseg will be wrong -- need to swap labels too
    date | tee -a $LF
    foreach vol ($vollist )
      echo "vol $vol --------------" | tee -a $LF
      set invol = mri/$vol.mgz
      if(! -e $invol) then
        set invol = mri/$vol.mgh
        if(! -e $invol) then
          echo "WARNING: cannot find $vol" | tee -a $LF
          echo " ... skipping ... "| tee -a $LF
          continue
          #exit 1;
        endif
      endif
      set cmd = (mri_vol2vol --mov $invol --targ $invol \
        --reg $reg --o $outdir/mri/$vol.mgz --keep-precision \
        --no-save-reg)
      echo $cmd  | tee -a $LF
      $cmd | tee -a $LF
      if($status) exit 1;

      # Left-Right Reverse label codes for segmentations
      if($vol == aseg || $vol == aparc+aseg || $vol == wmparc || \
         $vol == aparc.a2009s+aseg) then
        set cmd = (mri_convert $outdir/mri/$vol.mgz $outdir/mri/$vol.mgz --left-right-swap-label)
        echo $cmd  | tee -a $LF
        $cmd | tee -a $LF
        if($status) exit 1;
      endif

    end

    # Compute talairach reg in one of two ways
    set xxfm = $SUBJECTS_DIR/$subject/xhemi/mri/transforms/talairach.xfm
    set xlta = $SUBJECTS_DIR/$subject/xhemi/mri/transforms/talairach.xfm.lta
    # The LTA may be used to init the surfreg
    if($DoComputeTal) then
      # Recompute talairach registration by doing the registration
      set cmd = (recon-all -sb $subject/xhemi -talairach)
      date | tee -a $LF
      echo $cmd  | tee -a $LF
      $cmd | tee -a $LF
      if($status) exit 1;
      # Create the LTA
      set ud = `UpdateNeeded $xlta $xxfm`
      if($ud) then
        set mni305 = $FREESURFER_HOME/average/mni305.cor.mgz
        set cmd = (lta_convert --src $subjdir/mri/orig.mgz --trg $mni305 --inxfm $xxfm \
            --outlta $xlta --subject fsaverage --ltavox2vox)
        echo $cmd |& tee -a $LF 
        $cmd |& tee -a $LF
        if($status) exit 1
      endif

    endif

    if($DoEstimateTal) then 
      # Compute the reg based on the unflipped; then create xfm from
      # it.  This saves time, but it also better than recomputing it
      # when a non-recon-all process was needed to compute the
      # unflipped xfm.
      set subjdir = $SUBJECTS_DIR/$subject
      set xfm = $subjdir/mri/transforms/talairach.xfm # unflipped
      set lta = $subjdir/mri/transforms/talairach.xfm.lta
      set ud = `UpdateNeeded $lta $xfm`
      if($ud) then
        # Create the xfm.lta in the unflipped, if needed. 
        set mni305 = $FREESURFER_HOME/average/mni305.cor.mgz
        set cmd = (lta_convert --src $subjdir/mri/orig.mgz --trg $mni305 --inxfm $xfm \
            --outlta $lta --subject fsaverage --ltavox2vox)
        echo $cmd |& tee -a $LF 
        $cmd |& tee -a $LF
        if($status) exit 1
      endif
      set ud = (`UpdateNeeded $xlta $lta`)
      if($ud) then
        set cmd = (mri_coreg --lrrev $lta $xlta)
        echo $cmd |& tee -a $LF 
        $cmd |& tee -a $LF
        if($status) exit 1
        set cmd = (lta_convert --inlta $xlta --outmni $xxfm)
        echo $cmd |& tee -a $LF 
        $cmd |& tee -a $LF
        if($status) exit 1
      endif
    endif

    if($gcaprep) then
      # Create a new xfm from passed xfm
      # Convert passed xfm to lta
      set taltemplate = $FREESURFER_HOME/average/mni305.cor.mgz
      mkdir -p $outdir/mri/transforms
      set cmd = (lta_convert --trg $taltemplate --src mri/orig.mgz \
        --inmni mri/transforms/$gcaprepxfm \
        --outlta $outdir/mri/transforms/norev.$gcaprepxfm.lta)
      echo $cmd  | tee -a $LF
      $cmd | tee -a $LF
      if($status) exit 1;
      # Register the lrrev back to the norev. Note this is different than
      # the lrrev.register.dat above as this is an actual registration
      # of intensities
      set cmd = (mri_coreg --mov $outdir/mri/brainmask.mgz --targ mri/brainmask.mgz \
        --reg $outdir/mri/transforms/reg.lrrev-to-norev.lta --threads $CoRegThreads)
      echo $cmd  | tee -a $LF
      $cmd | tee -a $LF
      if($status) exit 1;
      # Concatenate the two to get the lrev reg to the taltemplate
      set cmd = (mri_concatenate_lta $outdir/mri/transforms/reg.lrrev-to-norev.lta \
        $outdir/mri/transforms/norev.$gcaprepxfm.lta \
        $outdir/mri/transforms/$gcaprepxfm.lta)
      echo $cmd  | tee -a $LF
      $cmd | tee -a $LF
      if($status) exit 1;
      # Convert registration to xfm
      set cmd = (lta_convert --inlta  $outdir/mri/transforms/$gcaprepxfm.lta \
        --outmni  $outdir/mri/transforms/$gcaprepxfm)
      echo $cmd  | tee -a $LF
      $cmd | tee -a $LF
      if($status) exit 1;
      # This command will need to be modified if an output folder has not been spec
      echo "To test new registration run:"  | tee -a $LF
      echo " tkregister2 --fstal --s `basename $outdir`"  | tee -a $LF
    endif

  endif # nthhemi == 1

  # Convert Surfaces (apply LR registration)
  foreach surf ($surflist)
    date | tee -a $LF
    set cmd = (mri_surf2surf --s $subject --sval-xyz $surf \
     --hemi $srchemi --tval-xyz $outdir/mri/$vollist[1].mgz \
     --tval $outdir/surf/$trghemi.$surf)
    if($surf == inflated || $surf == sphere) then
      set cmd = ($cmd --reg $regpure)
    else
      set cmd = ($cmd --reg $reg)
    endif
    echo $cmd  | tee -a $LF
    $cmd | tee -a $LF
    if($status) exit 1;
  end

  # Copy overlays (curvs, thicknesses, etc)
  foreach ov ($ovlist )
    if(! -e surf/$srchemi.$ov) continue;
    set cmd = (cp surf/$srchemi.$ov $outdir/surf/$trghemi.$ov)
    echo $cmd >> $LF
    $cmd | tee -a $LF
    if($status) exit 1;
  end

  # Copy labels (no need to convert, just swap lh rh)
  # May need to convert if xyz ever used.
  foreach label ($labellist )
    set srclabel = label/$srchemi.$label.label
    if(! -e $srclabel) continue;
    set trglabel = xhemi/label/$trghemi.$label.label
    cp $srclabel $trglabel 
    if($status) exit 1;
  end

  # Copy annotations (just swap lh rh)
  foreach annot ($annotlist)
    set fname = label/$srchemi.$annot.annot
    if(! -e $fname) continue;
    set cmd = (cp $fname $outdir/label/$trghemi.$annot.annot)
    echo $cmd  | tee -a $LF
    $cmd | tee -a $LF
    if($status) exit 1;
  end

  # Perform Surface-based registration to fsaverage
  if($DoReg) then
    set cmd = (recon-all -sb $xsubject -hemi $trghemi -surfreg)
    date | tee -a $LF
    echo $cmd  | tee -a $LF
    $cmd | tee -a $LF
    if($status) exit 1;
  endif

  echo "done hemi $srchemi"  | tee -a $LF

end # Loop over srchemi

#--------------------------------------------------------------
echo "" | tee -a $LF
date | tee -a $LF
echo "xhemireg done"  | tee -a $LF

exit 0

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

  set flag = $argv[1]; shift;

  switch($flag)

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

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

    case "--lh":
      set srchemilist = lh;
      breaksw

    case "--rh":
      set srchemilist = rh;
      breaksw

    case "--all-vol":
      set DoAllVol = 1;
      set DoReg = 0;
      breaksw

    case "--zilles":
      set vollist = (orig T1 brain brainmask nu norm)
      set DoReg = 0;
      breaksw

    case "--gcaprep":
      # This is for when creating a symmetric gca
      # Need to create mri/tranforms
      if( $#argv < 2) goto arg2err;
      set gcaprep = 1
      set gcaprepseg = $argv[1]; shift; # this should not have an extension
      set gcaprepxfm = $argv[1]; shift; # not sure what to do with this yet
      set vollist = (orig T1 nu brainmask $gcaprepseg)
      set ovlist = ()
      set surflist = ()
      set DoTal = 0;
      breaksw

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

    case "--avgerage-subject":
    case "--avgeragesubject":
    case "--avgsubject":
      set vollist = (orig T1 brain brainmask aseg aparc+aseg)
      set DoReg = 0;
      breaksw

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

    case "--noreg":
    case "--no-reg":
      set DoReg = 0;
      breaksw

    case "--reg":
      set DoReg = 1;
      breaksw

    case "--tal-compute":
      set DoComputeTal = 1;
      breaksw
    case "--no-tal":
    case "--no-tal-compute":
      set DoComputeTal = 0;
      breaksw

    case "--tal-estimate":
      set DoEstimateTal = 1;
      set DoComputeTal = 0;
      breaksw
    case "--no-tal-estimate":
      set DoEstimateTal = 0;
      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($#subject == 0) then
  echo "ERROR: must spec a subject id"
  exit 1;
endif

if(! -e $SUBJECTS_DIR/$subject) then
  echo "ERROR: cannot find $subject"
  exit 1;
endif

if($#srchemilist == 0) set srchemilist = (lh rh)

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 "xhemireg"
  echo ""
  echo "   --s subject"
  echo ""
  echo "   --lh : srchemi = lh (map from left to right) (default is to do lh and rh)"
  echo "   --reg       : perform reg to create sphere.reg"
  echo "   --rh : srchemi = rh (map from right to left) "
  echo "   --o outdir  : output dir (subject/xhemireg.hemi by def)"
  echo "   --tal-compute    : recompute talairach registration (default)"
  echo "   --no-tal-compute  : do not perform talairach registration"
  echo "   --tal-estimate    : compute estimate of talairach registration from unflipped reg"
  echo "   --no-tal-estimate "
  echo "   --gcaprep seg xfm : for training symmetrical GCA atlases (see gcaprepone)"
  echo "   --threads nthreads: only applies with --gcaprep"
  echo "   --version   : print version and exit"
  echo "   --help      : print help and exit"
  echo ""

  if($PrintHelp) \
  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