#!/usr/bin/env tcsh

@global_parse `basename $0` "$*" ; if ($status) exit 0

# Why use f00*, f01*, f02*, etc.? is iidx instead?

# written by PA Taylor (NIMH, NIH, USA)
# started March, 2017

# --------------------- revision history -------------------------
#
#set version   = "1.0";    set rev_dat   = "March 10, 2017"
#   + start: align t1 to t2
#
#set version   = "1.1";    set rev_dat   = "March 11, 2017"
#   + change name...
#   + help text added
#   + no-clean option added
#
#set version   = "1.2";    set rev_dat   = "Apr 12, 2017"
#   + standardize I/O, helpfile, option names
#
#set version   = "1.3";    set rev_dat   = "Apr 26, 2017"
#   + newer wdir convention
#
#set version   = "1.4";    set rev_dat   = "May 1, 2017"
#   + Align t2w -> t1w and reverse (better masking/weighting)
#
#set version   = "1.5";    set rev_dat   = "June 20, 2017"
#   + new opt: can have mask for t2w vol
#
#set version   = "1.6";    set rev_dat   = "July 25, 2017"
#   + new opt: can have more than warp = shift+rot only; default is
#   just that, still
#
#set version   = "1.7";    set rev_dat   = "Aug 11, 2017"
#   + switch a ">>" to ">" for general system compatibility
#
#set version   = "1.8";    set rev_dat   = "Sep 04, 2017"
#   + update to work with newer @chauffeur*, -prefix only
#
#set version   = "1.9";    set rev_dat   = "Jan 12, 2018"
#   + wt vol into working dir
#
#set version   = "1.93";    set rev_dat   = "Feb 8, 2018"
#   + new opt: can skullstrip intermediate t1w to help alignment, via
#   -do_ss_tmp_t1w ... and that file stays in wdir
#
#set version   = "1.0";    set rev_dat   = "Feb 19, 2018"
#   + new filename output structure, no postfix, make QC dir
#
#set version   = "1.1";  set rev_dat   = "July 25, 2018"
#   + match func_range_nz to updated @chauffeur
#
#set version   = "1.2";   set rev_dat   = "Feb 12, 2019"
#     + [PT] change "checks" to use '3dinfo -prefix ...' as a better
#            methodology
#
#set version   = "1.21";   set rev_dat   = "Oct 28, 2020"
#     + [PT] extra QC imaging: initial alignment
#
#set version   = "1.3";   set rev_dat   = "Aug 29, 2021"
#     + [PT] remove -no_fs_prep opt and "DO_PREP_FS" behavior,
#       which was neeever even needed, embarrassingly enough
#
#set version   = "1.31";   set rev_dat   = "Sep 27, 2021"
#     + [PT] chauffeur label_size 3 -> 4, bc imseq.c shifted all sizes
#       down one level
#
#set version   = "1.4";   set rev_dat   = "Feb 23, 2026"
#     + [PT] add new opt so output transform+dset is T2w->T1w
#            add new opt to turn off -cmass
#
set version   = "1.5";   set rev_dat   = "Mar 18, 2026"
#     + [PT] add new opt to open up parameter space for 3dAllineate
#
# ---------------------------------------------------------------

set this_prog = "fat_proc_align_anat_pair"
set tpname    = "${this_prog:gas/fat_proc_//}"
set here      = "$PWD"

# ----------------- find AFNI and set viewer ---------------------

# find AFNI binaries directory and viewer location
set adir      = ""
set my_viewer = ""
which afni >& /dev/null
if ( $status ) then
    echo "** Cannot find 'afni' (???)."
    goto BAD_EXIT
else
    set aa   = `which afni`
    set adir = `dirname $aa`
endif

# default location of viewer: user could modify
set my_viewer = "$adir/@chauffeur_afni"

# ----------------------- set defaults --------------------------

set it1w       = ""
set it2w       = ""
set imtrx      = ""              # user can input 12 DOF made elsewhere
set it2w_mask  = ""

# def: do prep t1w for freesurfing -> matrix has even numbered dims ->
# 1mm iso voxels, or other iso
set DO_SS_T1W  = 0

set RES        = ""              # def: 1mm output res, or smallest vox dim
set OUT_T2W_GRID = 0
set OUT_T1W_GRID = 0             # if doing T2w->T1w

set odir       = ""
set opref      = ""

set DO_VIEWER  = 1               # def: do viewing
set DO_CLEAN   = "1"
set wdir       = "__WORKING_$tpname"
set WDIR_EX  = "1"               # put opref on wdir (unless user names)
set output_cmd = 1               # def: output copy of this command
set cmd_file   = ""              # def: same name as viewer
set qc_prefix  = ""              # def: autoname; user can enter
set warp_opt   = "-warp shift_rotate" # default, but can be changed

set DO_T2W_TO_T1W = 0
set do_use_cmass  = 1

# to control max rot and max translation
set maxrot = ""
set maxshf = ""

# ------------------- process options, a la rr ----------------------

if ( $#argv == 0 ) goto SHOW_HELP

set ac = 1
while ( $ac <= $#argv )
    # terminal options
    if ( ("$argv[$ac]" == "-h" ) || ("$argv[$ac]" == "-help" )) then
        goto SHOW_HELP
    endif
    if ( "$argv[$ac]" == "-ver" ) then
        goto SHOW_VERSION
    endif

    if ( "$argv[$ac]" == "-echo" ) then
        set echo
    

    # --------------- input dset(s) ----------------

    else if ( "$argv[$ac]" == "-in_t1w" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set it1w = "$argv[$ac]"

    else if ( "$argv[$ac]" == "-in_t2w" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set it2w = "$argv[$ac]"

    else if ( "$argv[$ac]" == "-in_t2w_mask" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set it2w_mask = "$argv[$ac]"

    # [PT: July 25, 2017] works like opt of the same name in
    # 3dAllineate; takes shift_only, shift_rotate, etc.
    else if ( "$argv[$ac]" == "-warp" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set warp_opt = "-warp $argv[$ac]"

    # optional input: premade 3dAllineate matrix
    else if ( "$argv[$ac]" == "-matrix" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set imtrx = "$argv[$ac]"

    # ----------------- outputs ---------------------

    else if ( "$argv[$ac]" == "-prefix" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set opref = "$argv[$ac]"

    # ------------------- other opts ---------------

    # invert default alignment dir for output
    else if ( "$argv[$ac]" == "-align_t2w_to_t1w" ) then
        set DO_T2W_TO_T1W = 1

    else if ( "$argv[$ac]" == "-do_use_cmass" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set yn = "$argv[$ac]"
        if ( "$yn" == "Yes" || "$yn" == "1" ) then
            set do_use_cmass = 1
        else if ( "$yn" == "No" || "$yn" == "0" ) then
            set do_use_cmass = 0
        else
            echo "** ERROR, cannot interpret value after '-do_use_cmass':"
            echo "      ${yn}"
            echo "   It must be one of the following: Yes, 1, No, 0"
            goto BAD_EXIT
        endif

    else if ( "$argv[$ac]" == "-workdir" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set wdir = "$argv[$ac]"
        set WDIR_EX  = "0"

    else if ( "$argv[$ac]" == "-maxrot" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set maxrot = "-maxrot $argv[$ac]"

    else if ( "$argv[$ac]" == "-maxshf" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set maxshf = "-maxshf $argv[$ac]"

    else if ( "$argv[$ac]" == "-no_cmd_out" ) then
        set output_cmd = 0

    # set output resolution
    else if ( "$argv[$ac]" == "-newgrid" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set RES = "$argv[$ac]"

    else if ( "$argv[$ac]" == "-do_ss_tmp_t1w" ) then
        set DO_SS_T1W = "1"

    # set output resolution by using t2w grid
    else if ( "$argv[$ac]" == "-out_t2w_grid" ) then
        set OUT_T2W_GRID = 1

    # set output resolution by using t1w grid
    else if ( "$argv[$ac]" == "-out_t1w_grid" ) then
        set OUT_T1W_GRID = 1

    else if ( "$argv[$ac]" == "-no_clean" ) then
        set DO_CLEAN = "0"

    else
        echo "** unexpected option #$ac = '$argv[$ac]'"
        goto BAD_EXIT

    endif
    @ ac += 1
end

# =======================================================================
# ============================ ** SETUP ** ==============================
# =======================================================================

# ============================ input files ==============================

echo "++ Start script version: $version"

# NEED these two inputs
if ( "$it1w" == "" ) then
    echo "** ERROR: no t1w file input???"
    goto BAD_EXIT
else if ( "$it2w" == "" ) then
    echo "** ERROR: no t2w file input???"
    goto BAD_EXIT
endif

# make sure we can read volumes OK
foreach ff ( "$it1w" "$it2w" )
    set check = `3dinfo -prefix "$ff"`
    if ( "$check" == "NO-DSET" ) then
        echo "** ERROR: can't find input file:  $ff"
        goto BAD_EXIT
    else
        echo "++ Found input file:   $ff"
    endif
end

# ===================== output and working dirs =========================

if ( "$opref" == "" ) then
    echo "** ERROR: need '-prefix ...' option provided"
    echo "   See the helpfile for more information."
    goto BAD_EXIT
else
    set odir = `dirname $opref`
    set opref = `basename $opref`
    echo ""
    echo "++ Based on prefix, the output directory will be:"
    echo "     $odir"
    echo "++ Based on prefix, the output prefix will be:"
    echo "     $opref"
    echo ""
endif

# check output directory, use input one if nothing given
if ( ! -e "$odir" ) then
    echo "+* Output directory didn't exist.  Trying to make it now."
    \mkdir "$odir"
endif

set odirqc = "$odir/QC"
if ( "$DO_VIEWER" == "1" ) then
    if ( ! -e "$odirqc" ) then
        \mkdir "$odirqc"
    endif
endif

# and put working directory as subdirectory.
if ( "$WDIR_EX" == "1" ) then
    set wdir = $odir/${wdir}_$opref
else
    set wdir = $odir/$wdir
endif

# make the working directory
if ( ! -e $wdir ) then
    echo "++ Making working directory: $wdir"
    \mkdir $wdir
else
    echo "+* WARNING: Somehow found a premade working directory (?):"
    echo "      $wdir"

    # don't clean preexisting directories-- could be user mistake.
    echo "   NB: will *not* clean it afterwards."
    set DO_CLEAN = "0"
endif

if ( $do_use_cmass ) then
    set cmass = "-cmass"
else
    set cmass = ""
endif

# --------------------- output res ---------------------------

# If the user doesn't input a voxel resolution for output, then the
# spatial resolution will be 1mm iso, or the smallest dimension of an
# input voxel (so the user doesn't lose resolution by default).

if ( "$RES" == "" ) then

    set RES = "1.0"
    set dim = `3dinfo -ad3 "$it1w"`

    foreach dd ( $dim )
        set comp = `echo "$RES > $dd" | bc`
        if ( "$comp" != "0" ) then
            echo "++ Found smaller vox dim: $dd < ${RES} "
            set RES = "$dd"
        endif
    end

    echo "++ Note: final volume grid resolution will be $RES mm isotropic."

endif

# ========================= output fnames ==========================

set ocmd   = "${opref}_cmd.txt"      # name for output command
set oanat  = "${opref}.nii.gz"       # name for output anatomical
set omtrx  = "${opref}_map_anat.aff12.1D" # save matrix with other files

# =======================================================================
# =========================== ** PROCESS ** =============================
# =======================================================================

echo "\n-----> STARTING $this_prog ---->"

# ---------------------------- CMD ---------------------------------

echo "\n\nThis command:"
echo "$this_prog $argv\n\n"

if ( "$cmd_file" == "" ) then
    set cmd_file = "$odir/$ocmd"
endif

# copy original command:
# dump copy of command into workdir/..
if ( $output_cmd == 1 ) then
    echo "++ Echoing the command to: $cmd_file"

    set rec_afni_ver = `afni -ver`
    echo "### AFNI version:"  > $cmd_file
    echo "# $rec_afni_ver\n"            >> $cmd_file

    echo "### Executed from the directory location:"  >> $cmd_file
    echo "# $here\n"            >> $cmd_file
    echo "### The command was:" >> $cmd_file
    echo "# $this_prog $argv"   >> $cmd_file
    echo "\n"                   >> $cmd_file
endif

# ------------------------- earliest QC: check overlap
# [PT: Oct 28, 2020] added in

@djunct_overlap_check                             \
    -ulay   "$it1w"                               \
    -olay   "$it2w"                               \
    -prefix "${odirqc}/init_qc_00_overlap_ut1w_ot2w"

# --------------------- start proc ---------------------------

# copy the t1w file to wdir via resampling: match t2w-t1w orients
set orient_t2w = `3dinfo -orient "$it2w"`
set t1res      = $wdir/f00_t1res.nii
3dresample                          \
    -overwrite                      \
    -prefix $t1res                  \
    -orient $orient_t2w             \
    -inset "$it1w"

if ( $status ) then
    goto BAD_EXIT
endif

# prep this side for aligning with a bit of unifizing and brightness
# thresholding
echo "++ Unifize T1w proc vol."
set t1uni = $wdir/f01_t1uni.nii
3dUnifize                            \
    -GM                              \
    -prefix  $t1uni                  \
    -input   "$t1res"                \
    -overwrite 

echo "++ Threshold T1w proc vol."
set t1thr = $wdir/f02_t1thr.nii
set v1 = `3dBrickStat -non-zero -automask   \
            -percentile 95 1 95 ${t1uni}`
echo ${v1[2]}
3dcalc -a  ${t1uni}                         \
       -expr "maxbelow(${v1[2]},a)"         \
       -prefix $t1thr                       \
       -overwrite

if ( $status ) then
    goto BAD_EXIT
endif

set t1align   = $wdir/f05_t1align.nii
if ( "$imtrx" == "" ) then

    echo "++ Calculating alignment (t2w -> t1w) matrix."

    set mt_t2_t1 = $wdir/map_t2w_to_t1w.aff12.1D
    set mt_t1_t2 = $wdir/map_t1w_to_t2w.aff12.1D    

    # [PT: Feb 7, 2018] might need if T1w has lots of non-brain stuff
    if ( "$DO_SS_T1W" == "1" ) then

        3dSkullStrip                           \
            -overwrite                         \
            -orig_vol                          \
            -prefix $wdir/f03_t1w_ss.nii.gz    \
            -input "$t1thr"

        if ( $status ) then
            goto BAD_EXIT
        endif

        set t1thr = "$wdir/f03_t1w_ss.nii.gz"

    endif

    if ( "$it2w_mask" == "" ) then
        set it2w_mask = "-source_automask"
    else
        set it2w_mask = "-source_mask $it2w_mask"
    endif

    3dAllineate   -echo_edu               \
        -1Dmatrix_save  $mt_t2_t1         \
        -prefix         "$t1align"        \
        -source         "$it2w"           \
        -base           "$t1thr"          \
        -twopass                          \
        -cost lpc                         \
        $cmass                            \
        $warp_opt                         \
        $it2w_mask                        \
        -wtprefix      $wdir/wtvol.nii.gz \
        -autoweight                       \
        -final wsinc5                     \
        ${maxrot} ${maxshf}               \
        -overwrite
    
    if ( $status ) then
        goto BAD_EXIT
    endif

    cat_matvec -ONELINE                   \
        $mt_t2_t1 -I                      \
        > $mt_t1_t2

    if ( $DO_T2W_TO_T1W ) then
        echo "++ Calculated the alignment matrix (t2w -> t1w)."
        set imtrx = $mt_t2_t1
    else
        # get inverse
        echo "++ Calculated the alignment matrix (t1w -> t2w)."
        cat_matvec -ONELINE                   \
            $mt_t2_t1 -I                      \
            > $mt_t1_t2
        set imtrx = $mt_t1_t2
    endif
endif

echo "++ Apply alignment matrix."

if ( $DO_T2W_TO_T1W ) then
    # T2w -> T1w

    set anat_mapped = $wdir/f06_t2mapped.nii
    if ( $OUT_T1W_GRID == 1 ) then
        echo "++ Outputting final image to t1w space"
        3dAllineate   -echo_edu               \
            -1Dmatrix_apply   $imtrx          \
            -source           "$it2w"         \
            -master           "$t1res"        \
            -prefix           $anat_mapped    \
            -float                            \
            -final wsinc5
        if ( $status ) then
            goto BAD_EXIT
        endif
    else
        echo "++ Outputting final image to $RES mm iso voxels"
        3dAllineate   -echo_edu               \
            -1Dmatrix_apply   $imtrx          \
            -source           "$it2w"         \
            -prefix           $anat_mapped    \
            -base             "$t1res"        \
            -float                            \
            -newgrid         "$RES"           \
            -final wsinc5
        if ( $status ) then
            goto BAD_EXIT
        endif
    endif

    # prepare for making images: ulay=t2w, olay=t1w
    set ulay = "$anat_mapped"
    set olay = "$t1res"
else
    # T1w -> T2w

    set anat_mapped = $wdir/f06_t1mapped.nii
    if ( $OUT_T2W_GRID == 1 ) then
        echo "++ Outputting final image to t2w space"
        3dAllineate   -echo_edu               \
            -1Dmatrix_apply   $imtrx          \
            -source           "$t1res"        \
            -master           "$it2w"         \
            -prefix           $anat_mapped    \
            -float                            \
            -final wsinc5
        if ( $status ) then
            goto BAD_EXIT
        endif
    else
        echo "++ Outputting final image to $RES mm iso voxels"
        3dAllineate   -echo_edu               \
            -1Dmatrix_apply   $imtrx          \
            -source           "$t1res"        \
            -prefix           $anat_mapped    \
            -base             "$it2w"         \
            -float                            \
            -newgrid         "$RES"           \
            -final wsinc5
        if ( $status ) then
            goto BAD_EXIT
        endif
    endif

    # prepare for making images: ulay=t2w, olay=t1w
    set ulay = "$it2w"
    set olay = "$anat_mapped"
endif

# copy output matrix
\cp $imtrx $odir/$omtrx

# get info to report below
set fdims = `3dinfo -n4  ${anat_mapped}`
set fres  = `3dinfo -ad3 ${anat_mapped}`

echo "++ Done workin'.  Just need to copy the final dset."

3dcalc                         \
    -overwrite                 \
    -a       $anat_mapped      \
    -expr    'a'               \
    -prefix  $odir/$oanat

if ( $status ) then
    goto BAD_EXIT
endif

if ( $DO_VIEWER == 1 ) then

    set vedge = $wdir/f10_edges.nii

    echo "++ More QC images: b0 on initial ref."
    3dedge3                         \
        -overwrite                  \
        -prefix $vedge              \
        -input  $olay

    if ( $qc_prefix == "" ) then
        set vpref0 = ${opref}_qc_t2w_t1wE
        set vpref1 = ${opref}_qc_t2w_t1w
    else
        set vpref0 = ${qc_prefix}_qc_t2w_t1wE
        set vpref1 = ${qc_prefix}_qc_t2w_t1w
    endif

    echo "\n\n"
    echo "++ QC image 00 ($odir/$oanat edges on $it2w): $vpref0"
    echo "\n\n"
    # need to put '[0]' on $iref?
    $my_viewer                            \
        -ulay "$ulay"                     \
        -ulay_range "2%" "98%"            \
        -olay "$vedge"                    \
        -func_range_perc_nz 50            \
        -pbar_posonly                     \
        -cbar "red_monochrome"            \
        -opacity 6                        \
        -prefix "$odirqc/$vpref0"         \
        -montx 5 -monty 3                 \
        -set_xhairs OFF                   \
        -label_mode 1 -label_size 4       \
        -do_clean 

    echo "\n\n"
    echo "++ QC image 01 ($odir/$oanat olay on $it2w): $vpref1"
    echo "\n\n"

    $my_viewer                            \
        -ulay "$ulay"                     \
        -ulay_range "2%" "98%"            \
        -olay "$olay"                     \
        -pbar_posonly                     \
        -opacity 4                        \
        -prefix "$odirqc/$vpref1"         \
        -montx 5 -monty 3                 \
        -set_xhairs OFF                   \
        -label_mode 1 -label_size 4       \
        -do_clean 

endif

# clean, by default
if ( "$DO_CLEAN" == "1" ) then
    echo "\n++ Cleaning working directory\n"
    \rm -rf $wdir
else
    echo "\n++ NOT removing working directory '$wdir'.\n"
endif

# final messages
echo ""
echo "++ The final data set is here:  $odir/$oanat"
echo "++ The final spatial dims are:  $fdims[1] $fdims[2] $fdims[3]"
echo "++ The spatial resolution is:   $fres[1] $fres[2] $fres[3]"
echo ""

# ---------------------------------------------------------------------

goto GOOD_EXIT

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

SHOW_HELP:
cat << EOF
# -----------------------------------------------------------------------

 This program is for aligning a T1w anatomical to a T2w anatomical
 using solid body parameters (i.e., only translation and rotation);
 this program does not clean or alter the T1w's brightness values
 (beyond minor smoothing from regridding).  If one is going to be
 entering the T1w volume into the realms of FreeSurfer (FS), one might
 want to do this just *before* that step.  If one wants axialized (or
 AC-PC-ized) anatomical volumes, one could perform that step on the
 T2w volume *before* using this function.

 This program mainly assumes that the T1w and T2w volume come from the
 same subject, and have similar contrasts expected for standard
 sequences and healthy adult brains.  This might still work for other
 applications, but caveat emptor (even more than usual).  This would
 *not* be recommended for aligning brains that aren't from the same
 subject.

 As part of this alignment, the T1w volume will end up with the same
 orientation and a similar FOV as the T2w volume.  Additionally, by
 default, the anatomical will be prepped a bit with an eye toward
 using FS, to have properties favorable to using it: 

   + the T1w volume is resampled to isotropic spatial resolution of
     either 1 mm voxel edges or, if the input volume has any edge
     length smaller than this, to that value (i.e., resampled to 1 mm
     or the minimum input voxel edge length, whichever is less).  The
     user can adjust this with the '-newgrid ...' option, or decide to 
     match the grid of the T2w volume via '-out_t2w_grid'. 

   + the T1w will have a FOV matching or quite similar to the T2w
     volume (as well as matching orientation).

   + [Aug 9, 2021] no longer checking about all even row
     dimensions---turned out not to be necessary.

 Note that, if you are preparing to use FS afterwards, then make sure
 to use their current help files, archives, etc. for all options and
 settings.  For example, while at present (March, 2017) FS does seem
 to prefer isotropic voxels with 1 mm edge length by default, one can
 use high resolution options for data acquired at higher resolution.
 Anyways, you can read more about that there.

  Ver. $version (PA Taylor, ${rev_dat})

# ----------------------------------------------------------------------

  OUTPUT:

     + NIFTI file: aligned T1w volume

     + QC snapshots of the T1w volume overlaying the T2w volume, and
       also the T1w edges overlaying the T2w volume.

# ----------------------------------------------------------------------

  RUNNING:

    $this_prog \
      -in_t1w  T1W                          \
      -in_t2w  T2W                          \
      -prefix  PPP                          \
      {-newgrid RES}                        \
      {-out_t2w_grid}                       \
      {-in_t2w_mask MASK_T2W}               \
      {-do_ss_tmp_t1w}                      \
      {-matrix MMM}                         \
      {-workdir WWW}                        \
      {-warp WAR}                           \
      {-no_cmd_out}                         \
      {-no_clean} 

  where:

   -in_t1w  T1W   :T1w volume (required).
   -in_t2w  T2W   :T2w volume (required; preferably from same subject as
                   T1W).

   -prefix  PPP   :output prefix for files and snapshots (required).

   -newgrid RES   :specify output T1w volume's final resolution; will be
                   isotropic in this value (default: 1 mm, or smallest voxel
                   edge length of input T1W if that value is < 1 mm).

   -out_t2w_grid  :final T1w volume is on the T2w volume's grid

   -out_t1w_grid  :when using '-align_t2w_to_t1w', one can use this opt
                   so the final T2w volume is on the T1W volume's grid

   -in_t2w_mask MASK_T2W
                  :can input a mask to apply to the t2w volume for
                   alignment purposes; might help in times of aligning 
                   hardship.

   -do_use_cmass DUC :should '-cmass' be used during the 3dAllineate's 
                   alignment estimation?  Valid values for DUC are:
                      Yes, 1 (to use -cmass)
                      No, 0 (to not use -cmass)
                   In many cases nowadays, it might make sense to disable
                   using -cmass, esp. if the initial T1w-T2w overlap is good.
                   (def: ${do_use_cmass})

   -do_ss_tmp_t1w :during an intermediate step, apply skullstripping
                   to the T1w volume-- final output is *not*
                   skullstripped.  This might be useful if there is
                   lots of non-brain tissue still in the T1w volume.

   -align_t2w_to_t1w :by default, this program outputs the alignment between
                   the anatomicals in the direction T1w->T2w.  Using this 
                   option reverses that order, so T2w->T1w. Note that all
                   QC images keep the T1w and T2w datasets in their same
                   overlay/underlay roles even if this opt is used.

   -warp WAR      :can choose which of the possible affine degrees of freedom
                   are employed in the warping, selecting them in the same
                   manner described in 3dAllineate's help;  that is, WAR can
                   be any of shift_only, shift_rotate, shift_rotate_scale, or
                   affine_general.  Default: WAR = shift_rotate.

   -maxrot MR     :use this option to change the range of rotations
                   searched during the 3dAllineate alignment process;
                   this option is just passed along to 3dAllineate,
                   and if not used, the default range from 3dAllineate
                   is used (see its help for the details)

   -maxshf MS     :use this option to change the range of shifts/translation
                   searched during the 3dAllineate alignment process;
                   this option is just passed along to 3dAllineate,
                   and if not used, the default range from 3dAllineate
                   is used (see its help for the details)

   -matrix MMM    :one can apply a pre-made matrix that has been made by
                   3dAllineate previously.  With this option.  If you want.
                   Using this opt means that the alignment between the input
                   T1w and T2w dsets is not calculated anew here, but instead
                   MMM is used instead.                   

   -workdir WWW   :specify a working directory, which can be removed;
                   (default name = '$wdir')

   -no_cmd_out    :don't save the command line call of this program
                   and the location where it was run (otherwise, it is
                   saved by default in the ODIR/).                     
   -no_clean      :no not delete temporary working directory (default is to 
                   remove it to save disk space).

   -echo          :run very, very verbosely

# ----------------------------------------------------------------------

  EXAMPLE

    # 1) have isotropic 1x1x1 mm final anat:
    $this_prog  \
        -in_t1w    MPRAGE.nii.gz        \
        -in_t2w    T2w_anat.nii.gz      \
        -newgrid   1.0                  \
        -prefix    t1w_alnd

    # 2) match the final anat resolution to that of the t2w dset:
    $this_prog  \
        -in_t1w    MPRAGE.nii.gz        \
        -in_t2w    T2w_anat.nii.gz      \
        -out_t2w_grid                   \
        -prefix    t1w_alndb

    # 3) instead of aligning the t1w dset to the t2w, have the t2w dset
    #    align to the t1w one (and also match final resolution of the t1w);
    #    also, disable using '-cmass' during the alignment, which assumes
    #    and/or makes use of good initial alignment
    $this_prog  \
        -in_t1w    MPRAGE.nii.gz        \
        -in_t2w    T2w_anat.nii.gz      \
        -out_t1w_grid                   \
        -align_t2w_to_t1w               \
        -do_use_cmass  0                \
        -prefix    t2w_alndb

# -----------------------------------------------------------------------

EOF

    goto GOOD_EXIT

SHOW_VERSION:
    echo "version  $version (${rev_dat})"
    goto GOOD_EXIT

FAIL_MISSING_ARG:
    echo "** ERROR: Missing an argument after option flag: '$argv[$ac]'"
    goto BAD_EXIT

BAD_EXIT:
    exit 1

# send everyone here, in case there is any cleanup to do
GOOD_EXIT:
    exit 0
