#!/usr/bin/env tcsh

# Make AFNI's help readable in a text editor again!
@global_parse `basename $0` "$*" ; if ($status) exit 0

# --------------------- revision history -------------------------
#
#set version   = "1.0";
#set rev_dat   = "some prior event in the spacetime continuum"
#   + RWC started and developed this program
#
#set version   = "1.1"; set rev_dat   = "May 18, 2018"
#   + PT started optionizing this program
#
#set version   = "1.1"; set rev_dat   = "July 1, 2018"
#   + fixed searching for path of dset, use explicit check first
#
#set version   = "1.2"; set rev_dat   = "July 5, 2018"
#   + set check if OMP_NUM_THREADS is explicitly set, so echoing don't
#     break it -- thanks P. Molfese et al.
#
#set version   = "1.3"; set rev_dat   = "Jan 17, 2019"
#   + [PT] bug fix-- 3danisosmooth cmd was missing ${odir} on I/O
#
#set version   = "1.4"; set rev_dat   = "Feb 11, 2019"
#   + [PT] new opt: turn unifize off (e.g., it 'twere already done
#     before)
#
#set version   = "1.41"; set rev_dat   = "Feb 11, 2019"
#   + [PT] new opt: can turn skullstrip and/or anisosmooth off (e.g.,
#     it 'twere already done before)
#
#set version   = "1.42"; set rev_dat   = "Feb 11, 2019"
#   + [PT] rename the 'turn off skull strip' to a less potentially
#     confusing name: '-init_skullstr_off'.  (So it doesn't falsely
#     seem like *no* skullstripping at all would be done.)
#
#set version   = "1.43"; set rev_dat   = "Feb 21, 2019"
#   + [PT] include '-Urad 30' in unifizing part, for improving
#     final output-- thanks for careful code reading, Yoichi
#
#set version   = "1.44"; set rev_dat   = "Mar 27, 2019"
#   + [RWC] add '-SSopt' option for adding options to 3dSkullStrip
#           (per the request of Allison Nugent)
#
#set version   = "1.45"; set rev_dat   = "Mar 29, 2019"
#   + [RWC] if SubID is a dataset name, funny things happened;
#           so edit out suffixes like '.nii', '.HEAD', etc.`
#
#set version   = "1.5"; set rev_dat   = "April 15, 2019"
#   + [PT] add in a couple new opts
#          "-giant_move": larger opening angle, etc.
#          "-deoblique":  can deoblique
# 
#set version   = "1.51"; set rev_dat   = "May 13, 2019"
#   + [PT] fixed help file for sphinxification: got rid of some
#          wandering "+" symbols in subheading titles
# 
#set version   = "1.52"; set rev_dat   = "June 18, 2019"
#   + [RWC] add 3dAutomask step to clean up some of the
#           little junk at the edge of the brain
# 
#set version   = "1.6"; set rev_dat   = "Jan 7, 2020"
#   + [PT] put ceiling (98%ile of non-zero vals) after 3danisosmooth
#        + also, use lpa+ZZ and lpa cost function for all rounds of align
#          -> should give better results, but is also slower (slightly)...
#        + add in options for putting in cost function values (backwards
#          compatibility): -cost_*
#        + opt for turn ceiling off: -ceil_off (backward compat)
#
#set version   = "1.61"; set rev_dat   = "Jan 20, 2020"
#   + [PT] new opt for deobliquing a la 3drefit-- probably would be more
#          useful than 3dWarp-style?
# 
#set version   = "2.0"; set rev_dat   = "Jan 27, 2020"
#   + [PT] This is a majorly new version of @SSwarper.
#        + final set of options from testing a loooot of things with
#          the mixed groups of 178 subj from different studies. This
#          set of options performed the best: fewest weird things,
#          sharpest mean across groups, and smallest stdev (both
#          inside and outside the brain).
#        + several opts have been added for additional control, as well
# 
#set version   = "2.1"; set rev_dat   = "Feb 20, 2020"
#   + [PT] Extra QC output:  QC*jpg montages
#        + one to check skullstripping, one to view warping in more detail
#
#set version   = "2.11"; set rev_dat   = "Feb 20, 2020"
#   + [PT] fix paths in the ulay/olay of extra QC*jpg images
#
#set version   = "2.2"; set rev_dat   = "Feb 27, 2020"
#   + [PT] new warpscale opt for 3dQwarp, courtesy of RWC
#
#set version   = "2.21"; set rev_dat   = "Sep 1, 2020"
#   + [PT] if '-skipwarp' is used, now set DO_EXTRA_WC = 0.
#
#set version   = "2.3"; set rev_dat   = "Sep 23, 2020"
#   + [PT] put ${status} checks to terminate on ~first failures
#
#set version   = "2.31"; set rev_dat   = "Oct 19, 2020"
#   + [PT] new qc image of initial overlap (with/without applying obl)
#
#set version   = "2.33"; set rev_dat   = "Dec 1, 2020"
#   + [PT] minorly change default "junk" file names 
#          - amusing anecdote: one AFNI user got an unexpected failure
#            when a random string started with "1D", so a file name
#            looked like "junk.SSwarper.1DfZbuq78qk.nii".  When 3dcalc
#            got hold of such a file, it complained about not being
#            able to read a *.1D file, because that combination
#            appears in the middle of a filename (???), but it output
#            correctly; however, 3dNwarpApply wanted none of it, and
#            produced a fatal error.
#          - Anyways, we can avoid this one in a bazillion hassleby
#            changing a dot to an underscore, so did.
#          - And added new "-tmp_nice_name" opt flag to have to a
#            simpler temp file name (use esp. when odir holds single
#            subj output)
#
#set version   = "2.4"; set rev_dat   = "Feb 5, 2020"
#   + [PT] more QC along the processing: 2 new init*jpg files
#        + add "-echo" as an opt
#
#set version   = "2.5"; set rev_dat   = "Feb 10, 2020"
#   + [PT] allow input mask
#        + try Strategy 1 of just replacing 3dSkullStrip with it
#
#set version   = "2.51"; set rev_dat   = "June 15, 2021"
#   + [PT] minor, minor replacements to avoid errors on old tcsh:
#        + "if(" -> "if ("
#        + ")then" -> ") then"
#
#set version   = "2.52";   set rev_dat   = "Sep 27, 2021"
#     + [PT] chauffeur label_size 3 -> 4, bc imseq.c shifted all sizes
#       down one level
#
#set version   = "2.53";   set rev_dat   = "Sep 28, 2021"
#     + [PT] cp initial/raw anat to $odir 
#          + also cp any initial/raw mask_ss to $odir 
#
#set version   = "2.54";   set rev_dat   = "Feb 1, 2022"
#     + [PT] Put clearer ERROR messaging in a couple spots
#          + remove any non-programmatic exclamation marks. those are BAD
#            to have in t/csh scripts, even in comments.
#
#set version   = "2.55";   set rev_dat   = "May 10, 2022"
#     + [PT] for 3dWarp-deobliqueing, use NN interp for mask, not wsinc5
#            -> thanks, RCR
#
set version   = "2.6";   set rev_dat   = "OCt 31, 2022"
#     + [PT] removing use of -cost_aff opt (it will just produce a warning
#            message).  It had been ignored, anyways.  Affine cost just
#            matches the cost_nl_init choice.
#
# ----------------------------------------------------------------

# some AFNI environment variables

setenv AFNI_DONT_LOGFILE  YES
setenv AFNI_COMPRESSOR    NONE

# set number of threads if run via SLURM

if ( $?SLURM_CPUS_PER_TASK ) then
 setenv OMP_NUM_THREADS $SLURM_CPUS_PER_TASK
else if ( $?NSLOTS ) then
 setenv OMP_NUM_THREADS $NSLOTS
endif

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

set this_prog = "@SSwarper"

set Adataset = ""    # req/ input dataset
set SubID    = ""    # req/ the subject ID
set Basedset = ""    # req/ reference dset- must have 4 bricks
set odir     = ""    # opt/ output dir

set minp     = "11"  # opt/ the minimum warp patch size

set btemplate = '$btemplate'
set tpath     = '$tpath'
set subj      = '$subj'
set str_msg   = '`@FindAfniDsetPath $btemplate`'

set liteopt   = "-lite"
set tightness = 0
set doclean   = 1
set verbopt   = ""
set skipwarp  = 0
set warpscale = 1        # def;  [Feb, 2020] RWC added opt to 3dQwarp 

set DO_UNIFIZE = 1
set DO_SKULLST = 1
set DO_ANISO   = 1
set DO_CEIL    = 1       # have a ceiling value on anat
set DO_GIANT   = 0
set DO_DEOB    = 0       # 3dWarp-style
set DO_DEOB_REF = 0      # 3drefit-style
set DO_EXTRA_QC = 1      # more chauffeur output 
set JUMP_TO_EXTRA_QC = 0 # just for testing/internal purposes
set DO_RANDOM_TMP = 1    # use random chars in tmp fname; can use nicer

set USER_COST_AFF = 0       # did user use -cost_aff?
set cost_nli   = "-lpa"     # used in:  3dQwarp, initial rounds
set cost_nlf   = "-pcl"     # used in:  3dQwarp, final rounds

set gaus_wt    = 4.5    # def for 3dQwarp

set saveall    = ""         # for 3dQwarp (def: don't);  just for debugging

set SSopt      = " "

set mask_ss    = ""

set more_echo  = ""

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

    # --------- required ---------------

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

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

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

    # --------- optional ---------------

    # [PT: Feb 25, 2020] RWC added this opt to 3dQwarp
    # lower warpscale -> less flexible warps
    # may be useful if odd bumps occur.  vals: [0.1, 1.0]
    else if ( "$argv[$ac]" == "-warpscale" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set warpscale = "$argv[$ac]"

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

    # [Jan 2019 - RWC]
    else if ( "$argv[$ac]" == "-nolite" ) then
        set liteopt = "-nolite"

    # [Jan 2019 - RWC]

    else if ( "$argv[$ac]" == "-tight" ) then
        @ tightness ++

    # [Jan 2019 - RWC]

    else if ( "$argv[$ac]" == "-noclean" ) then
        set doclean = 0

    # [PT: April 15, 2019]
    else if ( "$argv[$ac]" == "-giant_move" ) then
        set DO_GIANT = 1

    # [PT: April 15, 2019]
    else if ( "$argv[$ac]" == "-deoblique" ) then
        set DO_DEOB = 1

    # [PT: Jan 14, 2020] be able to purge obliquity info, too
    else if ( "$argv[$ac]" == "-deoblique_refitly" ) then
        set DO_DEOB_REF = 1

    # [Jan 2019 - RWC]

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

    # [PT: Feb 11, 2019]
    else if ( "$argv[$ac]" == "-unifize_off" ) then
        set DO_UNIFIZE = 0

    # [PT: Feb 11, 2019]
    else if ( "$argv[$ac]" == "-init_skullstr_off" ) then
        set DO_SKULLST = 0

    # [PT: Feb 11, 2019]
    else if ( "$argv[$ac]" == "-aniso_off" ) then
        set DO_ANISO = 0

    # [PT: Feb 11, 2019]
    else if ( "$argv[$ac]" == "-mask_ss" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set mask_ss = "$argv[$ac]"
        set DO_SKULLST = 0

    # [PT: Jan 31, 2020]
    else if ( "$argv[$ac]" == "-extra_qc_off" ) then
        set DO_EXTRA_QC = 0

    # just for internal running/testing purposes
    else if ( "$argv[$ac]" == "-jump_to_extra_qc" ) then
        set JUMP_TO_EXTRA_QC = 1

    # [PT: Feb 11, 2019]
    else if ( "$argv[$ac]" == "-ceil_off" ) then
        set DO_CEIL = 0

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

    # [PT: Jan 14, 2020] probably only for backwards compatibility:
    # default will be lpa+ZZ (could be 'hel' for back-compat)
    else if ( "$argv[$ac]" == "-cost_aff" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        ##set cost_aff = "$argv[$ac]"
        echo "+* WARNING ********"
        echo "   The option '-cost_aff' is no longer used.  The affine"
        echo "   alignment cost function will simply match the"
        echo "   '-cost_nl_init ..' value."
        set USER_COST_AFF = 1

    # [PT: Jan 14, 2020] probably just for backwards compatibility:
    # default will be lpa hereafter  (could be '-pcl' for back-compat)
    else if ( "$argv[$ac]" == "-cost_nl_init" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set cost_nli = "-$argv[$ac]"

    # [PT: Jan 15, 2020] separate early and later rounds of NL
    # warping, because using LPA for later rounds can be slooow
    else if ( "$argv[$ac]" == "-cost_nl_final" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set cost_nlf = "-$argv[$ac]"

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

    # [Jan 2019 - RWC]
    else if ( "$argv[$ac]" == "-verb" ) then
        set verbopt = "-verb"

    # [PT: Dec 1, 2020] changing default tmp fnames; this reverts to
    # old/ugly style
    else if ( "$argv[$ac]" == "-tmp_name_nice" ) then
        set DO_RANDOM_TMP = 0

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

    # 3dSkullStrip options for Allison [RWC: 27 Mar 2019]
    else if ( "$argv[$ac]" == "-SSopt" ) then
        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
        @ ac += 1
        set SSopt = "$argv[$ac]"

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

    # ---------- fin ------------------

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

    endif
    @ ac += 1
end

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

echo "++ Starting: $this_prog v$version"
if ( ! $?OMP_NUM_THREADS ) then
    set sss = `dirname $0`
    if ( -e $sss/afni_check_omp ) then
      set nnn = `$sss/afni_check_omp`
      echo "++ Default OMP_NUM_THREADS is $nnn"
      unset nnn
    else
      echo "++ OMP_NUM_THREADS is: not set by user, "
      echo "   so possibly just using one CPU, depending on system config"
    endif
    unset sss
else
    echo "++ OMP_NUM_THREADS set to $OMP_NUM_THREADS"
endif

if ( "$Adataset" == "" ) then
  echo "** ERROR: No -input option? :("
  goto BAD_EXIT
endif

if ( "$Basedset" == "" ) then
  echo "** ERROR No -base option? :("
  goto BAD_EXIT
endif

if ( ! -f "$Adataset" ) then
  set chad = `@CheckForAfniDset "$Adataset"`
  if ( "$chad" == "0" || "$chad" == "1" ) then
    echo "** ERROR: ${this_prog} -- Not finding dset $Adataset -- exiting :(("
    goto BAD_EXIT
  endif
endif

if ( "$SubID" == "" ) then
  echo "** ERROR ${this_prog} -- no subject ID entered with -subid? :("
  goto BAD_EXIT
endif

if ( "${mask_ss}" != "" ) then
    set check = `3dinfo -prefix "${mask_ss}"`
    if ( "$check" == "NO-DSET" ) then
        echo "** ERROR: can't find inset file ${mask_ss}"
        goto BAD_EXIT
    endif
endif

# edit SubID to remove any dataset suffixes [29 Mar 2019]
set sss = `basename "$SubID" .nii.gz`
set sss = `basename "$sss"   .nii`
set sss = `basename "$sss"   .`
set sss = `basename "$sss"   .HEAD`
set sss = `basename "$sss"   .BRIK.gz`
set sss = `basename "$sss"   .BRIK`
set sss = `basename "$sss"   +orig`
set sss = `basename "$sss"   +tlrc`

if ( "$sss" != "$SubID" ) then
  echo "++ Editing subject ID (-subid) $SubID to $sss"
  set SubID = "$sss"
endif
unset sss

# output dir from $Adataset, if not set by user
if ( "$odir" == "" ) then
    set odir = `dirname "$Adataset"`
    echo ""
    echo "++ Based on input, the output directory will be:"
    echo "     $odir"
    echo ""
endif

# [PT: Sep 1, 2020] skipwarp basically means a lot of other stuff gets
# skipped...  Have it turn off "extra" QC:
if ( ${skipwarp} ) then
    echo "+* Since we are skipping the warp, there is no 'extra' QC"
    set DO_EXTRA_QC = 0
endif

\mkdir -p "${odir}"

# set random prefix for temp files

### [PT: Dec 1, 2020] doesn't seem a need to have giant random string
### chars in temporary "junk" file names
if ( $DO_RANDOM_TMP ) then
    # original style, slightly modified bc the purely random chars
    # after '.' caaaan cause mischief (see amusing anecdote in
    # comments at top for this date).
    set pppp = "`3dnewid -fun11`"
    set pref = "$odir/junk.SSwarper_${pppp}_"
else
    set pref = "$odir/junk_ssw_"
endif 

# ---------------------------- this will change -------------------
#### ----> also add check that it has multiple subbricks

## find the base template dataset: start by seeing if it has been
## given directly (not just needing @Find*Path):
if ( -e "$Basedset" ) then
    set tpath = `dirname "$Basedset"`
else
    set tpath = `@FindAfniDsetPath "$Basedset"`
    if ( "$tpath" == '' ) then
        echo "** ERROR: ${this_prog} -- Failed to find template $Basedset"
        echo "   -- exiting :("
        goto BAD_EXIT
    endif
    set Basedset = $tpath/$Basedset
endif

## Require it to have enough bricks to really be a ref
set nvolbase = `3dinfo -nv "$Basedset"`
if ( $nvolbase < 5 ) then
    echo "** ERROR: Base $Basedset only has $nvolbase volumes:"
    echo "   to serve as a reference for $this_prog, it needs 4."
    echo "   See '$this_prog -help' -> 'The Template Dataset' for more info"
    goto BAD_EXIT
endif

# [PT: Oct 19, 2020] Add in new QC: initial source-base overlap
@djunct_overlap_check  ${more_echo}               \
    -ulay   "${Adataset}"                         \
    -olay   "$Basedset"                           \
    -prefix "${odir}/init_qc_00_overlap_usrc_obase"

if ( ${JUMP_TO_EXTRA_QC} ) then
    echo "\n++ Just jumpin' to extra QC\n"
    goto JUMP_TO_EXTRA_QC
endif

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

## Step #-1: copy the raw anat (and perhaps mask_ss) into output
## dir---often useful to compare later and/or check header info

set dset_cp = "$odir/anat_cp.${SubID}.nii"
if ( ! -f "${dset_cp}" ) then
    3dTcat                        \
        -prefix "${dset_cp}"      \
        "$Adataset"

    if ( ${status} ) then
        echo "** ERROR: in getting dset"
        goto BAD_EXIT
    endif
endif

if ( "${mask_ss}" != "" ) then
    set mask_ss_cp = "$odir/mask_ss_cp.${SubID}.nii"
    if ( ! -f "${mask_ss_cp}" ) then
        3dTcat                           \
            -prefix "${mask_ss_cp}"      \
            "${mask_ss}"

        if ( ${status} ) then
            echo "** ERROR: in getting dset"
            goto BAD_EXIT
        endif
    endif
endif

## start the work

## Step #0: Deoblique, either in style of 3dWarp or 3drefit
if ( ${DO_DEOB} ) then
    set dset_do = "$odir/anatDO.${SubID}.nii"
    3dWarp                        \
        -deoblique                \
        -wsinc5                   \
        -prefix "${dset_do}"      \
        "$Adataset"
    # ... and this will be new "input" dset
    set Adataset = "${dset_do}"

    if ( "${mask_ss}" != "" ) then
        set mask_do = "$odir/maskDO.${SubID}.nii.gz"
        # [PT: 10 May 2022] kernel from wsinc5 -> NN, as noticed by RCR
        3dWarp                        \
            -deoblique                \
            -NN                       \
            -prefix "${mask_do}"      \
            "${mask_ss}"
        # ... and this will be new "input" dset
        set mask_ss = "${mask_do}"
    endif 

else if ( ${DO_DEOB_REF} ) then
    set dset_do = "$odir/anatDO.${SubID}.nii"
    # copy
    3dcalc                        \
        -a "$Adataset"            \
        -expr 'a'                 \
        -prefix "${dset_do}"
    # purge obliquity info
    3drefit                       \
        -deoblique                \
        "${dset_do}"

    if ( ${status} ) then
        echo "** ERROR: program failed in Step deob"
        goto BAD_EXIT
    endif

    # ... and this will be new "input" dset
    set Adataset = "${dset_do}"

    if ( "${mask_ss}" != "" ) then
        set mask_do = "$odir/maskDO.${SubID}.nii.gz"

        # copy
        3dcalc                        \
            -a "${mask_ss}"            \
            -expr 'a'                 \
            -prefix "${mask_do}"
        # purge obliquity info
        3drefit                       \
            -deoblique                \
            "${mask_do}"

        if ( ${status} ) then
            echo "** ERROR: program failed in Step deob"
            goto BAD_EXIT
        endif

        # ... and this will be new "input" dset
        set mask_ss = "${mask_do}"
    endif 

endif

## Step #1: Unifize the input T1
echo "++ SSW Step 1"

# [PT: Feb 11, 2019] now allows for option *not* to unifize
# [PT: Feb 21, 2019] use Urad=30 when unifizing (better for final output)
if ( ! -f "$odir/anatU.${SubID}.nii" ) then
    if ( ${DO_UNIFIZE} ) then
        3dUnifize                             \
            -GM                               \
            -clfrac 0.4                       \
            -Urad 30                          \
            -prefix "$odir/anatU.$SubID.nii"  \
            -input "$Adataset"

        if ( ${status} ) then
            echo "** ERROR: program failed in Step U"
            goto BAD_EXIT
        endif
    else
        echo "++ NOT unifizing dset."
        3dcopy                                \
            "$Adataset"                       \
            "$odir/anatU.$SubID.nii"

        if ( ${status} ) then
            echo "** ERROR: program failed in Step U"
            goto BAD_EXIT
        endif

    endif
endif



# [PT: Feb 11, 2019] now allows for option *not* to anismoothize
if ( ! -f "$odir/anatUA.${SubID}.nii" ) then
    if ( ${DO_ANISO} ) then
        3danisosmooth                           \
            -iters 1                            \
            -3D                                 \
            -automask                           \
            -noneg                              \
            -prefix "$odir/anatUA.${SubID}.nii" \
            "$odir/anatU.${SubID}.nii"

        if ( ${status} ) then
            echo "** ERROR: program failed in Step UA"
            goto BAD_EXIT
        endif
    else
        echo "++ NOT anisosmoothing dset."
        3dcopy                                  \
            "$odir/anatU.$SubID.nii"            \
            "$odir/anatUA.$SubID.nii"

        if ( ${status} ) then
            echo "** ERROR: program failed in Step UA"
            goto BAD_EXIT
        endif
    endif
endif

# [PT: Jan 6, 2020] get rid of potentially *very* high vals: take
# 98%ile of non-zero vox
if ( ! -f "$odir/anatUAC.${SubID}.nii" ) then
    if ( ${DO_CEIL} ) then
        set vvv = `3dBrickStat                  \
                    -percentile 98 1 98         \
                    -non-zero                   \
                    "$odir/anatUA.${SubID}.nii"`
        # apply ceiling
        3dcalc                                  \
            -a "$odir/anatUA.${SubID}.nii"      \
            -expr "maxbelow(${vvv[2]},a)"       \
            -prefix "$odir/anatUAC.${SubID}.nii"

        if ( ${status} ) then
            echo "** ERROR: program failed in Step UAC"
            goto BAD_EXIT
        endif

    else
        echo "++ NOT anisosmoothing dset."
        3dcopy                                   \
            "$odir/anatUA.$SubID.nii"            \
            "$odir/anatUAC.$SubID.nii"

        if ( ${status} ) then
            echo "** ERROR: program failed in Step UAC"
            goto BAD_EXIT
        endif
    endif
endif


# Step #2: Strip Skull (Ziad's way)
echo "++ SSW Step 2"

# [PT: Feb 11, 2019] now allows for option *not* to skullstrip 
if ( ! -f "$odir/anatS.${SubID}.nii" ) then
    if ( ${DO_SKULLST} ) then
        3dSkullStrip \
            -input  "$odir/anatUAC.$SubID.nii" \
            -prefix "$odir/anatS.$SubID.nii"   \
            -debug 1 -ld 33 -niter 777         \
            -shrink_fac_bot_lim 0.777          \
            -exp_frac 0.0666 -orig_vol         \
            $SSopt

        if ( ${status} ) then
            echo "** ERROR: program failed in Step 2"
            goto BAD_EXIT
        endif
    else if ( "${mask_ss}" != "" ) then
        echo "++ Instead of skullstripping, apply mask_ss dset."
        3dcalc                                 \
            -a "$odir/anatUAC.$SubID.nii"      \
            -b ${mask_ss}                      \
            -expr 'a*step(b)'                  \
            -prefix "$odir/anatS.$SubID.nii"

        if ( ${status} ) then
            echo "** ERROR: program failed in Step 2"
            goto BAD_EXIT
        endif
    else
        echo "++ NOT pre-skullstripping dset."
        3dcopy                                 \
            "$odir/anatUAC.$SubID.nii"         \
            "$odir/anatS.$SubID.nii"

        if ( ${status} ) then
            echo "** ERROR: program failed in Step 2"
            goto BAD_EXIT
        endif
    endif
endif

## Step #3: run 3dQwarp first time to a moderate level (skull on)
echo "++ SSW Step 3"

if ( ${DO_GIANT} ) then
    3dQwarp -echo_edu        ${saveall}                   \
        $liteopt $verbopt                                 \
        -base   "${Basedset}[1]"                          \
        ${cost_nli}                                       \
        -warpscale ${warpscale}                           \
        -source "$odir/anatUAC.$SubID.nii"                \
        -weight "${Basedset}[2]"                          \
        -allineate -noneg -maxlev 5 -iwarp -awarp         \
        -wtgaus ${gaus_wt}                                \
        -inedge                                           \
        -workhard:3:5 -nopenalty                          \
        -prefix "${pref}TAL5.nii"                         \
        -allopt '-twobest 11 -twopass -maxrot 45 -maxshf 40 -source_automask+2 -cmass'

    if ( ${status} ) then
        echo "** ERROR: program failed in Step 3"
        goto BAD_EXIT
    endif
else
    3dQwarp -echo_edu           ${saveall}                \
        $liteopt $verbopt                                 \
        -base   "${Basedset}[1]"                          \
        ${cost_nli}                                       \
        -warpscale ${warpscale}                           \
        -source "$odir/anatUAC.$SubID.nii"                \
        -weight "${Basedset}[2]"                          \
        -allineate -noneg -maxlev 5 -iwarp -awarp         \
        -wtgaus ${gaus_wt}                                \
        -inedge                                           \
        -workhard:3:5 -nopenalty                          \
        -prefix "${pref}TAL5.nii"

    if ( ${status} ) then
        echo "** ERROR: program failed in Step 3"
        goto BAD_EXIT
    endif
endif

# -------- check early NL alignment to ref vol space ---------

set ulay  = "${pref}TAL5.nii"
set olay  = "${Basedset}"
set opref_qc_nl0 = "${odir}/init_qc_01_nl0.${SubID}.jpg"

@djunct_edgy_align_check ${more_echo}  \
    -montx            8                \
    -monty            1                \
    -ulay             "${ulay}"        \
    -olay             "${olay}"        \
    -box_focus_slices "${olay}"        \
    -prefix           "${pref}MONT3"

2dcat                                  \
    -gap 5                             \
    -gap_col 255 204 153               \
    -nx 1                              \
    -ny 3                              \
    -prefix "${opref_qc_nl0}"          \
    "${pref}"MONT3*jpg

## Step #4: mask off the skull using the template (second skull-strip)
echo "++ SSW Step 4"

3dmask_tool                        \
    -input "${Basedset}[3]"        \
    -dilate_input 2                \
    -prefix "${pref}MASK.nii"

3dcalc                             \
    -a "${pref}MASK.nii"           \
    -b "${pref}TAL5.nii"           \
    -expr 'step(a)*b'              \
    -prefix "${pref}TAL5mm.nii"

if ( ${status} ) then
    echo "** ERROR: program failed in Step 4"
    goto BAD_EXIT
endif

## Step #5: warp this masked dataset back to original space
echo "++ SSW Step 5"

3dNwarpApply -echo_edu                 \
    -nwarp "${pref}TAL5_WARPINV.nii"   \
    -master "$odir/anatS.$SubID.nii"   \
    -source "${pref}TAL5mm.nii"        \
    -prefix "${pref}TAL5ww.nii"

if ( ${status} ) then
    echo "** ERROR: program failed in Step 5a"
    goto BAD_EXIT
endif

## warp the mask itself (dilated) back to orig space

\rm -f "${pref}MASK.nii"
3dmask_tool  -echo_edu           \
    -input "${Basedset}[3]"      \
    -dilate_input 3              \
    -prefix "${pref}MASK.nii"

3dNwarpApply \
    -nwarp "${pref}TAL5_WARPINV.nii"   \
    -master "$odir/anatS.$SubID.nii"   \
    -source "${pref}MASK.nii"          \
    -prefix "${pref}MASKO.nii"         \
    -ainterp NN

if ( ${status} ) then
    echo "** ERROR: program failed in Step 5b"
    goto BAD_EXIT
endif


## merge these backward warped datasets with the 3dSkullStrip
## output to get a better original skull-stripped result

3dcalc -a "$odir/anatS.${SubID}.nii"  \
       -b "${pref}TAL5ww.nii"         \
       -c "${pref}MASKO.nii"          \
       -expr 'step(c)*max(a,b)'       \
       -prefix "$odir/anatSS.${SubID}.nii"

if ( ${status} ) then
    echo "** ERROR: program failed in Step 5c"
    goto BAD_EXIT
endif

# DRG's erode-dilate trick for cleanup of little crap [16 Jan 2019]
3dmask_tool                           \
    -dilate_inputs -3 3               \
    -prefix "${pref}de3.nii"          \
    -input  "$odir/anatSS.${SubID}.nii"
3dcalc                                \
    -a "$odir/anatSS.${SubID}.nii"    \
    -b "${pref}de3.nii"               \
    -expr 'a*step(b)'                 \
    -prefix "$odir/anatSSc.${SubID}.nii"

if ( ${status} ) then
    echo "** ERROR: program failed in Step 5d"
    goto BAD_EXIT
endif

# Throw in an automask step to clean up the outer edges [18 Jun 2019]
\rm -f "${pref}de3.nii" "$odir/anatSS.${SubID}.nii"
3dAutomask                                            \
    -apply_prefix "$odir/anatSSd.${SubID}.nii"        \
    "$odir/anatSSc.${SubID}.nii"
\mv -f "$odir/anatSSd.${SubID}.nii" "$odir/anatSS.${SubID}.nii"
\rm -f "$odir/anatSSc.${SubID}.nii"

if ( ${status} ) then
    echo "** ERROR: program failed in Step 5e"
    goto BAD_EXIT
endif


# don't do any more if skipping the final (precision) warp

if ( $skipwarp ) goto CLEANUP

## Step #6: affine transform that result to template space
echo "++ SSW Step 6"

# [PT: Jan 7, 2019] use lpa+ZZ cost func by default; also add in -twopass
# [PT: Oct 31, 2022] um, $cost_aff has no effect here
echo "++ SSW Step 6 ... piece 0"
3dAllineate \
    -1Dmatrix_apply "${pref}TAL5_Allin.aff12.1D"   \
    -source "$odir/anatSS.${SubID}.nii"            \
    -master "${pref}TAL5mm.nii"                    \
    #-cost ${cost_aff}                              \
    #-twopass                                       \
    #-autoweight -source_automask                   \
    -final wsinc5                                  \
    -prefix "${pref}AffSS.nii"

if ( ${status} ) then
    echo "** ERROR: program failed in Step 6.0"
    goto BAD_EXIT
endif

# -------- check aff alignment to ref vol space ---------

set ulay  = "${pref}AffSS.nii"
set olay  = "${Basedset}"
set opref_qc_aff = "${odir}/init_qc_02_aff.${SubID}.jpg"

@djunct_edgy_align_check               \
    -montx            8                \
    -monty            1                \
    -ulay             "${ulay}"        \
    -olay             "${olay}"        \
    -box_focus_slices "${olay}"        \
    -prefix           "${pref}MONT6"

2dcat                                  \
    -gap 5                             \
    -gap_col 255 204 153               \
    -nx 1                              \
    -ny 3                              \
    -prefix "${opref_qc_aff}"          \
    "${pref}"MONT6*jpg



## warp to template space (skull off),
## initializing using the previous 3dQwarp -awarp output

# Run 3dQwarp in several segments, to avoid gcc OpenMP bug
#  where it freezes sometimes with inter-thread conflicts;
#  this happens only in long runs, so breaking the run
#  into pieces will (I hope) elide this annoyance - RWC :(

# Piece number 1
echo "++ SSW Step 6 ... piece 1"
3dQwarp -echo_edu        ${saveall}      \
    $liteopt $verbopt                    \
    -base "${Basedset}[0]"               \
    -source "${pref}AffSS.nii"           \
    -iniwarp "${pref}TAL5_AWARP.nii"     \
    -warpscale ${warpscale}              \
    ${cost_nli}                          \
    -inilev 1 -maxlev 5                  \
    -wtgaus ${gaus_wt}                   \
    -inedge                              \
    -pblur -workhard:5:5                 \
    -nodset -prefix "${pref}QQ5.nii"

if ( ${status} ) then
    echo "** ERROR: program failed in Step 6.1"
    goto BAD_EXIT
endif

# Piece number 2
echo "++ SSW Step 6 ... piece 2"
3dQwarp -echo_edu      ${saveall}        \
    $liteopt $verbopt                    \
    -base "${Basedset}[0]"               \
    -source "${pref}AffSS.nii"           \
    -iniwarp "${pref}QQ5_WARP.nii"       \
    -warpscale ${warpscale}              \
    ${cost_nlf}                          \
    -inilev 6 -maxlev 7 -workhard:6:7    \
    -wtgaus ${gaus_wt}                   \
    -inedge                              \
    -pblur                               \
    -nodset -prefix "${pref}QQ7.nii"

if ( ${status} ) then
    echo "** ERROR: program failed in Step 6.2"
    goto BAD_EXIT
endif

if ( $minp > 13 ) then
  # Final piece for coarse final patch size
  echo "++ SSW Step 6 ... piece 3"
  3dQwarp -echo_edu      ${saveall}        \
      $liteopt $verbopt                    \
      -base "${Basedset}[0]"               \
      -source "${pref}AffSS.nii"           \
      -iniwarp "${pref}QQ7_WARP.nii"       \
      -warpscale ${warpscale}              \
      ${cost_nlf}                          \
      -inilev 8                            \
      -wtgaus ${gaus_wt}                   \
      -inedge                              \
      -pblur -minpatch $minp               \
      -Qfinal                              \
      -prefix "$odir/anatQQ.${SubID}.nii"

  if ( ${status} ) then
      echo "** ERROR: program failed in Step 6.3a"
      goto BAD_EXIT
  endif

else
  # Penultimate piece for refined final patch size
  echo "++ SSW Step 6 ... piece 3a"
  @ mpp = $minp + 6
  3dQwarp  -echo_edu       ${saveall}      \
      $liteopt $verbopt                    \
      -base "${Basedset}[0]"               \
      -source "${pref}AffSS.nii"           \
      -iniwarp "${pref}QQ7_WARP.nii"       \
      -warpscale ${warpscale}              \
      ${cost_nlf}                          \
      -inilev 8 -minpatch $mpp             \
      -wtgaus ${gaus_wt}                   \
      -inedge                              \
      -pblur                               \
      -nodset -prefix "${pref}QQ9.nii"

  if ( ${status} ) then
      echo "** ERROR: program failed in Step 6.3b"
      goto BAD_EXIT
  endif

  # Ultimate piece for refined final patch size
  echo "++ SSW Step 6 ... piece 3b"
  3dQwarp  -echo_edu       ${saveall}      \
      $liteopt $verbopt                    \
      -base "${Basedset}[0]"               \
      -source "${pref}AffSS.nii"           \
      -iniwarp "${pref}QQ9_WARP.nii"       \
      -warpscale ${warpscale}              \
      ${cost_nlf}                          \
      -inilev 10                           \
      -wtgaus ${gaus_wt}                   \
      -inedge                              \
      -pblur -minpatch $minp               \
      -prefix "$odir/anatQQ.${SubID}.nii"

  if ( ${status} ) then
      echo "** ERROR: program failed in Step 6.3c"
      goto BAD_EXIT
  endif

endif

echo "++ SSW: done warping.  Finalize."

# DRG's erode-dilate trick for cleanup of little stuff [16 Jan 2019]
3dmask_tool                             \
    -dilate_inputs -1 1                 \
    -prefix "${pref}de3.nii"            \
    -input  "$odir/anatQQ.${SubID}.nii"
3dcalc                                  \
    -a "$odir/anatQQ.${SubID}.nii"      \
    -b "${pref}de3.nii"                 \
    -expr 'a*step(b)'                   \
    -prefix "$odir/anatQQc.${SubID}.nii"

if ( ${status} ) then
    echo "** ERROR: program failed in Finalize"
    goto BAD_EXIT
endif

\rm -f "${pref}de3.nii" "anatQQ.${SubID}.nii"
\mv -f "$odir/anatQQc.${SubID}.nii" "$odir/anatQQ.${SubID}.nii"

## Step #7: Make two pretty pictures, and scram
echo "++ SSW Step 7"

@snapshot_volreg "$odir/anatQQ.${SubID}.nii" "$Basedset" "$odir/AM${SubID}"
@snapshot_volreg "$Basedset" "$odir/anatQQ.${SubID}.nii" "$odir/MA${SubID}"

JUMP_TO_EXTRA_QC:

if ( ${DO_EXTRA_QC} ) then

    # -------- check skullstripping in orig space ---------

    set ulay  = "${odir}/anatU.${SubID}.nii"
    set olay0 = "${odir}/anatSS.${SubID}.nii"
    set olay  = "${pref}mask.nii"
    set opref_exqc1 = "${odir}/QC_anatSS.${SubID}.jpg"

    3dcalc                                 \
        -overwrite                         \
        -a ${olay0}                        \
        -expr 'bool(a)'                    \
        -prefix ${olay}                    \
        -byte
        
    @chauffeur_afni                        \
        -ulay ${ulay}                      \
        -ulay_range "2%" "98%"             \
        -olay ${olay}                      \
        -func_range_perc_nz 1              \
        -set_subbricks 0 0 0               \
        -box_focus_slices ${olay}          \
        -cbar "Reds_and_Blues_Inv"         \
        -opacity 4                         \
        -prefix  "${pref}MONT1"            \
        -save_ftype JPEG                   \
        -montx 9 -monty 1                  \
        -montgap 3                         \
        -set_xhairs OFF                    \
        -label_mode 1 -label_size 4      

    2dcat                                  \
        -gap 5                             \
        -gap_col 66 184 254                \
        -nx 1                              \
        -ny 3                              \
        -prefix "${opref_exqc1}"           \
        "${pref}"MONT1*jpg

    # -------- check alignment in ref vol space ---------

    set ulay  = "${odir}/anatQQ.${SubID}.nii"
    set olay  = "${Basedset}"
    set opref_exqc2 = "${odir}/QC_anatQQ.${SubID}.jpg"

    @djunct_edgy_align_check               \
        -montx            8                \
        -monty            1                \
        -ulay             "${ulay}"        \
        -olay             "${olay}"        \
        -box_focus_slices "${olay}"        \
        -prefix           "${pref}MONT2"

    2dcat                                  \
        -gap 5                             \
        -gap_col 66 184 254                \
        -nx 1                              \
        -ny 3                              \
        -prefix "${opref_exqc2}"           \
        "${pref}"MONT2*jpg

endif

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

## Clean up the junk

CLEANUP:

## Rename affine warp matrix
echo "++ SSW cleanup"

# [PT: Feb 19, 2020] Put the mv cmd in an IF condition, because with
# re-running extra QC images, might not have the first file name;
# now, avoid an unnecessary error msg.
if ( -e "${pref}TAL5_Allin.aff12.1D" ) then
    \mv -f "${pref}TAL5_Allin.aff12.1D" "$odir/anatQQ.${SubID}.aff12.1D"
endif

if ( $doclean ) then
  \rm -f "${pref}"*
endif

echo "--------------------------------------------------------"
echo "++ DONE with SSW.  Check QC images:"
printf "   $odir/AM${SubID}.jpg  $odir/MA${SubID}.jpg "
if ( ${DO_EXTRA_QC} ) then
printf " ${opref_exqc1}  ${opref_exqc2}"
endif 
echo ""
echo ""
echo "   In ref/base space (check warping):"
echo "     ulay=anat, olay=base edges  :  $odir/AM${SubID}.jpg"
echo "     ulay=base edges, olay=anat  :  $odir/MA${SubID}.jpg"
if ( ${DO_EXTRA_QC} ) then
    echo "     ulay=anat, olay=base edges  :  ${opref_exqc2}"
    echo ""
    echo "   In native space (check skullstripping):"
    echo "     ulay=anat, olay=SS mask     :  ${opref_exqc1}"
endif
echo ""
echo "--------------------------------------------------------"

goto GOOD_EXIT

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

SHOW_HELP:

cat <<EOF

OVERVIEW ~1~

This script has dual purposes for processing a given subject's
anatomical volume:
    + to skull-strip the brain, and
    + to calculate the warp to a reference template/standard space.
Automatic snapshots of the registration are created, as well, to help
the QC process.

This program cordially ties in directly with afni_proc.py, so you can
run it beforehand, check the results, and then provide both the
skull-stripped volume and the warps to the processing program.  That
is convenient.

*** This program has been superseded by the newer version, sswarper2,
    which has essentially the same syntax and usage as this one. ***

Current version = $version
Authorship      = Bob, Bob, there is one Bob, He spells it B-O-B.

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

USAGE ~1~

    ${this_prog}             \
        -input  AA           \
        -base   BB           \
        -subid  SS           \
       {-odir   OD}          \
       {-minp   MP}          \
       {-nolite}             \
       {-skipwarp}           \
       {-unifize_off}        \
       {-init_skullstr_off}  \
       {-extra_qc_off}       \
       {-jump_to_extra_qc}   \
       {-cost_nl_init CNI}   \
       {-cost_nl_final CNF}  \
       {-deoblique}          \
       {-deoblique_refitly}  \
       {-warpscale WS}       \
       {-SSopt 'strings'     \
       {-aniso_off}          \
       {-ceil_off}           \
       {-tmp_name_nice}      \
       {-echo}               \
       {-verb}               \
       {-noclean}

where (note: many of the options with 'no' and 'off' in their name are
really just included for backwards compatibility, as this program has
grown/improved over time):

  -input  AA :(req) an anatomical dataset, *not* skull-stripped, with
              resolution about 1 mm.

  -base   BB :(req) a base template dataset, with contrast similar to
              the input AA dset, probably from some kind of standard
              template.
              NB: this dataset is not *just* a standard template,
              because it is not a single volume-- read about its
              composition in the NOTES on the 'The Template Dataset',
              below.
              The program first checks if the dset BB exists as
              specified; if not, then if just the filename has been
              provided it searches the AFNI_GLOBAL_SESSION,
              AFNI_PLUGINPATH, and afni bin directory (in that order)
              for the named dataset.

  -subid  SS :(req) name code for output datasets (e.g., 'sub007').

  -odir   OD :(opt) output directory for all files from this program
              (def: directory of the '-input AA').

  -minp   MP :(opt) minimum patch size on final 3dQwarp (def: 11).

  -nolite    :(opt) Do not use the '-lite' option with 3dQwarp;
              This option is used for backward compatibility, if you want
              to run 3dQwarp the same way as older versions of @SSwarper.
              The new way (starting Jan 2019) is to use the '-lite'
              option with 3dQwarp to speed up the calculations.
              (def: use '-lite' for faster calculations).

  -skipwarp  :(opt) Do not compute past the output of anatSS.{subid}.nii.
              This option is used if you just want the skull-stripped
              result in original coordinates, without the warping
              to the template space (anatQQ). The script will run faster.

  -deoblique :(opt) apply obliquity information to deoblique the input
              volume ('3dWarp -deoblique -wsinc5 ...'), as an initial step.
              This might introduce the need to overcome a large rotation
              during the alignment, though

  -deoblique_refitly :(opt) purge obliquity information to deoblique
              the input volume (copy, and then '3drefit -deoblique ...'), 
              as an initial step.  This might help when data sets are
              very... oblique.

  -warpscale WS :(opt) opt to control flexibility of warps in 3dQwarp and
              how they adjust with patch size;  see 3dQwarp's help for 
              more info. Allowed values of WS are in range [0.1, 1.0].
              (def: 1.0)

  -giant_move :(opt) when starting the initial alignment to the template,
              apply the same parameter expansions to 3dAllineate that
              align_epi_anat.py does with the same option flag.  This
              might be useful if the brain has a very large angle away
              from "typical" ones, etc.

  -unifize_off :(opt) don't start with a 3dUnifize command to try reduce
              effects of brightness inhomogeneities.  Probably only
              useful if unifizing has been previously performed on the
              input dset.

  -aniso_off :(opt) don't preprocess with a 3danisosmooth command to
              try reduce effects of weird things (in a technical
              sense).  Possible that this will never be used in the
              history of running this program.

  -ceil_off  :(opt) by default, after anisosmoothing, this program
              will apply put a ceiling on values in the dset, to get rid
              of possible outliers (ceil = 98%ile of non-zero voxels in
              the whole volume).  This option will turn that off.

  -init_skullstr_off :(opt) don't preprocess with a 3dSkullstrip command
              to roughly isolated brain in the beginning.  This might
              be useful with macaque dsets.

  -extra_qc_off :(opt) don't make extra QC images QC*jpg (for some
              unknown reason).

  -jump_to_extra_qc :(opt) just make the two QC*jpg images from a
              previous run of @SSwarper.  These QC*jpg images are new
              QC output (as of late Feb, 2020), so this might be
              useful to add a quick check to previously run data.
              This command would just be tacked on to previously
              executed one.

  #-cost_aff CA :***no longer used.*** The affine cost function is only  
              set via cost_nl_init, since the affine alignment is just a
              'preliminary alignment' for that one.  So, what is specified
              for the cost_nl_init will be used for the affine.

  -cost_nl_init CNI 
             :(opt) specify cost function for initial nonlinear
              (3dQwarp) part of alignment.  Here, 'CNI' would be the
              cost function name to be provided (def: is now "lpa").
              This is probably only here for backwards compatibility
              to older @SSwarper (where def was 'pcl').

  -cost_nl_final CNF 
             :(opt) specify cost function for final nonlinear
              (3dQwarp) parts of alignment.  Here, 'CNF' would be the
              cost function to be provided (def: is now "pcl").  This
              is separate from the initial nonlinear warp cost values
              '-cost_nl_init ..', because using those here might be
              pretty slow; however, using "lpa" here might help
              results.

  -SSopt 'strings' :(opt) The content of 'strings' (which should be
              in quotes if there are any blanks) is copied to the
              end of the 3dSkullStrip command line. Example:
                -SSopt '-o_ply Fred.Is.Wonderful'
              to have 3dSkullStrip produce a .ply surface file
              as an additional output.

  -mask_ss MSS :(opt) as an alternative to skullstripping at an early
              stage, you can provide a mask to be used before the
              initial affine alignment.  The mask MSS can come from
              anywhere, but @SUMA_Make_Spec_FS now makes a convenient
              one from the FS parcellation (though it would have to be
              resampled to the input anatomical's grid).

  -tmp_name_nice :(opt) default temporary "junk.*" filenames include
              a large, random char string.  This is ugly, but useful
              if outputting several different SSW runs into the same
              directory that we intermediate files (very likely) don't
              get overwritten.  However, if you prefer, you can use a
              nicer, non-random intermediate file prefix: "junk_ssw".
              I would use this when the output dir ("-odir ..")
              doesn't contain multiple SSW outputs.

  -verb      :(opt) Apply the '-verb' option to 3dQwarp, to get more
              verbose progress information - mostly used for debugging.

  -echo      :(opt) Run the script with "set echo", for extra verbosity
              in the terminal output.  Mainly for debugging times. 

  -noclean   :(opt) Do not delete the 'junk' files at the end of
              computations - mostly used for debugging and testing.

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

REFERENCE DATASETS ~1~

If you are reading this message, then several reference data sets
(base volumes) for @SSwarper now exist within the AFNI realm. Oh, what
a time it is to be alive.  A current list includes:

+ MNI152_2009_template_SSW.nii.gz
+ TT_N27_SSW.nii.gz
+ HaskinsPeds_NL_template1.0_SSW.nii.gz

Some of these are distributed with the AFNI binaries, and other may be
found online. You can make other reference base templates in whatever
space you prefer, but note that it must have several subvolumes of
information included-- see NOTES on the 'The Template Dataset', below
(which also contains a link to the @SSwarper template tutorial online
help).

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

OUTPUTS ~1~

Datasets ~2~

  Suppose the -prefix is 'sub007' (because you scanned Bond, JamesBond?).
  Then the outputs from this script will be"

  anatDO.sub007.nii       = deobliqued version of original dataset;
                            (*only if* using '-deoblique' opt);
  anatU.sub007.nii        = intensity uniform-ized original dataset
                            (or, if '-unifize_off' used, a copy of orig dset);
  anatUA.sub007.nii       = anisotropically smoothed version of the above
                            (or, if '-aniso_off' used, a copy of anatU.*.nii)
  anatUAC.sub007.nii      = ceiling-capped ver of the above (at 98%ile of 
                            non-zero values)
                            (or, if '-ceil_off' used, a copy of anatUA.*.nii)

  anatS.sub007.nii        = first pass skull-stripped original dataset
                            (or, if '-init_skullstr_off' used, a copy of
                            anatUAC.*.nii);
  anatSS.sub007.nii       = second pass skull-stripped original dataset;
                            * note that anatS and anatSS are 'original'
                              in the sense that they are aligned with
                              the input dataset - however, they have been
                              unifized and weakly smoothed: they are
                              stripped versions of anatUAC; if you want
                              a skull-stripped copy of the input with
                              no other processing, use a command like
                                3dcalc -a INPUTDATASET        \
                                       -b anatSS.sub007.nii   \
                                       -expr 'a*step(b)'      \
                                       -prefix anatSSorig.sub007.nii

  anatQQ.sub007.nii       = skull-stripped dataset nonlinearly warped to
                            the base template space;
  anatQQ.sub007.aff12.1D  = affine matrix to transform original dataset
                            to base template space;
  anatQQ.sub007_WARP.nii  = incremental warp from affine transformation
                            to nonlinearly aligned dataset;

  * The .aff12.1D and _WARP.nii transformations need to be catenated to get
    the full warp from original space to the base space; example:

      3dNwarpApply -nwarp 'anatQQ.sub007_WARP.nii anatQQ.sub007.aff12.1D' ...

QC images ~2~

  AMsub007.jpg            = 3x3 snapshot image of the anatQQ.sub007.nii
                            dataset with the edges from the base template
                            overlaid -- to check the alignment;
  MAsub007.jpg            = similar to the above, with the roles of the
                            template and the anatomical datasets reversed.
  QC_anatQQ.sub007.jpg    = like AM*.jpg, but 3 rows of 8 slices
  QC_anatSS.sub007.jpg    = check skullstripping in orig space: ulay is
                            input dset, and olay is mask of
                            skullstripped output (anatSS* dset)

  init_qc_00_overlap_uinp_obase.jpg
    o [ulay] original source dset
      [olay] original base dset
    o single image montage to check initial overlap of source and base,
      ignoring any obliquity that might be present (i.e., the way AFNI 
      GUI does by default, and also how alignment starts)
    o if initial overlap is not strong, alignment can fail or
      produce weirdness
    o *if* either dset has obliquity, then an image of both after 
      deobliquing with 3dWarp is created (*DEOB.jpg), and a text file
      about obliquity is also created (*DEOB.txt).

* It is important to examine (at least) the two .jpg snapshot images to
  make sure that the skull-stripping and nonlinear warping worked well.

USING SSW WITH AFNI_PROC.PY ~1~

  When B-O-B uses @SSwarper for skull-stripping plus warping, He gives
  afni_proc.py these options (among others, hence the ellipses), after
  running ${this_prog} successfully. Here, 'subj' is the subject
  identifier:

  |    set template = MNI152_2009_template_SSW.nii.gz
  |    
  |    afni_proc.py                                                  \
  |      ...                                                         \
  |      -copy_anat             anatSS.\${subj}.nii                   \
  |      -anat_has_skull        no                                   \
  |      -align_opts_aea        -cost lpc+ZZ -giant_move             \
  |                             -check_flip                          \
  |      -volreg_align_to       MIN_OUTLIER                          \
  |      -volreg_align_e2a                                           \
  |      -volreg_tlrc_warp                                           \
  |      -tlrc_base             \${template}                          \
  |      -tlrc_NL_warp                                               \
  |      -tlrc_NL_warped_dsets  anatQQ.\${subj}.nii                   \
  |                             anatQQ.\${subj}.aff12.1D              \
  |                             anatQQ.\${subj}_WARP.nii
  |      ...

NOTES ~1~

The Template dataset ~2~

  Any reference base template dataset, such as
  MNI152_2009_template_SSW.nii.gz, must have the first *4* volumes here
  (and can have the optional 5th for later uses, as described):
    [0] = skull-stripped template brain volume
    [1] = skull-on template brain volume
    [2] = weight mask for nonlinear registration, with the
          brain given greater weight than the skull
    [3] = binary mask for the brain
    [4] = binary mask for gray matter plus some CSF (slightly dilated)
          ++ this volume is not used in this script
          ++ it is intended for use in restricting FMRI analyses
             to the 'interesting' parts of the brain
          ++ this mask should be resampled to your EPI spatial
             resolution (see program 3dfractionize), and then
             combined with a mask from your experiment reflecting
             your EPI brain coverage (see program 3dmask_tool).

  More information about making these (with scripts) is provided on
  the Interweb:
    https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/template_atlas/sswarper_base.html

The steps being run ~2~

  You Know My Methods, Watson...

  #1: Uniform-ize the input dataset's intensity via 3dUnifize.
       ==> anatU.sub007.nii
  #2: Strip the skull with 3dSkullStrip, with mildly aggressive settings.
       ==> anatS.sub007.nii
  #3: Nonlinearly warp (3dQwarp) the result from #1 to the skull-on
      template, driving the warping to a medium level of refinement.
  #4: Use a slightly dilated brain mask from the template to
      crop off the non-brain tissue resulting from #3 (3dcalc).
  #5: Warp the output of #4 back to original anatomical space,
      along with the template brain mask, and combine those
      with the output of #2 to get a better skull-stripped
      result in original space (3dNwarpApply and 3dcalc).
       ==> anatSS.sub007.nii
  #6  Restart the nonlinear warping, registering the output
      of #5 to the skull-off template brain volume (3dQwarp).
       ==> anatQQ.sub007.nii (et cetera)
  #7  Use @snapshot_volreg to make the pretty pictures.
       ==> AMsub007.jpg and MAsub007.jpg

Temporary files ~2~

  If the script crashes for some reason, it might leave behind files
  whose names start with 'junk' -- you should delete these files
  manually.

WHAT TO DO IF RESULTS ARE WAY OFF? ~1~

  The importance of initial dset overlap ~2~

  Always, always, always check the initial image made by SSW when it
  runs:  

      init_qc_00_overlap_uinp_obase.jpg

  This image tells you how well your datasets overlap initially before
  the alignment work begins. **The better the overlap, the lower the
  chance that something weird happens in your output.** All the SSW
  templates have reasonable coordinates, meaning that (x, y, z) = (0,
  0, 0) is in a good spot for it.  If there is poor overlap, probably
  your input dataset has weird/bad coordinates for some reason.  

  You can use @Align_Centers to put your anatomical dset in a better
  spot (though note, if you are going to be processing EPI data
  afterwards, you will want to move that along, as well, perhaps as a
  "child" dataset).

  By far the most common problem leading to obviously bad outputs is
  that the initial datasets are waaay far apart when they start, and
  the program gets stuck in a false minimum of solutions.

  Other issues ~2~

  Sometimes, it can be hard to separate the brain from dura and/or
  skull surrounding the brain.  If little bits are left around in the
  masking images, then perhaps adding one of the following options for
  will help (this can help the initial skullstripping):

    -SSopt '-blur_fwhm 2'
    -SSopt '-blur_fwhm 3'

  Any other questions/oddities, please don't hesitate to inquire on
  the AFNI Message Board

EXAMPLES ~1~

*** This program has been superseded by the newer version, sswarper2,
    which has essentially the same syntax and usage as this one. ***


  1) Run the program, deciding what the main output directory will be
  called (e.g., based on the subject ID):

    @SSwarper                                    \
        -input  anat_t1w.nii.gz                  \
        -base   MNI152_2009_template_SSW.nii.gz  \
        -subid  sub-001                          \
        -odir   group/o.aw_sub-001

  2) Same as above, but since we are using one outdir per subject, use
  more aesthetically pleasing names of temporary files (which get
  deleted, anyways):

    @SSwarper                                    \
        -tmp_name_nice                           \
        -input  anat_t1w.nii.gz                  \
        -base   MNI152_2009_template_SSW.nii.gz  \
        -subid  sub-001                          \
        -odir   group/o.aw_sub-001

  3) As of version 2.5, you can input a mask to be used instead of
  skullstripping.  For example, a good one might be the
  parcellation-derived (but filled in) mask from @SUMA_Make_Spec_FS
  after running FS's recon-all (though you will have to resample it
  from the FS output grid to that of your input anatomical):

    @SSwarper                                     \
        -tmp_name_nice                            \
        -input   anat_t1w.nii.gz                  \
        -mask_ss fs_parc_wb_mask_RES.nii.gz       \
        -base    MNI152_2009_template_SSW.nii.gz  \
        -subid   sub-001                          \
        -odir    group/o.aw_sub-001

EOF

if ( ${USER_COST_AFF} == 1 ) then
    echo "+* WARNING: the user used the '-cost_aff ..' option,"
    echo "   but it is no longer used within the program---it has no effect."
    echo "   The affine cost function just follows what is set by"
    echo "   The '-cost_nl_init ..' option value (or its default)."
endif


   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

GOOD_EXIT:
    exit 0