1#!/bin/tcsh
2
3# Make AFNI's help readable in a text editor again!
4@global_parse `basename $0` "$*" ; if ($status) exit 0
5
6# --------------------- revision history -------------------------
7#
8#set version   = "1.0";
9#set rev_dat   = "some prior event in the spacetime continuum"
10#   + RWC started and developed this program
11#
12#set version   = "1.1"; set rev_dat   = "May 18, 2018"
13#   + PT started optionizing this program
14#
15#set version   = "1.1"; set rev_dat   = "July 1, 2018"
16#   + fixed searching for path of dset, use explicit check first
17#
18#set version   = "1.2"; set rev_dat   = "July 5, 2018"
19#   + set check if OMP_NUM_THREADS is explicitly set, so echoing don't
20#     break it -- thanks P. Molfese et al.!
21#
22#set version   = "1.3"; set rev_dat   = "Jan 17, 2019"
23#   + [PT] bug fix-- 3danisosmooth cmd was missing ${odir} on I/O
24#
25#set version   = "1.4"; set rev_dat   = "Feb 11, 2019"
26#   + [PT] new opt: turn unifize off (e.g., it 'twere already done
27#     before)
28#
29#set version   = "1.41"; set rev_dat   = "Feb 11, 2019"
30#   + [PT] new opt: can turn skullstrip and/or anisosmooth off (e.g.,
31#     it 'twere already done before)
32#
33#set version   = "1.42"; set rev_dat   = "Feb 11, 2019"
34#   + [PT] rename the 'turn off skull strip' to a less potentially
35#     confusing name: '-init_skullstr_off'.  (So it doesn't falsely
36#     seem like *no* skullstripping at all would be done.)
37#
38#set version   = "1.43"; set rev_dat   = "Feb 21, 2019"
39#   + [PT] include '-Urad 30' in unifizing part, for improving
40#     final output-- thanks for careful code reading, Yoichi!
41#
42#set version   = "1.44"; set rev_dat   = "Mar 27, 2019"
43#   + [RWC] add '-SSopt' option for adding options to 3dSkullStrip
44#           (per the request of Allison Nugent)
45#
46#set version   = "1.45"; set rev_dat   = "Mar 29, 2019"
47#   + [RWC] if SubID is a dataset name, funny things happened;
48#           so edit out suffixes like '.nii', '.HEAD', etc.`
49#
50#set version   = "1.5"; set rev_dat   = "April 15, 2019"
51#   + [PT] add in a couple new opts
52#          "-giant_move": larger opening angle, etc.
53#          "-deoblique":  can deoblique
54#
55#set version   = "1.51"; set rev_dat   = "May 13, 2019"
56#   + [PT] fixed help file for sphinxification: got rid of some
57#          wandering "+" symbols in subheading titles
58#
59#set version   = "1.52"; set rev_dat   = "June 18, 2019"
60#   + [RWC] add 3dAutomask step to clean up some of the
61#           little junk at the edge of the brain
62#
63#set version   = "1.6"; set rev_dat   = "Jan 7, 2020"
64#   + [PT] put ceiling (98%ile of non-zero vals) after 3danisosmooth
65#        + also, use lpa+ZZ and lpa cost function for all rounds of align
66#          -> should give better results, but is also slower (slightly)...
67#        + add in options for putting in cost function values (backwards
68#          compatibility): -cost_*
69#        + opt for turn ceiling off: -ceil_off (backward compat)
70#
71#set version   = "1.61"; set rev_dat   = "Jan 20, 2020"
72#   + [PT] new opt for deobliquing a la 3drefit-- probably would be more
73#          useful than 3dWarp-style?
74#
75#set version   = "2.0"; set rev_dat   = "Jan 27, 2020"
76#   + [PT] This is a majorly new version of @SSwarper.
77#        + final set of options from testing a loooot of things with
78#          the mixed groups of 178 subj from different studies. This
79#          set of options performed the best: fewest weird things,
80#          sharpest mean across groups, and smallest stdev (both
81#          inside and outside the brain).
82#        + several opts have been added for additional control, as well
83#
84#set version   = "2.1"; set rev_dat   = "Feb 20, 2020"
85#   + [PT] Extra QC output:  QC*jpg montages
86#        + one to check skullstripping, one to view warping in more detail
87#
88#set version   = "2.11"; set rev_dat   = "Feb 20, 2020"
89#   + [PT] fix paths in the ulay/olay of extra QC*jpg images
90#
91#set version   = "2.2"; set rev_dat   = "Feb 27, 2020"
92#   + [PT] new warpscale opt for 3dQwarp, courtesy of RWC
93#
94#set version   = "2.21"; set rev_dat   = "Sep 1, 2020"
95#   + [PT] if '-skipwarp' is used, now set DO_EXTRA_WC = 0.
96#
97#set version   = "2.3"; set rev_dat   = "Sep 23, 2020"
98#   + [PT] put ${status} checks to terminate on ~first failures
99#
100#set version   = "2.31"; set rev_dat   = "Oct 19, 2020"
101#   + [PT] new qc image of initial overlap (with/without applying obl)
102#
103#set version   = "2.33"; set rev_dat   = "Dec 1, 2020"
104#   + [PT] minorly change default "junk" file names
105#          - amusing anecdote: one AFNI user got an unexpected failure
106#            when a random string started with "1D", so a file name
107#            looked like "junk.SSwarper.1DfZbuq78qk.nii".  When 3dcalc
108#            got hold of such a file, it complained about not being
109#            able to read a *.1D file, because that combination
110#            appears in the middle of a filename (!?!?), but it output
111#            correctly; however, 3dNwarpApply wanted none of it, and
112#            produced a fatal error.
113#          - Anyways, we can avoid this one in a bazillion hassleby
114#            changing a dot to an underscore, so did.
115#          - And added new "-tmp_nice_name" opt flag to have to a
116#            simpler temp file name (use esp. when odir holds single
117#            subj output)
118#
119#set version   = "2.4"; set rev_dat   = "Feb 5, 2020"
120#   + [PT] more QC along the processing: 2 new init*jpg files
121#        + add "-echo" as an opt
122#
123#set version   = "2.5"; set rev_dat   = "Feb 10, 2020"
124#   + [PT] allow input mask
125#        + try Strategy 1 of just replacing 3dSkullStrip with it
126#
127#set version   = "2.51"; set rev_dat   = "June 15, 2021"
128#   + [PT] minor, minor replacements to avoid errors on old tcsh:
129#        + "if(" -> "if ("
130#        + ")then" -> ") then"
131#
132#set version   = "2.52";   set rev_dat   = "Sep 27, 2021"
133#     + [PT] chauffeur label_size 3 -> 4, bc imseq.c shifted all sizes
134#       down one level
135#
136set version   = "2.53";   set rev_dat   = "Sep 28, 2021"
137#     + [PT] cp initial/raw anat to $odir
138#          + also cp any initial/raw mask_ss to $odir
139#
140# ----------------------------------------------------------------
141
142# some AFNI environment variables
143
144setenv AFNI_DONT_LOGFILE  YES
145setenv AFNI_COMPRESSOR    NONE
146
147# set number of threads if run via SLURM
148
149if ( $?SLURM_CPUS_PER_TASK ) then
150 setenv OMP_NUM_THREADS $SLURM_CPUS_PER_TASK
151else if ( $?NSLOTS ) then
152 setenv OMP_NUM_THREADS $NSLOTS
153endif
154
155# ===================================================================
156
157set this_prog = "@SSwarper"
158
159set Adataset = ""    # req/ input dataset
160set SubID    = ""    # req/ the subject ID
161set Basedset = ""    # req/ reference dset- must have 4 bricks
162set odir     = ""    # opt/ output dir
163
164set minp     = "11"  # opt/ the minimum warp patch size
165
166set btemplate = '$btemplate'
167set tpath     = '$tpath'
168set subj      = '$subj'
169set str_msg   = '`@FindAfniDsetPath $btemplate`'
170
171set liteopt   = "-lite"
172set tightness = 0
173set doclean   = 1
174set verbopt   = ""
175set skipwarp  = 0
176set warpscale = 1        # def;  [Feb, 2020] RWC added opt to 3dQwarp
177
178set DO_UNIFIZE = 1
179set DO_SKULLST = 1
180set DO_ANISO   = 1
181set DO_CEIL    = 1       # have a ceiling value on anat
182set DO_GIANT   = 0
183set DO_DEOB    = 0       # 3dWarp-style
184set DO_DEOB_REF = 0      # 3drefit-style
185set DO_EXTRA_QC = 1      # more chauffeur output
186set JUMP_TO_EXTRA_QC = 0 # just for testing/internal purposes
187set DO_RANDOM_TMP = 1    # use random chars in tmp fname; can use nicer
188
189set cost_aff   = "lpa+ZZ"   # used in:  3dAllineate -cost ...
190set cost_nli   = "-lpa"     # used in:  3dQwarp, initial rounds
191set cost_nlf   = "-pcl"     # used in:  3dQwarp, final rounds
192
193set gaus_wt    = 4.5    # def for 3dQwarp
194
195set saveall    = ""         # for 3dQwarp (def: don't);  just for debugging
196
197set SSopt      = " "
198
199set mask_ss    = ""
200
201# ------------------- process options, a la rr ----------------------
202
203if ( $#argv == 0 ) goto SHOW_HELP
204
205set ac = 1
206while ( $ac <= $#argv )
207    # terminal options
208    if ( ("$argv[$ac]" == "-h" ) || ("$argv[$ac]" == "-help" )) then
209        goto SHOW_HELP
210    endif
211    if ( "$argv[$ac]" == "-ver" ) then
212        goto SHOW_VERSION
213    endif
214
215    # --------- required ---------------
216
217    if ( "$argv[$ac]" == "-input" ) then
218        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
219        @ ac += 1
220        set Adataset = "$argv[$ac]"
221
222    else if ( "$argv[$ac]" == "-subid" ) then
223        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
224        @ ac += 1
225        set SubID = "$argv[$ac]"
226
227    else if ( "$argv[$ac]" == "-base" ) then
228        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
229        @ ac += 1
230        set Basedset = "$argv[$ac]"
231
232    # --------- optional ---------------
233
234    # [PT: Feb 25, 2020] RWC added this opt to 3dQwarp
235    # lower warpscale -> less flexible warps
236    # may be useful if odd bumps occur.  vals: [0.1, 1.0]
237    else if ( "$argv[$ac]" == "-warpscale" ) then
238        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
239        @ ac += 1
240        set warpscale = "$argv[$ac]"
241
242    # min patch
243    else if ( "$argv[$ac]" == "-minp" ) then
244        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
245        @ ac += 1
246        set minp = "$argv[$ac]"
247
248    # [Jan 2019 - RWC]
249    else if ( "$argv[$ac]" == "-nolite" ) then
250        set liteopt = "-nolite"
251
252    # [Jan 2019 - RWC]
253
254    else if ( "$argv[$ac]" == "-tight" ) then
255        @ tightness ++
256
257    # [Jan 2019 - RWC]
258
259    else if ( "$argv[$ac]" == "-noclean" ) then
260        set doclean = 0
261
262    # [PT: April 15, 2019]
263    else if ( "$argv[$ac]" == "-giant_move" ) then
264        set DO_GIANT = 1
265
266    # [PT: April 15, 2019]
267    else if ( "$argv[$ac]" == "-deoblique" ) then
268        set DO_DEOB = 1
269
270    # [PT: Jan 14, 2020] be able to purge obliquity info, too
271    else if ( "$argv[$ac]" == "-deoblique_refitly" ) then
272        set DO_DEOB_REF = 1
273
274    # [Jan 2019 - RWC]
275
276    else if ( "$argv[$ac]" == "-skipwarp" ) then
277        set skipwarp = 1
278
279    # [PT: Feb 11, 2019]
280    else if ( "$argv[$ac]" == "-unifize_off" ) then
281        set DO_UNIFIZE = 0
282
283    # [PT: Feb 11, 2019]
284    else if ( "$argv[$ac]" == "-init_skullstr_off" ) then
285        set DO_SKULLST = 0
286
287    # [PT: Feb 11, 2019]
288    else if ( "$argv[$ac]" == "-aniso_off" ) then
289        set DO_ANISO = 0
290
291    # [PT: Feb 11, 2019]
292    else if ( "$argv[$ac]" == "-mask_ss" ) then
293        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
294        @ ac += 1
295        set mask_ss = "$argv[$ac]"
296        set DO_SKULLST = 0
297
298    # [PT: Jan 31, 2020]
299    else if ( "$argv[$ac]" == "-extra_qc_off" ) then
300        set DO_EXTRA_QC = 0
301
302    # just for internal running/testing purposes
303    else if ( "$argv[$ac]" == "-jump_to_extra_qc" ) then
304        set JUMP_TO_EXTRA_QC = 1
305
306    # [PT: Feb 11, 2019]
307    else if ( "$argv[$ac]" == "-ceil_off" ) then
308        set DO_CEIL = 0
309
310    else if ( "$argv[$ac]" == "-saveall" ) then
311        set saveall = "-saveall"
312
313    # [PT: Jan 14, 2020] probably only for backwards compatibility:
314    # default will be lpa+ZZ (could be 'hel' for back-compat)
315    else if ( "$argv[$ac]" == "-cost_aff" ) then
316        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
317        @ ac += 1
318        set cost_aff = "$argv[$ac]"
319
320    # [PT: Jan 14, 2020] probably just for backwards compatability:
321    # default will be lpa hereafter  (could be '-pcl' for back-compat)
322    else if ( "$argv[$ac]" == "-cost_nl_init" ) then
323        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
324        @ ac += 1
325        set cost_nli = "-$argv[$ac]"
326
327    # [PT: Jan 15, 2020] separate early and later rounds of NL
328    # warping, because using LPA for later rounds can be slooow
329    else if ( "$argv[$ac]" == "-cost_nl_final" ) then
330        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
331        @ ac += 1
332        set cost_nlf = "-$argv[$ac]"
333
334    else if ( "$argv[$ac]" == "-wtgaus" ) then
335        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
336        @ ac += 1
337        set gaus_wt = "$argv[$ac]"
338
339    # [Jan 2019 - RWC]
340    else if ( "$argv[$ac]" == "-verb" ) then
341        set verbopt = "-verb"
342
343    # [PT: Dec 1, 2020] changing default tmp fnames; this reverts to
344    # old/ugly style
345    else if ( "$argv[$ac]" == "-tmp_name_nice" ) then
346        set DO_RANDOM_TMP = 0
347
348    # output dir
349    else if ( "$argv[$ac]" == "-odir" ) then
350        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
351        @ ac += 1
352        set odir = "$argv[$ac]"
353
354    # 3dSkullStrip options for Allison [RWC: 27 Mar 2019]
355    else if ( "$argv[$ac]" == "-SSopt" ) then
356        if ( $ac >= $#argv ) goto FAIL_MISSING_ARG
357        @ ac += 1
358        set SSopt = "$argv[$ac]"
359
360    else if ( "$argv[$ac]" == "-echo" ) then
361        set echo
362
363    # ---------- fin ------------------
364
365    else
366        echo "** unexpected option #$ac = '$argv[$ac]'"
367        goto BAD_EXIT
368
369    endif
370    @ ac += 1
371end
372
373# ====================================================================
374
375echo "++ Starting: $this_prog v$version"
376if ( ! $?OMP_NUM_THREADS ) then
377    set sss = `dirname $0`
378    if ( -e $sss/afni_check_omp ) then
379      set nnn = `$sss/afni_check_omp`
380      echo "++ Default OMP_NUM_THREADS is $nnn"
381      unset nnn
382    else
383      echo "++ OMP_NUM_THREADS is: not set by user, "
384      echo "   so possibly just using one CPU, depending on system config"
385    endif
386    unset sss
387else
388    echo "++ OMP_NUM_THREADS set to $OMP_NUM_THREADS"
389endif
390
391if ( "$Adataset" == "" ) then
392  echo "***** No -input option? :("
393  goto BAD_EXIT
394endif
395
396if ( ! -f "$Adataset" ) then
397  set chad = `@CheckForAfniDset "$Adataset"`
398  if ( "$chad" == "0" || "$chad" == "1" ) then
399    echo "***** ${this_prog} -- Not finding dset $Adataset -- exiting :(("
400    goto BAD_EXIT
401  endif
402endif
403
404if ( "$SubID" == "" ) then
405  echo "***** ${this_prog} -- no subject ID entered with -subid? :("
406  goto BAD_EXIT
407endif
408
409if ( "${mask_ss}" != "" ) then
410    set check = `3dinfo -prefix "${mask_ss}"`
411    if ( "$check" == "NO-DSET" ) then
412        echo "** ERROR: can't find inset file ${mask_ss}"
413        goto BAD_EXIT
414    endif
415endif
416
417# edit SubID to remove any dataset suffixes [29 Mar 2019]
418set sss = `basename "$SubID" .nii.gz`
419set sss = `basename "$sss"   .nii`
420set sss = `basename "$sss"   .`
421set sss = `basename "$sss"   .HEAD`
422set sss = `basename "$sss"   .BRIK.gz`
423set sss = `basename "$sss"   .BRIK`
424set sss = `basename "$sss"   +orig`
425set sss = `basename "$sss"   +tlrc`
426
427if ( "$sss" != "$SubID" ) then
428  echo "++ Editing subject ID (-subid) $SubID to $sss"
429  set SubID = "$sss"
430endif
431unset sss
432
433# output dir from $Adataset, if not set by user
434if ( "$odir" == "" ) then
435    set odir = `dirname $Adataset`
436    echo ""
437    echo "++ Based on input, the output directory will be:"
438    echo "     $odir"
439    echo ""
440endif
441
442# [PT: Sep 1, 2020] skipwarp basically means a lot of other stuff gets
443# skipped...  Have it turn off "extra" QC:
444if ( ${skipwarp} ) then
445    echo "+* Since we are skipping the warp, there is no 'extra' QC"
446    set DO_EXTRA_QC = 0
447endif
448
449\mkdir -p ${odir}
450
451# set random prefix for temp files
452
453### [PT: Dec 1, 2020] doesn't seem a need to have giant random string
454### chars in temporary "junk" file names
455if ( $DO_RANDOM_TMP ) then
456    # original style, slightly modified bc the purely random chars
457    # after '.' caaaan cause mischief (see amusing anecdote in
458    # comments at top for this date).
459    set pppp = "`3dnewid -fun11`"
460    set pref = $odir/junk.SSwarper_${pppp}_
461else
462    set pref = $odir/junk_ssw_
463endif
464
465# ---------------------------- !! this will change !! -------------------
466#### ----> also add check that it has multiple subbricks
467
468## find the base template dataset: start by seeing if it has been
469## given directly (not just needing @Find*Path):
470if ( -e "$Basedset" ) then
471    set tpath = `dirname "$Basedset"`
472else
473    set tpath = `@FindAfniDsetPath "$Basedset"`
474    if ( "$tpath" == '' ) then
475        echo "** ERROR: ${this_prog} -- Failed to find template $Basedset"
476        echo "   -- exiting :("
477        goto BAD_EXIT
478    endif
479    set Basedset = $tpath/$Basedset
480endif
481
482## Require it to have enough bricks to really be a ref
483set nvolbase = `3dinfo -nv "$Basedset"`
484if ( $nvolbase < 5 ) then
485    echo "** Base $Basedset only has $nvolbase volumes:"
486    echo "   to serve as a reference for $this_prog, it needs 4!"
487    echo "   See '$this_prog -help' -> 'The Template Dataset' for more info"
488    goto BAD_EXIT
489endif
490
491# [PT: Oct 19, 2020] Add in new QC: initial source-base overlap
492@djunct_overlap_check                             \
493    -ulay   "${Adataset}"                         \
494    -olay   "$Basedset"                           \
495    -prefix ${odir}/init_qc_00_overlap_usrc_obase
496
497if ( ${JUMP_TO_EXTRA_QC} ) then
498    echo "\n++ Just jumpin' to extra QC\n"
499    goto JUMP_TO_EXTRA_QC
500endif
501
502# --------------------------------------------------------------------
503
504## Step #-1: copy the raw anat (and perhaps mask_ss) into output
505## dir---often useful to compare later and/or check header info
506
507set dset_cp = $odir/anat_cp."${SubID}".nii
508if ( ! -f "${dset_cp}" ) then
509    3dTcat                        \
510        -prefix "${dset_cp}"      \
511        "$Adataset"
512endif
513
514if ( "${mask_ss}" != "" ) then
515    set mask_ss_cp = $odir/mask_ss_cp."${SubID}".nii
516    if ( ! -f "${mask_ss_cp}" ) then
517        3dTcat                           \
518            -prefix "${mask_ss_cp}"      \
519            "${mask_ss}"
520    endif
521endif
522
523## start the work
524
525## Step #0: Deoblique, either in style of 3dWarp or 3drefit
526if ( ${DO_DEOB} ) then
527    set dset_do = $odir/anatDO."${SubID}".nii
528    3dWarp                        \
529        -deoblique                \
530        -wsinc5                   \
531        -prefix "${dset_do}"      \
532        "$Adataset"
533    # ... and this will be new "input" dset
534    set Adataset = "${dset_do}"
535
536    if ( "${mask_ss}" != "" ) then
537        set mask_do = $odir/maskDO."${SubID}".nii.gz
538        3dWarp                        \
539            -deoblique                \
540            -wsinc5                   \
541            -prefix "${mask_do}"      \
542            "${mask_ss}"
543        # ... and this will be new "input" dset
544        set mask_ss = "${mask_do}"
545    endif
546
547else if ( ${DO_DEOB_REF} ) then
548    set dset_do = $odir/anatDO."${SubID}".nii
549    # copy
550    3dcalc                        \
551        -a "$Adataset"            \
552        -expr 'a'                 \
553        -prefix "${dset_do}"
554    # purge obliquity info
555    3drefit                       \
556        -deoblique                \
557        "${dset_do}"
558    # ... and this will be new "input" dset
559    set Adataset = "${dset_do}"
560
561    if ( "${mask_ss}" != "" ) then
562        set mask_do = $odir/maskDO."${SubID}".nii.gz
563
564        # copy
565        3dcalc                        \
566            -a "${mask_ss}"            \
567            -expr 'a'                 \
568            -prefix "${mask_do}"
569        # purge obliquity info
570        3drefit                       \
571            -deoblique                \
572            "${mask_do}"
573        # ... and this will be new "input" dset
574        set mask_ss = "${mask_do}"
575    endif
576
577endif
578
579if ( ${status} ) then
580    echo "** ERROR: program failed in Step deob"
581    goto BAD_EXIT
582endif
583
584## Step #1: Unifize the input T1
585echo "++ SSW Step 1"
586
587# [PT: Feb 11, 2019] now allows for option *not* to unifize
588# [PT: Feb 21, 2019] use Urad=30 when unifizing (better for final output)
589if ( ! -f $odir/anatU."${SubID}".nii ) then
590    if ( ${DO_UNIFIZE} ) then
591        3dUnifize                             \
592            -GM                               \
593            -clfrac 0.4                       \
594            -Urad 30                          \
595            -prefix $odir/anatU."$SubID".nii  \
596            -input "$Adataset"
597    else
598        echo "++ NOT unifizing dset."
599        3dcopy                                \
600            "$Adataset"                       \
601            $odir/anatU."$SubID".nii
602    endif
603endif
604
605if ( ${status} ) then
606    echo "** ERROR: program failed in Step U"
607    goto BAD_EXIT
608endif
609
610
611# [PT: Feb 11, 2019] now allows for option *not* to anismoothize
612if ( ! -f $odir/anatUA."${SubID}".nii ) then
613    if ( ${DO_ANISO} ) then
614        3danisosmooth                           \
615            -iters 1                            \
616            -3D                                 \
617            -automask                           \
618            -noneg                              \
619            -prefix $odir/anatUA."${SubID}".nii \
620            $odir/anatU."${SubID}".nii
621    else
622        echo "++ NOT anisosmoothing dset."
623        3dcopy                                  \
624            $odir/anatU."$SubID".nii            \
625            $odir/anatUA."$SubID".nii
626    endif
627endif
628
629if ( ${status} ) then
630    echo "** ERROR: program failed in Step UA"
631    goto BAD_EXIT
632endif
633
634# [PT: Jan 6, 2020] get rid of potentially *very* high vals: take
635# 98%ile of non-zero vox
636if ( ! -f $odir/anatUAC."${SubID}".nii ) then
637    if ( ${DO_CEIL} ) then
638        set vvv = `3dBrickStat                  \
639                    -percentile 98 1 98         \
640                    -non-zero                   \
641                    $odir/anatUA."${SubID}".nii`
642        # apply ceiling
643        3dcalc                                  \
644            -a $odir/anatUA."${SubID}".nii      \
645            -expr "maxbelow(${vvv[2]},a)"       \
646            -prefix $odir/anatUAC."${SubID}".nii
647
648    else
649        echo "++ NOT anisosmoothing dset."
650        3dcopy                                   \
651            $odir/anatUA."$SubID".nii            \
652            $odir/anatUAC."$SubID".nii
653    endif
654endif
655
656if ( ${status} ) then
657    echo "** ERROR: program failed in Step UAC"
658    goto BAD_EXIT
659endif
660
661
662# Step #2: Strip Skull (Ziad's way)
663echo "++ SSW Step 2"
664
665# [PT: Feb 11, 2019] now allows for option *not* to skullstrip
666if ( ! -f $odir/anatS."${SubID}".nii ) then
667    if ( ${DO_SKULLST} ) then
668        3dSkullStrip \
669            -input  $odir/anatUAC."$SubID".nii \
670            -prefix $odir/anatS."$SubID".nii   \
671            -debug 1 -ld 33 -niter 777         \
672            -shrink_fac_bot_lim 0.777          \
673            -exp_frac 0.0666 -orig_vol         \
674            $SSopt
675    else if ( "${mask_ss}" != "" ) then
676        echo "++ Instead of skullstripping, apply mask_ss dset."
677        3dcalc                                 \
678            -a $odir/anatUAC."$SubID".nii      \
679            -b ${mask_ss}                      \
680            -expr 'a*step(b)'                  \
681            -prefix $odir/anatS."$SubID".nii
682    else
683        echo "++ NOT pre-skullstripping dset."
684        3dcopy                                 \
685            $odir/anatUAC."$SubID".nii         \
686            $odir/anatS."$SubID".nii
687    endif
688endif
689
690if ( ${status} ) then
691    echo "** ERROR: program failed in Step 2"
692    goto BAD_EXIT
693endif
694
695## Step #3: run 3dQwarp first time to a moderate level (skull on)
696echo "++ SSW Step 3"
697
698if ( ${DO_GIANT} ) then
699    3dQwarp -echo_edu        ${saveall}                   \
700        $liteopt $verbopt                                 \
701        -base   "${Basedset}[1]"                          \
702        ${cost_nli}                                       \
703        -warpscale ${warpscale}                           \
704        -source $odir/anatUAC."$SubID".nii                \
705        -weight "${Basedset}[2]"                          \
706        -allineate -noneg -maxlev 5 -iwarp -awarp         \
707        -wtgaus ${gaus_wt}                                \
708        -inedge                                           \
709        -workhard:3:5 -nopenalty                          \
710        -prefix "${pref}TAL5.nii"                         \
711        -allopt '-twobest 11 -twopass -maxrot 45 -maxshf 40 -source_automask+2 -cmass'
712else
713    3dQwarp -echo_edu           ${saveall}                \
714        $liteopt $verbopt                                 \
715        -base   "${Basedset}[1]"                          \
716        ${cost_nli}                                       \
717        -warpscale ${warpscale}                           \
718        -source $odir/anatUAC."$SubID".nii                \
719        -weight "${Basedset}[2]"                          \
720        -allineate -noneg -maxlev 5 -iwarp -awarp         \
721        -wtgaus ${gaus_wt}                                \
722        -inedge                                           \
723        -workhard:3:5 -nopenalty                          \
724        -prefix "${pref}TAL5.nii"
725endif
726
727if ( ${status} ) then
728    echo "** ERROR: program failed in Step 3"
729    goto BAD_EXIT
730endif
731
732# -------- check early NL alignment to ref vol space ---------
733
734set ulay  = "${pref}TAL5.nii"
735set olay  = "${Basedset}"
736set opref_qc_nl0 = "${odir}/init_qc_01_nl0.${SubID}.jpg"
737
738@djunct_edgy_align_check               \
739    -montx            8                \
740    -monty            1                \
741    -ulay             "${ulay}"        \
742    -olay             "${olay}"        \
743    -box_focus_slices "${olay}"        \
744    -prefix           ${pref}MONT3
745
7462dcat                                  \
747    -gap 5                             \
748    -gap_col 255 204 153               \
749    -nx 1                              \
750    -ny 3                              \
751    -prefix ${opref_qc_nl0}            \
752    ${pref}MONT3*jpg
753
754## Step #4: mask off the skull using the template (second skull-strip)
755echo "++ SSW Step 4"
756
7573dmask_tool                        \
758    -input "${Basedset}[3]"        \
759    -dilate_input 2                \
760    -prefix "${pref}MASK.nii"
761
7623dcalc                             \
763    -a "${pref}MASK.nii"           \
764    -b "${pref}TAL5.nii"           \
765    -expr 'step(a)*b'              \
766    -prefix "${pref}TAL5mm.nii"
767
768if ( ${status} ) then
769    echo "** ERROR: program failed in Step 4"
770    goto BAD_EXIT
771endif
772
773## Step #5: warp this masked dataset back to original space
774echo "++ SSW Step 5"
775
7763dNwarpApply -echo_edu                 \
777    -nwarp "${pref}TAL5_WARPINV.nii"   \
778    -master $odir/anatS."$SubID".nii   \
779    -source "${pref}TAL5mm.nii"        \
780    -prefix "${pref}TAL5ww.nii"
781
782if ( ${status} ) then
783    echo "** ERROR: program failed in Step 5a"
784    goto BAD_EXIT
785endif
786
787## warp the mask itself (dilated) back to orig space
788
789\rm -f "${pref}MASK.nii"
7903dmask_tool  -echo_edu           \
791    -input "${Basedset}[3]"      \
792    -dilate_input 3              \
793    -prefix "${pref}MASK.nii"
7943dNwarpApply \
795    -nwarp "${pref}TAL5_WARPINV.nii"   \
796    -master $odir/anatS."$SubID".nii   \
797    -source "${pref}MASK.nii"          \
798    -prefix "${pref}MASKO.nii"         \
799    -ainterp NN
800
801if ( ${status} ) then
802    echo "** ERROR: program failed in Step 5b"
803    goto BAD_EXIT
804endif
805
806
807## merge these backward warped datasets with the 3dSkullStrip
808## output to get a better original skull-stripped result
809
8103dcalc -a "$odir/anatS.${SubID}.nii"  \
811       -b "${pref}TAL5ww.nii"         \
812       -c "${pref}MASKO.nii"          \
813       -expr 'step(c)*max(a,b)'       \
814       -prefix "$odir/anatSS.${SubID}.nii"
815
816if ( ${status} ) then
817    echo "** ERROR: program failed in Step 5c"
818    goto BAD_EXIT
819endif
820
821# DRG's erode-dilate trick for cleanup of little crap [16 Jan 2019]
8223dmask_tool                           \
823    -dilate_inputs -3 3               \
824    -prefix ${pref}de3.nii            \
825    -input  "$odir/anatSS.${SubID}.nii"
8263dcalc                                \
827    -a "$odir/anatSS.${SubID}.nii"    \
828    -b ${pref}de3.nii                 \
829    -expr 'a*step(b)'                 \
830    -prefix "$odir/anatSSc.${SubID}.nii"
831
832if ( ${status} ) then
833    echo "** ERROR: program failed in Step 5d"
834    goto BAD_EXIT
835endif
836
837# Throw in an automask step to clean up the outer edges [18 Jun 2019]
838\rm -f ${pref}de3.nii "$odir/anatSS.${SubID}.nii"
8393dAutomask                                            \
840    -apply_prefix "$odir/anatSSd.${SubID}.nii"        \
841    "$odir/anatSSc.${SubID}.nii"
842\mv -f "$odir/anatSSd.${SubID}.nii" "$odir/anatSS.${SubID}.nii"
843\rm -f "$odir/anatSSc.${SubID}.nii"
844
845if ( ${status} ) then
846    echo "** ERROR: program failed in Step 5e"
847    goto BAD_EXIT
848endif
849
850
851# don't do any more if skipping the final (precision) warp
852
853if ( $skipwarp ) goto CLEANUP
854
855## Step #6: affine transform that result to template space
856echo "++ SSW Step 6"
857
858# [PT: Jan 7, 2019] use lpa+ZZ cost func by default; also add in -twopass
859echo "++ SSW Step 6 ... piece 0"
8603dAllineate \
861    -1Dmatrix_apply "${pref}TAL5_Allin.aff12.1D"   \
862    -source $odir/anatSS."${SubID}".nii            \
863    -master "${pref}TAL5mm.nii"                    \
864    -cost ${cost_aff}                              \
865    -twopass                                       \
866    -autoweight -source_automask                   \
867    -final wsinc5                                  \
868    -prefix "${pref}AffSS.nii"
869
870if ( ${status} ) then
871    echo "** ERROR: program failed in Step 6.0"
872    goto BAD_EXIT
873endif
874
875# -------- check aff alignment to ref vol space ---------
876
877set ulay  = "${pref}AffSS.nii"
878set olay  = "${Basedset}"
879set opref_qc_aff = "${odir}/init_qc_02_aff.${SubID}.jpg"
880
881@djunct_edgy_align_check               \
882    -montx            8                \
883    -monty            1                \
884    -ulay             "${ulay}"        \
885    -olay             "${olay}"        \
886    -box_focus_slices "${olay}"        \
887    -prefix           ${pref}MONT6
888
8892dcat                                  \
890    -gap 5                             \
891    -gap_col 255 204 153               \
892    -nx 1                              \
893    -ny 3                              \
894    -prefix ${opref_qc_aff}            \
895    ${pref}MONT6*jpg
896
897
898
899## warp to template space (skull off),
900## initializing using the previous 3dQwarp -awarp output
901
902# Run 3dQwarp in several segments, to avoid gcc OpenMP bug
903#  where it freezes sometimes with inter-thread conflicts;
904#  this happens only in long runs, so breaking the run
905#  into pieces will (I hope) elide this annoyance - RWC :(
906
907# Piece number 1
908echo "++ SSW Step 6 ... piece 1"
9093dQwarp -echo_edu        ${saveall}      \
910    $liteopt $verbopt                    \
911    -base "${Basedset}[0]"               \
912    -source "${pref}AffSS.nii"           \
913    -iniwarp "${pref}TAL5_AWARP.nii"     \
914    -warpscale ${warpscale}              \
915    ${cost_nli}                          \
916    -inilev 1 -maxlev 5                  \
917    -wtgaus ${gaus_wt}                   \
918    -inedge                              \
919    -pblur -workhard:5:5                 \
920    -nodset -prefix "${pref}QQ5.nii"
921
922if ( ${status} ) then
923    echo "** ERROR: program failed in Step 6.1"
924    goto BAD_EXIT
925endif
926
927# Piece number 2
928echo "++ SSW Step 6 ... piece 2"
9293dQwarp -echo_edu      ${saveall}        \
930    $liteopt $verbopt                    \
931    -base "${Basedset}[0]"               \
932    -source "${pref}AffSS.nii"           \
933    -iniwarp "${pref}QQ5_WARP.nii"       \
934    -warpscale ${warpscale}              \
935    ${cost_nlf}                          \
936    -inilev 6 -maxlev 7 -workhard:6:7    \
937    -wtgaus ${gaus_wt}                   \
938    -inedge                              \
939    -pblur                               \
940    -nodset -prefix "${pref}QQ7.nii"
941
942if ( ${status} ) then
943    echo "** ERROR: program failed in Step 6.2"
944    goto BAD_EXIT
945endif
946
947if ( $minp > 13 ) then
948  # Final piece for coarse final patch size
949  echo "++ SSW Step 6 ... piece 3"
950  3dQwarp -echo_edu      ${saveall}        \
951      $liteopt $verbopt                    \
952      -base "${Basedset}[0]"               \
953      -source "${pref}AffSS.nii"           \
954      -iniwarp "${pref}QQ7_WARP.nii"       \
955      -warpscale ${warpscale}              \
956      ${cost_nlf}                          \
957      -inilev 8                            \
958      -wtgaus ${gaus_wt}                   \
959      -inedge                              \
960      -pblur -minpatch $minp               \
961      -Qfinal                              \
962      -prefix $odir/anatQQ."${SubID}".nii
963
964  if ( ${status} ) then
965      echo "** ERROR: program failed in Step 6.3a"
966      goto BAD_EXIT
967  endif
968
969else
970  # Penultimate piece for refined final patch size
971  echo "++ SSW Step 6 ... piece 3a"
972  @ mpp = $minp + 6
973  3dQwarp  -echo_edu       ${saveall}      \
974      $liteopt $verbopt                    \
975      -base "${Basedset}[0]"               \
976      -source "${pref}AffSS.nii"           \
977      -iniwarp "${pref}QQ7_WARP.nii"       \
978      -warpscale ${warpscale}              \
979      ${cost_nlf}                          \
980      -inilev 8 -minpatch $mpp             \
981      -wtgaus ${gaus_wt}                   \
982      -inedge                              \
983      -pblur                               \
984      -nodset -prefix "${pref}QQ9.nii"
985
986  if ( ${status} ) then
987      echo "** ERROR: program failed in Step 6.3b"
988      goto BAD_EXIT
989  endif
990
991  # Ultimate piece for refined final patch size
992  echo "++ SSW Step 6 ... piece 3b"
993  3dQwarp  -echo_edu       ${saveall}      \
994      $liteopt $verbopt                    \
995      -base "${Basedset}[0]"               \
996      -source "${pref}AffSS.nii"           \
997      -iniwarp "${pref}QQ9_WARP.nii"       \
998      -warpscale ${warpscale}              \
999      ${cost_nlf}                          \
1000      -inilev 10                           \
1001      -wtgaus ${gaus_wt}                   \
1002      -inedge                              \
1003      -pblur -minpatch $minp               \
1004      -prefix $odir/anatQQ."${SubID}".nii
1005
1006  if ( ${status} ) then
1007      echo "** ERROR: program failed in Step 6.3c"
1008      goto BAD_EXIT
1009  endif
1010
1011endif
1012
1013echo "++ SSW: done warping.  Finalize."
1014
1015# DRG's erode-dilate trick for cleanup of little stuff [16 Jan 2019]
10163dmask_tool                             \
1017    -dilate_inputs -1 1                 \
1018    -prefix ${pref}de3.nii              \
1019    -input  $odir/anatQQ."${SubID}".nii
10203dcalc                                  \
1021    -a $odir/anatQQ."${SubID}".nii      \
1022    -b ${pref}de3.nii                   \
1023    -expr 'a*step(b)'                   \
1024    -prefix $odir/anatQQc."${SubID}".nii
1025
1026if ( ${status} ) then
1027    echo "** ERROR: program failed in Finalize"
1028    goto BAD_EXIT
1029endif
1030
1031\rm -f ${pref}de3.nii anatQQ."${SubID}".nii
1032\mv -f $odir/anatQQc."${SubID}".nii $odir/anatQQ."${SubID}".nii
1033
1034## Step #7: Make two pretty pictures, and scram
1035echo "++ SSW Step 7"
1036
1037@snapshot_volreg $odir/anatQQ."${SubID}".nii $Basedset $odir/AM"${SubID}"
1038@snapshot_volreg $Basedset $odir/anatQQ."${SubID}".nii $odir/MA"${SubID}"
1039
1040JUMP_TO_EXTRA_QC:
1041
1042if ( ${DO_EXTRA_QC} ) then
1043
1044    # -------- check skullstripping in orig space ---------
1045
1046    set ulay  = "${odir}/anatU.${SubID}.nii"
1047    set olay0 = "${odir}/anatSS.${SubID}.nii"
1048    set olay  = "${pref}mask.nii"
1049    set opref_exqc1 = "${odir}/QC_anatSS.${SubID}.jpg"
1050
1051    3dcalc                                 \
1052        -overwrite                         \
1053        -a ${olay0}                        \
1054        -expr 'bool(a)'                    \
1055        -prefix ${olay}                    \
1056        -byte
1057
1058    @chauffeur_afni                        \
1059        -ulay ${ulay}                      \
1060        -ulay_range "2%" "98%"             \
1061        -olay ${olay}                      \
1062        -func_range_perc_nz 1              \
1063        -set_subbricks 0 0 0               \
1064        -box_focus_slices ${olay}          \
1065        -cbar "Reds_and_Blues_Inv"         \
1066        -opacity 4                         \
1067        -prefix  ${pref}MONT1              \
1068        -save_ftype JPEG                   \
1069        -montx 9 -monty 1                  \
1070        -montgap 3                         \
1071        -set_xhairs OFF                    \
1072        -label_mode 1 -label_size 4
1073
1074    2dcat                                  \
1075        -gap 5                             \
1076        -gap_col 66 184 254                \
1077        -nx 1                              \
1078        -ny 3                              \
1079        -prefix ${opref_exqc1}             \
1080        ${pref}MONT1*jpg
1081
1082    # -------- check alignment in ref vol space ---------
1083
1084    set ulay  = "${odir}/anatQQ.${SubID}.nii"
1085    set olay  = "${Basedset}"
1086    set opref_exqc2 = "${odir}/QC_anatQQ.${SubID}.jpg"
1087
1088    @djunct_edgy_align_check               \
1089        -montx            8                \
1090        -monty            1                \
1091        -ulay             "${ulay}"        \
1092        -olay             "${olay}"        \
1093        -box_focus_slices "${olay}"        \
1094        -prefix           ${pref}MONT2
1095
1096    2dcat                                  \
1097        -gap 5                             \
1098        -gap_col 66 184 254                \
1099        -nx 1                              \
1100        -ny 3                              \
1101        -prefix ${opref_exqc2}             \
1102        ${pref}MONT2*jpg
1103
1104endif
1105
1106# ------------------------------------------------------------------
1107
1108## Clean up the junk
1109
1110CLEANUP:
1111
1112## Rename affine warp matrix
1113echo "++ SSW cleanup"
1114
1115# [PT: Feb 19, 2020] Put the mv cmd in an IF condition, because with
1116# re-running extra QC images, might not have the first file name;
1117# now, avoid an unnecessary error msg.
1118if ( -e "${pref}TAL5_Allin.aff12.1D" ) then
1119    \mv -f "${pref}TAL5_Allin.aff12.1D" $odir/anatQQ."${SubID}".aff12.1D
1120endif
1121
1122if ( $doclean ) then
1123  \rm -f "${pref}"*
1124endif
1125
1126echo "--------------------------------------------------------"
1127echo "++ DONE with SSW.  Check QC images:"
1128printf "   $odir/AM${SubID}.jpg  $odir/MA${SubID}.jpg "
1129if ( ${DO_EXTRA_QC} ) then
1130printf " ${opref_exqc1}  ${opref_exqc2}"
1131endif
1132echo ""
1133echo ""
1134echo "   In ref/base space (check warping):"
1135echo "     ulay=anat, olay=base edges  :  $odir/AM${SubID}.jpg"
1136echo "     ulay=base edges, olay=anat  :  $odir/MA${SubID}.jpg"
1137if ( ${DO_EXTRA_QC} ) then
1138    echo "     ulay=anat, olay=base edges  :  ${opref_exqc2}"
1139    echo ""
1140    echo "   In native space (check skullstripping):"
1141    echo "     ulay=anat, olay=SS mask     :  ${opref_exqc1}"
1142endif
1143echo ""
1144echo "--------------------------------------------------------"
1145
1146goto GOOD_EXIT
1147
1148# ========================================================================
1149# ========================================================================
1150
1151SHOW_HELP:
1152
1153cat <<EOF
1154
1155OVERVIEW ~1~
1156
1157This script has dual purposes for processing a given subject's
1158anatomical volume:
1159    + to skull-strip the brain, and
1160    + to calculate the warp to a reference template/standard space.
1161Automatic snapshots of the registration are created, as well, to help
1162the QC process.
1163
1164This program cordially ties in directly with afni_proc.py, so you can
1165run it beforehand, check the results, and then provide both the
1166skull-stripped volume and the warps to the processing program.  That
1167is convenient!
1168
1169Current version = $version
1170Authorship      = Bob, Bob, there is one Bob, He spells it B-O-B.
1171
1172# -----------------------------------------------------------------
1173
1174USAGE ~1~
1175
1176    ${this_prog}             \
1177        -input  AA           \
1178        -base   BB           \
1179        -subid  SS           \
1180       {-odir   OD}          \
1181       {-minp   MP}          \
1182       {-nolite}             \
1183       {-skipwarp}           \
1184       {-unifize_off}        \
1185       {-init_skullstr_off}  \
1186       {-extra_qc_off}       \
1187       {-jump_to_extra_qc}   \
1188       {-cost_aff CA}        \
1189       {-cost_nl_init CNI}   \
1190       {-cost_nl_final CNF}  \
1191       {-deoblique}          \
1192       {-deoblique_refitly}  \
1193       {-warpscale WS}       \
1194       {-SSopt 'strings'     \
1195       {-aniso_off}          \
1196       {-ceil_off}           \
1197       {-tmp_name_nice}      \
1198       {-echo}               \
1199       {-verb}               \
1200       {-noclean}
1201
1202where (note: many of the options with 'no' and 'off' in their name are
1203really just included for backwards compatibility, as this program has
1204grown/improved over time):
1205
1206  -input  AA :(req) an anatomical dataset, *not* skull-stripped, with
1207              resolution about 1 mm.
1208
1209  -base   BB :(req) a base template dataset, with contrast similar to
1210              the input AA dset, probably from some kind of standard
1211              template.
1212              NB: this dataset is not *just* a standard template,
1213              because it is not a single volume-- read about its
1214              composition in the NOTES on the 'The Template Dataset',
1215              below.
1216              The program first checks if the dset BB exists as
1217              specified; if not, then if just the filename has been
1218              provided it searches the AFNI_GLOBAL_SESSION,
1219              AFNI_PLUGINPATH, and afni bin directory (in that order)
1220              for the named dataset.
1221
1222  -subid  SS :(req) name code for output datasets (e.g., 'sub007').
1223
1224  -odir   OD :(opt) output directory for all files from this program
1225              (def: directory of the '-input AA').
1226
1227  -minp   MP :(opt) minimum patch size on final 3dQwarp (def: 11).
1228
1229  -nolite    :(opt) Do not use the '-lite' option with 3dQwarp;
1230              This option is used for backward compatibility, if you want
1231              to run 3dQwarp the same way as older versions of @SSwarper.
1232              The new way (starting Jan 2019) is to use the '-lite'
1233              option with 3dQwarp to speed up the calculations.
1234              (def: use '-lite' for faster calculations).
1235
1236  -skipwarp  :(opt) Do not compute past the output of anatSS.{subid}.nii.
1237              This option is used if you just want the skull-stripped
1238              result in original coordinates, without the warping
1239              to the template space (anatQQ). The script will run faster.
1240
1241  -deoblique :(opt) apply obliquity information to deoblique the input
1242              volume ('3dWarp -deoblique -wsinc5 ...'), as an initial step.
1243              This might introduce the need to overcome a large rotation
1244              during the alignment, though!
1245
1246  -deoblique_refitly :(opt) purge obliquity information to deoblique
1247              the input volume (copy, and then '3drefit -deoblique ...'), 
1248              as an initial step.  This might help when data sets are
1249              very... oblique.
1250
1251  -warpscale WS :(opt) opt to control flexibility of warps in 3dQwarp and
1252              how they adjust with patch size;  see 3dQwarp's help for 
1253              more info. Allowed values of WS are in range [0.1, 1.0].
1254              (def: 1.0)
1255
1256  -giant_move :(opt) when starting the initial alignment to the template,
1257              apply the same parameter expansions to 3dAllineate that
1258              align_epi_anat.py does with the same option flag.  This
1259              might be useful if the brain has a very large angle away
1260              from "typical" ones, etc.
1261
1262  -unifize_off :(opt) don't start with a 3dUnifize command to try reduce
1263              effects of brightness inhomogeneities.  Probably only
1264              useful if unifizing has been previously performed on the
1265              input dset.
1266
1267  -aniso_off :(opt) don't preprocess with a 3danisosmooth command to
1268              try reduce effects of weird things (in a technical
1269              sense).  Possible that this will never be used in the
1270              history of running this program.
1271
1272  -ceil_off  :(opt) by default, after anisosmoothing, this program
1273              will apply put a ceiling on values in the dset, to get rid
1274              of possible outliers (ceil = 98%ile of non-zero voxels in
1275              the whole volume).  This option will turn that off.
1276
1277  -init_skullstr_off :(opt) don't preprocess with a 3dSkullstrip command
1278              to roughly isolated brain in the beginning.  This might
1279              be useful with macaque dsets.
1280
1281  -extra_qc_off :(opt) don't make extra QC images QC*jpg (for some
1282              unknown reason).
1283
1284  -jump_to_extra_qc :(opt) just make the two QC*jpg images from a
1285              previous run of @SSwarper.  These QC*jpg images are new
1286              QC output (as of late Feb, 2020), so this might be
1287              useful to add a quick check to previously run data.
1288              This command would just be tacked on to previously
1289              executed one.
1290
1291  -cost_aff CA :(opt) specify cost function for affine (3dAllineate)
1292              part of alignment.  Here, 'CA' would be just the name of
1293              the cost function to be provided after '-cost ..' (def:
1294              is now "lpa+ZZ").  This is probably only here for
1295              backwards compatability to older @SSwarper (where def
1296              was 'hel').
1297
1298  -cost_nl_init CNI 
1299             :(opt) specify cost function for initial nonlinear
1300              (3dQwarp) part of alignment.  Here, 'CNI' would be the
1301              cost function name to be provided (def: is now "lpa").
1302              This is probably only here for backwards compatability
1303              to older @SSwarper (where def was 'pcl').
1304
1305  -cost_nl_final CNF 
1306             :(opt) specify cost function for final nonlinear
1307              (3dQwarp) parts of alignment.  Here, 'CNF' would be the
1308              cost function to be provided (def: is now "pcl").  This
1309              is separate from the initial nonlinear warp cost values
1310              '-cost_nl_init ..', because using those here might be
1311              pretty slow; however, using "lpa" here might help
1312              results.
1313
1314  -SSopt 'strings' :(opt) The content of 'strings' (which should be
1315              in quotes if there are any blanks) is copied to the
1316              end of the 3dSkullStrip command line. Example:
1317                -SSopt '-o_ply Fred.Is.Wonderful'
1318              to have 3dSkullStrip produce a .ply surface file
1319              as an additional output.
1320
1321  -mask_ss MSS :(opt) as an alternative to skullstripping at an early
1322              stage, you can provide a mask to be used before the
1323              initial affine alignment.  The mask MSS can come from
1324              anywhere, but @SUMA_Make_Spec_FS now makes a convenient
1325              one from the FS parcellation (though it would have to be
1326              resampled to the input anatomical's grid).
1327
1328  -tmp_name_nice :(opt) default temporary "junk.*" filenames include
1329              a large, random char string.  This is ugly, but useful
1330              if outputting several different SSW runs into the same
1331              directory that we intermediate files (very likely) don't
1332              get overwritten.  However, if you prefer, you can use a
1333              nicer, non-random intermediate file prefix: "junk_ssw".
1334              I would use this when the output dir ("-odir ..")
1335              doesn't contain multiple SSW outputs.
1336
1337  -verb      :(opt) Apply the '-verb' option to 3dQwarp, to get more
1338              verbose progress information - mostly used for debugging.
1339
1340  -echo      :(opt) Run the script with "set echo", for extra verbosity
1341              in the terminal output.  Mainly for debugging times. 
1342
1343  -noclean   :(opt) Do not delete the 'junk' files at the end of
1344              computations - mostly used for debugging and testing.
1345
1346# -----------------------------------------------------------------
1347
1348REFERENCE DATASETS ~1~
1349
1350If you are reading this message, then several reference data sets
1351(base volumes) for @SSwarper now exist within the AFNI realm. Oh, what
1352a time it is to be alive.  A current list includes:
1353
1354+ MNI152_2009_template_SSW.nii.gz
1355+ TT_N27_SSW.nii.gz
1356+ HaskinsPeds_NL_template1.0_SSW.nii.gz
1357
1358Some of these are distributed with the AFNI binaries, and other may be
1359found online. You can make other reference base templates in whatever
1360space you prefer, but note that it must have several subvolumes of
1361information included-- see NOTES on the 'The Template Dataset', below
1362(which also contains a link to the @SSwarper template tutorial online
1363help).
1364
1365# ----------------------------------------------------------------------
1366
1367OUTPUTS ~1~
1368
1369Datasets ~2~
1370
1371  Suppose the -prefix is 'sub007' (because you scanned Bond, JamesBond?).
1372  Then the outputs from this script will be"
1373
1374  anatDO.sub007.nii       = deobliqued version of original dataset;
1375                            (*only if* using '-deoblique' opt);
1376  anatU.sub007.nii        = intensity uniform-ized original dataset
1377                            (or, if '-unifize_off' used, a copy of orig dset);
1378  anatUA.sub007.nii       = anisotropically smoothed version of the above
1379                            (or, if '-aniso_off' used, a copy of anatU.*.nii)
1380  anatUAC.sub007.nii      = ceiling-capped ver of the above (at 98%ile of 
1381                            non-zero values)
1382                            (or, if '-ceil_off' used, a copy of anatUA.*.nii)
1383
1384  anatS.sub007.nii        = first pass skull-stripped original dataset
1385                            (or, if '-init_skullstr_off' used, a copy of
1386                            anatUAC.*.nii);
1387  anatSS.sub007.nii       = second pass skull-stripped original dataset;
1388                            * note that anatS and anatSS are 'original'
1389                              in the sense that they are aligned with
1390                              the input dataset - however, they have been
1391                              unifized and weakly smoothed: they are
1392                              stripped versions of anatUAC; if you want
1393                              a skull-stripped copy of the input with
1394                              no other processing, use a command like
1395                                3dcalc -a INPUTDATASET        \
1396                                       -b anatSS.sub007.nii   \
1397                                       -expr 'a*step(b)'      \
1398                                       -prefix anatSSorig.sub007.nii
1399
1400  anatQQ.sub007.nii       = skull-stripped dataset nonlinearly warped to
1401                            the base template space;
1402  anatQQ.sub007.aff12.1D  = affine matrix to transform original dataset
1403                            to base template space;
1404  anatQQ.sub007_WARP.nii  = incremental warp from affine transformation
1405                            to nonlinearly aligned dataset;
1406
1407  * The .aff12.1D and _WARP.nii transformations need to be catenated to get
1408    the full warp from orginal space to the base space; example:
1409
1410      3dNwarpApply -nwarp 'anatQQ.sub007_WARP.nii anatQQ.sub007.aff12.1D' ...
1411
1412QC images ~2~
1413
1414  AMsub007.jpg            = 3x3 snapshot image of the anatQQ.sub007.nii
1415                            dataset with the edges from the base template
1416                            overlaid -- to check the alignment;
1417  MAsub007.jpg            = similar to the above, with the roles of the
1418                            template and the anatomical datasets reversed.
1419  QC_anatQQ.sub007.jpg    = like AM*.jpg, but 3 rows of 8 slices
1420  QC_anatSS.sub007.jpg    = check skullstripping in orig space: ulay is
1421                            input dset, and olay is mask of
1422                            skullstripped output (anatSS* dset)
1423
1424  init_qc_00_overlap_uinp_obase.jpg
1425    o [ulay] original source dset
1426      [olay] original base dset
1427    o single image montage to check initial overlap of source and base,
1428      ignoring any obliquity that might be present (i.e., the way AFNI 
1429      GUI does by default, and also how alignment starts)
1430    o if initial overlap is not strong, alignment can fail or
1431      produce weirdness
1432    o *if* either dset has obliquity, then an image of both after 
1433      deobliquing with 3dWarp is created (*DEOB.jpg), and a text file
1434      about obliquity is also created (*DEOB.txt).
1435
1436* It is important to examine (at least) the two .jpg snapshot images to
1437  make sure that the skull-stripping and nonlinear warping worked well.
1438
1439USING SSW WITH AFNI_PROC.PY ~1~
1440
1441  When B-O-B uses @SSwarper for skull-stripping plus warping, He
1442  gives afni_proc.py these options (among others), after running
1443  ${this_prog} successfully -- here, 'subj' is the subject
1444  identifier:
1445
1446  |  set btemplate = MNI152_2009_template_SSW.nii.gz
1447  |  set tpath = \`@FindAfniDsetPath \${btemplate}\`
1448  |  if ( "$tpath" == "" ) exit 1
1449  |
1450  |  afni_proc.py                                                  \
1451  |    [...other stuff here: processing blocks, options...]        \
1452  |    -copy_anat anatSS.\${subj}.nii                               \
1453  |    -anat_has_skull no                                          \
1454  |    -align_opts_aea -ginormous_move -deoblique on -cost lpc+ZZ  \
1455  |    -volreg_align_to MIN_OUTLIER                                \
1456  |    -volreg_align_e2a                                           \
1457  |    -volreg_tlrc_warp -tlrc_base $tpath/$btemplate              \
1458  |    -tlrc_NL_warp                                               \
1459  |    -tlrc_NL_warped_dsets                                       \
1460  |       anatQQ.\${subj}.nii                                       \
1461  |       anatQQ.\${subj}.aff12.1D                                  \
1462  |       anatQQ.\${subj}_WARP.nii
1463
1464
1465NOTES ~1~
1466
1467The Template dataset ~2~
1468
1469  Any reference base template dataset, such as
1470  MNI152_2009_template_SSW.nii.gz, must have the first *4* volumes here
1471  (and can have the optional 5th for later uses, as described):
1472    [0] = skull-stripped template brain volume
1473    [1] = skull-on template brain volume
1474    [2] = weight mask for nonlinear registration, with the
1475          brain given greater weight than the skull
1476    [3] = binary mask for the brain
1477    [4] = binary mask for gray matter plus some CSF (slightly dilated)
1478          ++ this volume is not used in this script
1479          ++ it is intended for use in restricting FMRI analyses
1480             to the 'interesting' parts of the brain
1481          ++ this mask should be resampled to your EPI spatial
1482             resolution (see program 3dfractionize), and then
1483             combined with a mask from your experiment reflecting
1484             your EPI brain coverage (see program 3dmask_tool).
1485
1486  More information about making these (with scripts) is provided on
1487  the Interweb:
1488    https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/template_atlas/sswarper_base.html
1489
1490The steps being run ~2~
1491
1492  You Know My Methods, Watson...
1493
1494  #1: Uniform-ize the input dataset's intensity via 3dUnifize.
1495       ==> anatU.sub007.nii
1496  #2: Strip the skull with 3dSkullStrip, with mildly agressive settings.
1497       ==> anatS.sub007.nii
1498  #3: Nonlinearly warp (3dQwarp) the result from #1 to the skull-on
1499      template, driving the warping to a medium level of refinement.
1500  #4: Use a slightly dilated brain mask from the template to
1501      crop off the non-brain tissue resulting from #3 (3dcalc).
1502  #5: Warp the output of #4 back to original anatomical space,
1503      along with the template brain mask, and combine those
1504      with the output of #2 to get a better skull-stripped
1505      result in original space (3dNwarpApply and 3dcalc).
1506       ==> anatSS.sub007.nii
1507  #6  Restart the nonlinear warping, registering the output
1508      of #5 to the skull-off template brain volume (3dQwarp).
1509       ==> anatQQ.sub007.nii (et cetera)
1510  #7  Use @snapshot_volreg to make the pretty pictures.
1511       ==> AMsub007.jpg and MAsub007.jpg
1512
1513Temporary files ~2~
1514
1515  If the script crashes for some reason, it might leave behind files
1516  whose names start with 'junk' -- you should delete these files
1517  manually.
1518
1519WHAT TO DO IF RESULTS ARE WAY OFF? ~1~
1520
1521  The importance of initial dset overlap ~2~
1522
1523  Always, always, always check the initial image made by SSW when it
1524  runs:  
1525
1526      init_qc_00_overlap_uinp_obase.jpg
1527
1528  This image tells you how well your datasets overlap initially before
1529  the alignment work begins. **The better the overlap, the lower the
1530  chance that something weird happens in your output.** All the SSW
1531  templates have reasonable coordinates, meaning that (x, y, z) = (0,
1532  0, 0) is in a good spot for it.  If there is poor overlap, probably
1533  your input dataset has weird/bad coordinates for some reason.  
1534
1535  You can use @Align_Centers to put your anatomical dset in a better
1536  spot (though note, if you are going to be processing EPI data
1537  afterwards, you will want to move that along, as well, perhaps as a
1538  "child" dataset).
1539
1540  By far the most common problem leading to obviously bad outputs is
1541  that the initial datasets are waaay far apart when they start, and
1542  the program gets stuck in a false minimum of solutions.
1543
1544  Other issues ~2~
1545
1546  Sometimes, it can be hard to separate the brain from dura and/or
1547  skull surrounding the brain.  If little bits are left around in the
1548  masking images, then perhaps adding one of the following options for
1549  will help (this can help the initial skullstripping):
1550
1551    -SSopt '-blur_fwhm 2'
1552    -SSopt '-blur_fwhm 3'
1553
1554  Any other questions/oddities, please don't hesitate to inquire on
1555  the AFNI Message Board!
1556
1557EXAMPLES ~1~
1558
1559  1) Run the program, deciding what the main output directory will be
1560  called (e.g., based on the subject ID):
1561
1562    @SSwarper                                    \
1563        -input  anat_t1w.nii.gz                  \
1564        -base   MNI152_2009_template_SSW.nii.gz  \
1565        -subid  sub-001                          \
1566        -odir   group/o.aw_sub-001
1567
1568  2) Same as above, but since we are using one outdir per subject, use
1569  more aesthetically pleasing names of temporary files (which get
1570  deleted, anyways):
1571
1572    @SSwarper                                    \
1573        -tmp_name_nice                           \
1574        -input  anat_t1w.nii.gz                  \
1575        -base   MNI152_2009_template_SSW.nii.gz  \
1576        -subid  sub-001                          \
1577        -odir   group/o.aw_sub-001
1578
1579  3) As of version 2.5, you can input a mask to be used instead of
1580  skullstripping.  For example, a good one might be the
1581  parcellation-derived (but filled in) mask from @SUMA_Make_Spec_FS
1582  after running FS's recon-all (though you will have to resample it
1583  from the FS output grid to that of your input anatomical):
1584
1585    @SSwarper                                     \
1586        -tmp_name_nice                            \
1587        -input   anat_t1w.nii.gz                  \
1588        -mask_ss fs_parc_wb_mask_RES.nii.gz       \
1589        -base    MNI152_2009_template_SSW.nii.gz  \
1590        -subid   sub-001                          \
1591        -odir    group/o.aw_sub-001
1592
1593EOF
1594
1595   goto GOOD_EXIT
1596
1597SHOW_VERSION:
1598   echo "version  $version (${rev_dat})"
1599   goto GOOD_EXIT
1600
1601FAIL_MISSING_ARG:
1602    echo "** ERROR! Missing an argument after option flag: '$argv[$ac]'"
1603    goto BAD_EXIT
1604
1605BAD_EXIT:
1606    exit 1
1607
1608GOOD_EXIT:
1609    exit 0
1610