#! /bin/tcsh -f

#
# bbregister
#
# 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
#
#

#
set VERSION = 'bbregister dev';
set inputargs = ($argv);

set subject = ();
set movvol = ();
set intvol = (); # intermediate volume
set frame = ();
set midframe = 0;
set outreg = ();
set outfsl = ();
set outlta = ();
set InitFSL = 0;
set InitSPM = 0;
set InitCoreg = 1;
set InitRR = 0;
set SPMUseNII = 0;
set InitHeader = 0;
set InitReg = ();
set InitRegOut = ();
set VSM = ();
set pedir = 2
set VSMmul = ()
set Brute1Max = 4;
set Brute1Delta = 4;
set DoBrute2 = 1;
set SubSamp1 = 100;
set DoAbs = 0;
set SaveCurReg = 1;
set featdir = ();
set InitCostFile = ();
set CostFailThresh = ();
set fwhm = ();
set CoRegRefMask = 1;

set LHOnly = 0;
set RHOnly = 0;
set Slope1 = 0.5;
set Slope2 = 0.5;
set Offset2 = 0;
set Contrast = ()
set TolF = 1e-8;
set Tol1D = 1e-3;
set nPowellMax = 36;
set DOF = 6;
set coregDOF = 6;
set FSLDOF = 6;
set FSLSwapTrans = 0;
set FSLBETMov = 0;
set DoPass1 = 1;
set InitSurfCostOnly = 0;

set WMProjAbs = 2;
set GMProjFrac = 0.5;
set GMProjAbs = ();
set IncludeZeroVoxels = 0

set nSubSamp = 1;
set Interp = trilinear
set UseEPIMask = 0;
set UseCortexLabel = 1;

set RMSFile = ();
set RMSFile0 = ();
set templateout = ();
set OutVol = ();
set fsvol  = brainmask;
set initsurfcost = ();
set surfcost = ();
set surfcon  = ();
set surfname = "white";
set MaskLabel = ();
set lhmask = ();
set rhmask = ();
set TargetVol = (); # Full path, instead of SUBJETS_DIR/mri/orig.mgz

set RandInitMax = ();

set debug = 0;
set tmpdir = ();
set cleanup = 1;
set PrintHelp = 0;
set nolog = 0;

set DoInitBest = 0;
set InitBestList = ()

if($?FS_OMP_NUM_THREADS) then
  setenv OMP_NUM_THREADS $FS_OMP_NUM_THREADS 
else
  setenv OMP_NUM_THREADS 1
endif

if($#argv == 0) goto usage_exit;
set n = `echo $argv | egrep -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'`;

if($#outreg)   set outdir = `dirname $outreg`;
if(! $#outreg) set outdir = `dirname $initsurfcost`;

if($#tmpdir == 0) set tmpdir = $outdir/tmp.bbregister.$$
mkdir -p $tmpdir
echo tmp $tmpdir

if(! $nolog) then
  set outregdir  = `dirname $outreg`
  set outregbase = `basename $outreg .dat`
  set outregbase = `basename $outregbase .lta`
  set LF = $outregdir/$outregbase.log
  if(-e $LF) mv $LF $LF.old
else
  set LF = /dev/null
endif

echo "Log file is $LF"

echo "Logfile for bbregister" >> $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
echo $VERSION |& tee -a $LF
uname -a |& tee -a $LF
echo "FREESURFER_HOME $FREESURFER_HOME" |& tee -a $LF

# Loop through init methods
if($DoInitBest) then
  echo "Finding best from init methods $InitBestList" | tee -a $LF
  rm -f $tmpdir/mincost.dat
  foreach init ($InitBestList)
    echo "#I# $init `date`" | tee -a $LF
    set outregbase = `basename $outreg`;
    set regB = $tmpdir/init.$init/$outregbase
    set cmd = (bbregister $inputargs --no-init-best --init-$init --reg $regB --nolog \
      --tmpdir $tmpdir/init.$init/tmp.bbregister)
    echo "$cmd" | tee -a $LF
    $cmd | tee -a $LF
    if($status) exit 1;
    set mincost = `cat $regB.mincost`
    echo "$init $mincost[1]" >> $tmpdir/mincost.dat
  end
  cat $tmpdir/mincost.dat | sort -n -k 2 > $tmpdir/mincost.sort.dat
  set winner = `head -n 1 $tmpdir/mincost.sort.dat | awk '{print $1}'`
  echo "MM mincost $winner -------------------------" | tee -a $LF
  cat $tmpdir/mincost.sort.dat| tee -a $LF
  echo "MM mincost -------------------------" | tee -a $LF
  set cmd = (cp $tmpdir/init.$winner/$outregbase"*" `dirname $outreg`)
  echo $cmd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit 1;
  goto final;
  set tSecEnd = `date '+%s'`;
  @ tSecRun = $tSecEnd - $tSecStart;
  echo "" | tee -a $LF
  echo " " |& tee -a $LF
  echo "Started at $StartTime " |& tee -a $LF
  echo "Ended   at `date`" |& tee -a $LF
  echo "BBR-Run-Time-Sec $tSecRun" |& tee -a $LF
  echo " " |& tee -a $LF
  echo "bbregister Done" |& tee -a $LF

  echo "To check results, run:" |& tee -a $LF
  if($#templateout) then
    echo "tkregister2 --mov $templateout --reg $outreg --surf" |& tee -a $LF
  else
    echo "tkregister2 --mov $movvol --reg $outreg --surf" |& tee -a $LF
  endif
  echo " "
  exit 0;
endif

# Create template
if($#templateout) then
  set template = $templateout
else
  set template = $tmpdir/template.nii
endif
set cmd = (mri_convert $movvol $template)
if($#frame != 0) set cmd = ($cmd --frame $frame)
if($midframe) set cmd = ($cmd --mid-frame)
echo $cmd | tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;

if($#fwhm) then
  # This can be done for testing purposes
  echo "Smoothing template by $fwhm" | tee -a $LF
  set cmd = (mri_fwhm --so --fwhm $fwhm --i $template  --o $template)
  echo $cmd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit 1;
endif

if($#InitReg == 0) then
  set regvol = $template
  set InitReg = $tmpdir/reg.init.dat
  if($#intvol) then
    set regvol = $intvol
    set InitReg = $tmpdir/reg.intermediate.dat
  endif
  if($InitFSL) then
    set cmd = (fslregister --s $subject --mov $regvol --reg $InitReg \
      --niters 1 --maxangle 90 --tmp $tmpdir/fslregister \
      --dof $FSLDOF --fsvol $fsvol.mgz)
    if($FSLBETMov)   set cmd = ($cmd --betmov)
    if(! $FSLBETMov) set cmd = ($cmd --nobetmov)
    if($FSLSwapTrans) set cmd = ($cmd --allow-swap --trans)
    if($nolog) set cmd = ($cmd --nolog)
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
  endif
  if($InitSPM) then
    set cmd = (spmregister --s $subject --mov $regvol --reg $InitReg \
      --tmp $tmpdir/spmregister)
    if($nolog) set cmd = ($cmd --nolog)
    if($SPMUseNII) set cmd = ($cmd --nii)
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
  endif
  if($InitCoreg) then
    set cmd = (mri_coreg --s $subject --mov $regvol --regdat $InitReg \
      --reg $tmpdir/mri_coreg.lta --nthreads $OMP_NUM_THREADS --dof $coregDOF \
      --sep 4 --ftol .0001 --linmintol .01)
    set apas = $SUBJECTS_DIR/$subject/mri/aparc+aseg.mgz
    if($CoRegRefMask == 0 || ! -e $apas) set cmd = ($cmd --no-ref-mask)
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
  endif
  if($InitRR) then
    set rrlta = $tmpdir/rr.lta
    set cmd = (mri_robust_register --mov $regvol \
     --dst $SUBJECTS_DIR/$subject/mri/$fsvol.mgz \
     --lta $rrlta --cost NMI --nosym);
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
    set cmd = (tkregister2_cmdl --mov $regvol --s $subject --lta $rrlta --reg $InitReg)
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
  endif
  if($InitHeader) then
    set cmd = (tkregister2_cmdl --s $subject --mov $regvol \
      --regheader --reg $InitReg --noedit)
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
  endif
  if($#intvol) then
    # Intermediate volume
    set InitReg0 = $tmpdir/reg.init.dat
    set cmd = (tkregister2_cmdl --s $subject --mov $template \
      --int $intvol $InitReg  --noedit --reg $InitReg0)
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
    set InitReg = $InitReg0;
  endif
else
  IsLTA --r $InitReg --o $tmpdir/islta.$$ | tee -a $LF
  if($status) exit 1;
  set InitRegIsLTA = `cat $tmpdir/islta.$$`
  if($InitRegIsLTA) then
    echo "Changing from LTA to register.dat" | tee -a $LF
    set InitRegTKR = $tmpdir/init.reg.$$.dat
    set cmd = (lta_convert --inlta $InitReg --outreg $InitRegTKR)
    echo $cmd | tee -a $LF
    $cmd | tee -a $LF
    if($status) exit 1;
    set InitReg = $InitRegTKR;
 endif
endif

if($#InitRegOut) cp $InitReg $InitRegOut

# Multiply the VSM by the given value
if($#VSM && $#VSMmul) then
  set vsmmulvol = $tmpdir/vsm.mul.nii.gz
  set ud = `UpdateNeeded $vsmmulvol $VSM`
  if($ud) then
    set cmd = (mri_concat --mul $VSMmul $VSM --o $vsmmulvol)
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
  endif
  set VSM = $vsmmulvol
endif

# Pass 1
set Pass1Reg = $tmpdir/bbr.pass1.dat
if($DoPass1) then
  set cmd = (mri_segreg --mov $template --init-reg $InitReg --out-reg $Pass1Reg\
    --subsamp-brute $SubSamp1 --subsamp $SubSamp1 --tol 1e-4 --tol1d 1e-3 \
    --brute -$Brute1Max $Brute1Max $Brute1Delta --surf $surfname);
  if($#GMProjFrac) set cmd = ($cmd --gm-proj-frac $GMProjFrac)
  if($#GMProjAbs)  set cmd = ($cmd --gm-proj-abs  $GMProjAbs)
  if("$Contrast" == "-1") set cmd = ($cmd --wm-gt-gm $Slope1);
  if("$Contrast" == "+1") set cmd = ($cmd --gm-gt-wm $Slope1);
  if($#MaskLabel)  set cmd = ($cmd --label $MaskLabel);
  if($#lhmask)  set cmd = ($cmd --lh-mask $lhmask);
  if($#rhmask)  set cmd = ($cmd --rh-mask $rhmask);
  if($LHOnly) set cmd = ($cmd --lh-only);
  if($RHOnly) set cmd = ($cmd --rh-only);
  if($#initsurfcost) set cmd = ($cmd --init-surf-cost $initsurfcost);
  if($InitSurfCostOnly) set cmd = ($cmd --init-surf-cost-only);
  if($UseEPIMask) set cmd = ($cmd --mask);
  if(! $UseCortexLabel) set cmd = ($cmd --no-cortex-label);
  if($#VSM) set cmd = ($cmd --vsm $VSM $pedir);
  if($IncludeZeroVoxels) set cmd = ($cmd --include-zero-voxels)
  if($#TargetVol) set cmd = ($cmd --target-volume $TargetVol)
  if($#RandInitMax) then
    # Only for testing
    set cmd = ($cmd --trans-rand $RandInitMax --rot-rand $RandInitMax);
  endif
  echo $cmd | tee -a $LF
  $cmd |& tee -a $LF
  if($status) exit 1;
  if($InitSurfCostOnly) then
    if($cleanup) rm -r $tmpdir
    echo "InitSurfCostOnly requested, so exiting now" |& tee -a $LF
    echo "bbregister Done" |& tee -a $LF
    exit 0;
  endif
else
  cp $InitReg $Pass1Reg
endif

# Pass 2
set MinCostFile = $outreg.mincost; # This will be the final cost
set ParamFile = $outreg.param; # This will be the final parameters
set cmd = (mri_segreg --mov $template --init-reg $Pass1Reg \
  --out-reg $outreg --interp $Interp \
  --wm-proj-abs $WMProjAbs --tol $TolF --tol1d $Tol1D\
  --c0 $Offset2 --mincost $MinCostFile --dof $DOF \
  --nmax $nPowellMax --param $ParamFile --surf $surfname);
if($DoBrute2) set cmd = ($cmd --brute -0.1 0.1 0.1)
if($SaveCurReg) set cmd = ($cmd --cur-reg $tmpdir/reg.curopt.dat)
if($#GMProjFrac) set cmd = ($cmd --gm-proj-frac $GMProjFrac)
if($#GMProjAbs)  set cmd = ($cmd --gm-proj-abs  $GMProjAbs)
if($#nSubSamp) set cmd = ($cmd --nsub $nSubSamp);
if($#lhmask)  set cmd = ($cmd --lh-mask $lhmask);
if($#rhmask)  set cmd = ($cmd --rh-mask $rhmask);
if($LHOnly) set cmd = ($cmd --lh-only);
if($RHOnly) set cmd = ($cmd --rh-only);
if($UseEPIMask) set cmd = ($cmd --mask);
if(! $UseCortexLabel) set cmd = ($cmd --no-cortex-label);
if($#VSM) set cmd = ($cmd --vsm $VSM $pedir);
if($#OutVol) set cmd = ($cmd --o $OutVol);
if($#MaskLabel)  set cmd = ($cmd --label $MaskLabel);
if("$Contrast" == "-1") set cmd = ($cmd --wm-gt-gm $Slope2);
if("$Contrast" == "+1") set cmd = ($cmd --gm-gt-wm $Slope2);
if($#initsurfcost && ! $DoPass1)  set cmd = ($cmd --init-surf-cost $initsurfcost);
if($InitSurfCostOnly) set cmd = ($cmd --init-surf-cost-only);
if($#surfcost)  set cmd = ($cmd --surf-cost $surfcost);
if($#surfcon)   set cmd = ($cmd --surf-con  $surfcon);
if($DoAbs) set cmd = ($cmd --penalty-abs)
if($#InitCostFile) set cmd = ($cmd --initcost $InitCostFile);
if($#RMSFile0) set cmd = ($cmd --rms $RMSFile0);
if($IncludeZeroVoxels) set cmd = ($cmd --include-zero-voxels)
if($#TargetVol) set cmd = ($cmd --target-volume $TargetVol)
echo $cmd | tee -a $LF
$cmd |& tee -a $LF
if($status) exit 1;
if($InitSurfCostOnly) then
  if($cleanup) rm -r $tmpdir
  echo "InitSurfCostOnly requested, so exiting now" |& tee -a $LF
  echo "bbregister Done" |& tee -a $LF
  exit 0;
endif

echo "MinCost: `cat $MinCostFile`" | tee -a $LF

# Compute RMS wrt the initial reg (however it was created)
if($#RMSFile) then
  set rmsdat = ();
  if(! $RHOnly) then
    set lhrms = $tmpdir/rh.rms
    set cmd = (mri_surf2surf --reg $InitReg --reg-diff $outreg \
     --sval-xyz white --rms $lhrms --s $subject --hemi lh)
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
    set rmsdat = (`cat $lhrms`);
  endif
  if(! $LHOnly) then
    set rhrms = $tmpdir/rh.rms
    set cmd = (mri_surf2surf --reg $InitReg --reg-diff $outreg \
     --sval-xyz white --rms $rhrms --s $subject --hemi rh)
    echo $cmd | tee -a $LF
    $cmd |& tee -a $LF
    if($status) exit 1;
    set rmsdat = ($rmsdat `cat $rhrms`);
  endif
  echo $rmsdat > $RMSFile
  echo Final RMS $rmsdat | tee -a $LF
endif

final:

if($#outfsl) then
  set cmd = (tkregister2_cmdl --mov $movvol --reg $outreg \
    --noedit --fslregout $outfsl);
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) exit 1;
endif

if($#outlta) then
  set cmd = (tkregister2_cmdl --mov $movvol --reg $outreg \
    --noedit --ltaout $outlta);
  if($#TargetVol) set cmd = ($cmd --targ $TargetVol)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) exit 1;
endif

if($#featdir) then
  cp $outreg $featdir/reg/freesurfer/register.dat
  set reg152 = $featdir/reg/freesurfer/anat2std.register.dat
  set cmd = (mni152reg --s $subject --o $reg152)
  echo $cmd | tee -a $LF
  $cmd | tee -a $LF
  if($status) exit 1;
endif

# Cleanup
if($cleanup) then
  echo "Cleaning up" |& tee -a $LF
  rm -r $tmpdir
endif

set tSecEnd = `date '+%s'`;
@ tSecRun = $tSecEnd - $tSecStart;

echo " " |& tee -a $LF
echo "Started at $StartTime " |& tee -a $LF
echo "Ended   at `date`" |& tee -a $LF
echo "BBR-Run-Time-Sec $tSecRun" |& tee -a $LF
echo " " |& tee -a $LF
echo "bbregister Done" |& tee -a $LF

echo "To check results, run:" |& tee -a $LF
set tkrreg = $outreg
if($#outlta) set tkrreg = $outlta
set hemi = ""
if($LHOnly) set hemi = "--lh-only"
if($RHOnly) set hemi = "--rh-only"
if($#templateout) then
  echo "tkregisterfv --mov $templateout --reg $tkrreg --surfs $hemi --sd $SUBJECTS_DIR" |& tee -a $LF
else
  echo "tkregisterfv --mov $movvol --reg $tkrreg --surfs $hemi --sd $SUBJECTS_DIR" |& tee -a $LF
endif
echo " "

if($#CostFailThresh) then
  set FinalCost = `cat  $MinCostFile | awk '{print $1}'`
  set Fail = `echo "$FinalCost > $CostFailThresh" | bc -l`
  if($Fail) then
    echo "ERROR: bbregister ran to completion without error, but the " |& tee -a $LF
    echo "  final cost ($FinalCost) exceed the given threshold ($CostFailThresh)." |& tee -a $LF
    echo "  Manually check the registration. Look for gross misregistrations,"|& tee -a $LF
    echo "  left-right flips, or different subject identities for mov and targ."|& tee -a $LF
    echo "  A gross misregistration probably means the initialization failed."|& tee -a $LF
    echo "  If the registration is ok, rerun with a higher or no threshold."|& tee -a $LF
    echo " "
    exit 1;
  endif
endif

exit 0;
###############################################

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

  set flag = $argv[1]; shift;

  switch($flag)

    case "-h"
    case "-u"
    case "-usage"
    case "--usage"
    case "-help"
    case "--help"
      set PrintHelp = 1;
      goto usage_exit;
      breaksw

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

    case "--feat":
      if ( $#argv < 1) goto arg1err;
      set featdir = $argv[1]; shift;
      if(! -e $featdir) then
        echo "ERROR: cannot find $featdir"
        exit 1;
      endif
      set movvol = `stem2fname $featdir/example_func`
      if($status) then
        echo "$movvol"
        exit 1;
      endif
      mkdir -p $featdir/reg/freesurfer
      if($status) exit 1;
      set outreg = $featdir/reg/freesurfer/anat2exf.register.dat
      set InitFSL = 1;
      set InitSPM = 0;
      set InitHeader = 0;
      set InitReg = ();
      set Contrast = +1;
      breaksw

    case "--s-from-reg":
      if($#argv < 1) goto arg1err;
      set tmp = $argv[1]; shift;
      if(! -e $tmp) then
        echo "ERROR: cannot find $tmp"
        exit 1;
      endif
      #set subject = `head -n 1 $tmp`
      set subject = `reg2subject --r $tmp`;
      breaksw

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

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

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

    case "--vsm":
      if ( $#argv < 1) goto arg1err;
      set VSM = $argv[1]; shift;
      if(! -e $VSM) then
        echo "ERROR: cannot find $VSM"
        exit 1;
      endif
      set UseEPIMask = 1; # Turn on EPI B0 mask by default
      breaksw

    case "--vsm-pedir":
      if ( $#argv < 1) goto arg1err;
      set pedir = $argv[1]; shift;
      if($pedir != 1 && $pedir != 2 && $pedir != 3 && \
         $pedir != +1 && $pedir != +2 && $pedir != +3 && \
         $pedir != -1 && $pedir != -2 && $pedir != -3 ) then
        echo "ERROR: pedir = $pedir, must be +/-1, +/-2, or +/-3"
        exit 1;
      endif
      breaksw

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

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

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

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

    case "--fsl-swap-trans":
      set FSLSwapTrans = 1;
      breaksw

    case "--init-coreg":
      set InitFSL = 0;
      set InitSPM = 0;
      set InitCoreg = 1;
      set InitRR = 0;
      set InitHeader = 0;
      set InitReg = ();
      breaksw
 
    case "--init-rr":
      set InitFSL = 0;
      set InitSPM = 0;
      set InitCoreg = 0;
      set InitRR = 1;
      set InitHeader = 0;
      set InitReg = ();
      breaksw
 
    case "--init-fsl":
      set InitFSL = 1;
      set InitSPM = 0;
      set InitCoreg = 0;
      set InitRR = 0;
      set InitHeader = 0;
      set InitReg = ();
      breaksw
 
    case "--init-spm":
      set InitFSL = 0;
      set InitSPM = 1;
      set InitCoreg = 0;
      set InitRR = 0;
      set InitHeader = 0;
      set InitReg = ();
      breaksw
 
    case "--spm-nii":
      set SPMUseNII = 1;
      breaksw
 
    case "--regheader":
    case "--reg-header":
    case "--init-header":
      set InitFSL = 0;
      set InitSPM = 0;
      set InitCoreg = 0;
      set InitHeader = 1;
      set InitReg = ();
      breaksw
 
    case "--init-reg":
      if ( $#argv < 1) goto arg1err;
      set InitReg = $argv[1]; shift;
      if(! -e $InitReg) then
        echo "ERROR: cannot find $InitReg"
        exit 1;
      endif
      set subject = `reg2subject --r $InitReg`
      set InitFSL = 0;
      set InitSPM = 0;
      set InitCoreg = 0;
      set InitHeader = 0;
      breaksw

    case "--no-coreg-ref-mask":
      set CoRegRefMask = 0;
      breaksw
    case "--coreg-ref-mask":
      set CoRegRefMask = 1;
      breaksw

    case "--init-reg-out":
      if( $#argv < 1) goto arg1err;
      set InitRegOut = $argv[1]; shift;
      breaksw

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

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

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

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

    case "--abs":
      # Experimental
      set Contrast = +1;
      set DoAbs = 1;
      breaksw

    case "--bold":
    case "--dti":
    case "--T2":
    case "--t2":
      set Contrast = +1;
      breaksw

    case "--T1":
    case "--t1":
      set Contrast = -1;
      breaksw

    case "--s-from-reg":
      if ( $#argv < 1) goto arg1err;
      set tmpfile = $argv[1]; shift;
      set subject = `head -n 1 $tmpfile`;
      if($status) then
        echo "$subject"
        exit 1;
      endif
      breaksw

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

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

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

    case "--mid-frame":
      set midframe = 1;
      breaksw

    case "--target-volume"
      # Full path, instead of SUBJETS_DIR/mri/orig.mgz
      set TargetVol = $argv[1]; shift;
      breaksw

    case "--lh-only":
      set LHOnly = 1;
      set RHOnly = 0;
      breaksw

    case "--rh-only":
      set LHOnly = 0;
      set RHOnly = 1;
      breaksw

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

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

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

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

    case "--no-brute2":
      set DoBrute2 = 0;
      breaksw

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

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

    case "--init-surf-cost":
      if ( $#argv < 1) goto arg1err;
      set initsurfcost = $argv[1]; shift;
      breaksw

    case "--init-surf-cost-only":
      set InitSurfCostOnly = 1;
      set nolog = 1;
      breaksw

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

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

    case "--include-zero-voxels":
      set IncludeZeroVoxels = 1
      breaksw
    case "--no-include-zero-voxels":
      set IncludeZeroVoxels = 1
      breaksw

    case "--tol":
      if ( $#argv < 1) goto arg1err;
      set Tol = $argv[1]; shift;
      set TolF = $Tol
      # set Tol1D = $Tol; only set tolf!
      breaksw

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

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

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

    case "--gm-proj-frac":
      if($#argv < 1) goto arg1err;
      set GMProjFrac = $argv[1]; shift;
      set GMProjAbs = ();
      breaksw

    case "--gm-proj-abs":
      if($#argv < 1) goto arg1err;
      set GMProjAbs = $argv[1]; shift;
      set GMProjFrac = ();
      breaksw

    case "--wm-proj-abs":
      if($#argv < 1) goto arg1err;
      set WMProjAbs = $argv[1]; shift;
      breaksw

    case "--proj-abs":
      if($#argv < 1) goto arg1err;
      set ProjAbs = $argv[1]; shift;
      set WMProjAbs = $ProjAbs;
      set GMProjAbs = $ProjAbs;
      set GMProjFrac = ();
      breaksw

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

    case "--6":
      # Default anyway
      set DOF = 6;
      set coregDOF = 6;
      breaksw

    case "--9":
      set DOF = 9;
      set coregDOF = 9;
      breaksw

    case "--12":
      set DOF = 12;
      set coregDOF = 12;
      breaksw

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

    case "--no-pass1":
      set DoPass1 = 0;
      breaksw

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

    case "--fsl-bet-mov":
      set FSLBETMov = 1;
      breaksw
    case "--no-fsl-bet-mov":
      set FSLBETMov = 0;
      breaksw

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

    case "--rms0":
      # only works with --no-pass1
      if($#argv < 1) goto arg1err;
      set RMSFile0 = $argv[1]; shift;
      breaksw

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

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

    case "--nearest":
      set Interp = nearest
      breaksw

    case "--trilin":
    case "--trilinear":
      set Interp = trilinear
      breaksw

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

    case "--mask":
      if($#argv < 2) goto arg2err;
      set lhmask = $argv[1]; shift;
      set rhmask = $argv[1]; shift;
      if(! -e $lhmask) then
        echo "ERROR: cannot find $lhmask"
        exit 1;
      endif
      if(! -e $rhmask) then
        echo "ERROR: cannot find $rhmask"
        exit 1;
      endif
      breaksw

    case "--epi-mask":
      set UseEPIMask = 1;
      breaksw
    case "--no-epi-mask":
      set UseEPIMask = 0;
      breaksw

    case "--no-cortex-label":
      set UseCortexLabel = 0;
      breaksw

    case "--no-save-cur-reg":
      set SaveCurReg = 0;
      breaksw

    case "--init-best":
      set DoInitBest = 1;
      set InitFSL = 0;
      set InitSPM = 0;
      set InitCoreg = 0;
      set InitRR = 0;
      set InitHeader = 0;
      set InitReg = ();
      breaksw
    case "--no-init-best":
      set DoInitBest = 0;
      breaksw
    case "--init-best-fsl":
      set InitBestList = ($InitBestList fsl); 
      #set DoInitBest = 1; Do not do this!
      breaksw
    case "--init-best-spm":
      set InitBestList = ($InitBestList spm); 
      breaksw
    case "--init-best-header":
      set InitBestList = ($InitBestList header); 
      breaksw
    case "--init-best-rr":
      set InitBestList = ($InitBestList rr); 
      breaksw
    case "--init-best-coreg":
      set InitBestList = ($InitBestList coreg); 
      breaksw

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

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

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

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

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

    case "--nolog":
      set nolog = 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 && $#TargetVol == 0) then
    # Adding this stuff about the TargetVol is a bit of a hack
    echo "ERROR: cannot find $subject in $SUBJECTS_DIR"
    exit 1;
  endif

  if($#movvol == 0) then
    echo "ERROR: must spec an movput vol"
    exit 1;
  endif

  if($#Contrast == 0) then
    echo "ERROR: you must specify a contrast."
    echo " use --bold, --dti, --t2, or --t1,"
    echo " which ever is most appropriate"
    exit 1;
  endif

  if($InitSurfCostOnly && $#initsurfcost == 0) then
    echo "ERROR: must supply --init-surf-cost with --init-surf-cost-only"
    exit 1;
  endif

  # Check whether output registration file is really lta format
  if($#outreg != 0 && $#outlta == 0) then
    set tmp1 = `basename $outreg `
    set tmp2 = `basename $outreg .lta`
    if("$tmp1" != "$tmp2") then
      set outlta = $outreg
      set outreg = ()
    endif
  endif

  if($#outreg == 0 && ! $InitSurfCostOnly) then
    if ($#outlta == 0) then
      echo "ERROR: must spec an output reg file (using either --reg or --lta)"
      exit 1;
    else
      set outreg = ${outlta:r}.dat
    endif
  endif

  if($#frame && $midframe) then
    echo "ERROR: cannot --frame AND --mid-frame"
    exit 1;
  endif

  if(! $DoInitBest) then
    if(! $InitFSL && ! $InitSPM && ! $InitCoreg && ! $InitRR && ! $InitHeader && ! $#InitReg) then
      echo "ERROR: must spec an init method"
      exit 1;
    endif
  else
    if($InitFSL || $InitSPM ||  $InitCoreg || $InitRR || $InitHeader || $#InitReg) then
      echo "ERROR: cannot spec an init method with --init-best"
      exit 1;
    endif
    if($#InitBestList == 0) set InitBestList = (fsl spm coreg rr header)
    if($#InitBestList == 1) then
      echo "ERROR: must use more than one init method with --init-best"
      exit 1;
    endif
  endif

  if($?BBR_TEST_TOLERANCE) then
    # This allows a program at a much higher level to set the tolerances.
    # This is useful in the case where only the execution needs to be
    # tested (and the output itself is not important). Setting the tolerances
    # to be high makes BBR run faster.
    set TolF  = $BBR_TEST_TOLERANCE
    set Tol1D = $BBR_TEST_TOLERANCE
    echo ""
    echo "WARNING: using BBR_TEST_TOLERANCE $BBR_TEST_TOLERANCE"
    echo ""
    sleep 2
  endif

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:

if($PrintHelp) then
  cat $0 | \
    awk 'BEGIN{prt=0}{if(prt) print $0; if($1 == "BEGINHELP") prt = 1 }'
else
  echo "Usage: bbregister --s <subj> --mov <volid> --reg <regfile> --<contrast>"
  echo "Help:  bbregister --help"
endif

exit 1;

#---- Everything below is printed out as part of help -----#
#-- During make install, the output of 'fsPrintHelp bbregister.help.xml' -- #
#-- is concatenated to this file --#
BEGINHELP
				Help

NAME
	bbregister

SYNOPSIS
	bbregister --s <subj> --mov <volid> --reg <regfile> \
		--<contrast>

DESCRIPTION
	This program performs within-subject, cross-modal registration using a
	boundary-based cost function. The registration is constrained to be 6 
	DOF (rigid). It is required that you have an anatomical scan of the 
	subject that has been analyzed in freesurfer.

POSITIONAL ARGUMENTS
	None

REQUIRED FLAGGED ARGUMENTS
	--s subject
		FreeSurfer subject name as found in $SUBJECTS_DIR.

	--mov volid
		"Moveable" volume. This is the template for the cross-modal 
		volume. Eg, for fMRI, it is the volume used for motion 
		correction.

	--reg register.{dat,lta}
		Output FreeSurfer registration file (tkregister-style or LTA 
		format).

	--vsm VoxelShiftMap 
		Include B0 distortion correction while performing 
		registration. A voxel shift map can be created from a B0 map 
		with the epidewarp.fsl script (in FreeSurfer). Automatically 
		turns on epi B0 mask; so run with --no-epi-mask if you do not 
		want that.

	--vsm-pedir PEdirection 
		Phase encode direction. +/-1=+/-x, +/-2=+/-y, +/-3=+/-z 
		(default is +2)

	--vsm-mul MultFactor 
		Multiply VSM by the given factor

	Initialization Arguments

	--init-coreg (Default)
		Initialize the FreeSurfer mri_coreg program.

	--no-coreg-ref-mask
		Do NOT use aparc+aseg.mgz as a reference mask. It will not 
		attempt to use aparc+aseg if it is not there. 

	--init-header
		Assume that the geometry information in the cross-modal and 
		anatomical are sufficient to get a close voxel-to-voxel 
		registration. This usually is only the case if they were 
		acquired in the same session.

	--init-reg initregfile
		Supply an initial registration matrix; can be LTA.

	Contrast Arguments (one required)

	--t1
		Assume t1 contrast, ie, White Matter brighter than Grey Matter

	--t2
		Assume t2 contrast, ie, Gray Matter brighter than White Matter

	--bold
		Same as --t2

	--dti
		Same as --t2

OPTIONAL FLAGGED ARGUMENTS
	--int intvol
		Supply a volume to use an an intermediate volume when 
		performing registration. This is useful for when the 
		cross-modal is volume is a partial field-of-view (FoV). If you
		acquire in the same session a whole-head FoV, then pass the 
		whole-head as the intermediate and the partial as the 
		moveable.

	--mid-frame
		reg to middle frame (not with --frame)

	--frame frameno
		reg to frameno (default 0=1st)

	--template-out template
		save template (good with --frame)

	--o outvol
		resample mov and save as outvol

	--s-from-reg reg
		get subject name from regfile

	--rms rmsfile
		RMS change in cortical surface position

	--fslmat flirt.mtx
		output an FSL FLIRT matrix

	--lta output.lta
		output an LTA registration matrix (This flag can be used along
		with or instead of --reg!) 

	--lh-only
		only use left hemi

	--rh-only
		only use right hemi

	--slope1 slope1
		cost slope for 1st stage (default is 0.5)

	--slope2  slope2
		cost slope for 2nd stage (default is 0.5)

	--offset2 offset2
		cost offset for 2nd stage (default is 0)

	--tol1d tol1d
		2nd stage 1D tolerance 

	--tol tol
		2nd stage loop tolerance (same as --tolf)

	--tolf  tolf
		2nd stage loop tolerance (same as --tol)

	Be careful making these more stringent as they can cause underflows 
	and NaNs.

	--nmax nPowellMax
		set max number of iterations (default 36)

	--rand-init randmax
		randomly change input to 1st stage reg

	--gm-proj-frac frac
		default is 0.5

	--gm-proj-abs abs
		use absolute instead of relative

	--wm-proj-abs  dist
		2nd stage, default is 2mm

	--proj-abs  dist
		use wm and gm proj abs in 2nd stage

	--subsamp nsub
		2nd stage vertex subsampling, default is 1

	--nearest
		2nd stage, use nearest neighbor interp (defalt is trilinear)

	--epi-mask
		mask out brain edge and B0 regions (1st and 2nd stages)

	--no-epi-mask
		Turn off epi B0 mask

	--no-cortex-label
		Do not use ?h.cortex.label to mask. 

	--label labelfile
		Use label to mask. Only vertices in the label are used.

	--brute1max max
		pass 1 search -max to +max (default 4)

	--brute1delta delta
		pass 1 search -max to +max step delta (default 4)

	--subsamp1 nsubsamp
		pass 1 vertex subsampling (default 1)

	--no-pass1
		turn off pass 1

	--no-brute2
		turn off brute force search on pass 2

	--include-zero-voxels
		Included zero-valued voxels in cost, including those outside 
		of the FoV. 

	--surf surfname
		change surface to surfname from ?h.white 

	--surf-cost basename
		saves final cost as basename.?h.mgh

	--surf-con  basename
		saves final contrast as basename.?h.mgh

	--init-surf-cost basename
		saves initial cost as basename.?h.mgh

	--init-reg-out outinitreg
		save initial reg

	--initcost initcostfile
		save initial cost

	--cost-fail costvalue
		exit with error if the final cost exceeds costvalue 

	--12, --9, --6
		Use 12 or 9 or 6 DOF (6 is default)

	--coreg-dof coregDOF
		Set dof to use for mri_coreg

	--spm-nii
		Use NIFTI format as input to SPM when using --init-spm 
		(spmregister). Ordinarily, it uses ANALYZE images to be 
		compatible with older versions of SPM, but this has caused 
		some left-right reversals in SPM8.

	--feat featdir
		FSL FEAT directory. Sets mov to featdir/example_func, uses 
		--init-fsl, --bold, sets reg to featdir/reg/freesurfer/
		anat2exf.register.dat. This replaces reg-feat2anat.

	--fsl-bet-mov
		When using FSL, apply BET to mov volume (turned by default).

	--init-rr
		Initialize using NMI as implemented in mri_robust_register.

	--init-fsl
		Initialize using FSL FLIRT (requires that FSL be installed).

	--init-spm
		Initialize using SPM spm_coreg (requires that SPM and matlab 
		be installed).

	--tmp tmpdir
		temporary dir (implies --nocleanup)

	--nocleanup
		do not delete temporary files

	--version
		print version and exit

	--help
		print help and exit

EXAMPLE 1
	Intialize with FLIRT, view result with tkregister2:
	bbregister --s bert --mov func.nii --init-fsl --reg 
	register.dat
	tkregister2 --mov func.nii --reg register.dat --surf 

EXAMPLE 2
	Intialize with SPM, view result with tkregister2:
	bbregister --s bert --mov func.nii --init-spm --reg 
	register.dat
	tkregister2 --mov func.nii --reg register.dat --surf

EXAMPLE 3
	Register a partial FoV using whole FoV as intermdediate:
	bbregister --s bert --mov partial.nii --init-fsl --int whole.nii --reg
	register.partial.dat
	tkregister2 --mov partial.nii --reg register.partial.dat --surf

REPORTING
	Report bugs to <freesurfer@nmr.mgh.harvard.edu>