1#!/bin/tcsh
2# @ROI_modal_grow
3# dilate ROIs in volumetric dataset with non-zero modal smoothing
4
5
6@global_parse `basename $0` "$*" ; if ($status) exit 0
7
8set progname = @ROI_modal_grow
9
10# --------------------- version history with changes -----------------------
11#
12set version = "1.00"
13#
14
15set dset = ""
16set outdir = "rmgrow"
17set niters = ""
18set maskset = ""
19set outprefix = "rmg"
20set keep_rm_files = ""
21set NN = "-1"
22
23set here = `pwd`
24# ------------------------ process user options --------------------------
25
26if ( $#argv == 0 ) goto HELP
27
28set ac = 1
29while ($ac <= $#argv)
30    if ("$argv[$ac]" == "-help" || "$argv[$ac]" == "-h") then
31        goto HELP
32
33    else if ("$argv[$ac]" == "-ver") then
34        echo $version
35        exit 0
36
37    # -------------------
38
39    else if ("$argv[$ac]" == "-input") then
40        set this_opt = "$argv[$ac]"
41        @ ac ++
42        if ( $ac > $#argv ) then
43            echo "** missing parameter for option '${this_opt}'"
44            exit 1
45        endif
46        set dset =  $argv[$ac]
47    else if ("$argv[$ac]" == "-mask") then
48        set this_opt = "$argv[$ac]"
49        @ ac ++
50        if ( $ac > $#argv ) then
51            echo "** missing parameter for option '${this_opt}'"
52            exit 1
53        endif
54        set maskset =  $argv[$ac]
55    else if ("$argv[$ac]" == "-outdir") then
56        set this_opt = "$argv[$ac]"
57        @ ac ++
58        if ( $ac > $#argv ) then
59            echo "** missing parameter for option '${this_opt}'"
60            exit 1
61        endif
62        set outdir =  $argv[$ac]
63    else if ("$argv[$ac]" == "-niters") then
64        set this_opt = "$argv[$ac]"
65        @ ac ++
66        if ( $ac > $#argv ) then
67            echo "** missing parameter for option '${this_opt}'"
68            exit 1
69        endif
70        set niters =  $argv[$ac]
71    else if ("$argv[$ac]" == "-prefix") then
72        set this_opt = "$argv[$ac]"
73        @ ac ++
74        if ( $ac > $#argv ) then
75            echo "** missing parameter for option '${this_opt}'"
76            exit 1
77        endif
78        set outprefix =  $argv[$ac]
79
80    else if ("$argv[$ac]" == "-NN") then
81        set this_opt = "$argv[$ac]"
82        @ ac ++
83        if ( $ac > $#argv ) then
84            echo "** missing parameter for option '${this_opt}'"
85            exit 1
86        endif
87        if ("$argv[$ac]" == "1") then
88           # facing neighbors only - default
89           set NN = "-1"
90        else if ("$argv[$ac]" == "2") then
91           # >sqrt 2 - 2D diagonal neighbors
92           set NN = "-1.415"
93        else if ("$argv[$ac]" == "3") then
94           # >sqrt 3 - 3D diagonal neighbors
95           set NN = "-1.74"
96        else
97           echo "Only values of 1,2,3 allowed"
98           exit 1
99        endif
100    endif
101
102    else if ( "$argv[$ac]" == "-keep_rm_files" ) then
103        set keep_rm_files = 1
104
105    else if ( "$argv[$ac]" == "-echo" ) then
106        set echo
107
108    # ---------- fin ------------------
109
110    else
111        echo "** unknown option $argv[$ac]"
112        exit 1
113    endif
114    @ ac ++
115end
116
117# ------------------------- check for required inputs
118if ($dset == "") then
119    echo "No input dataset provided"
120    exit 1
121endif
122
123if ("${niters}" == "") then
124    echo "Must select number of iterations"
125    exit 1
126endif
127
128# ------------------------- make OUTDIR
129# We don't allow output directory to be $PWD -- must be a subdir; too
130# complicated with copying fnames
131
132\mkdir -p "${outdir}"
133cd "${outdir}"
134if ( "${PWD}" == "${here}" ) then
135    echo "** ERROR: Output directory cannot be same current working dir."
136    echo "          Please specify new output directory"
137    exit 1
138endif
139cd -
140
141# copy the datasets into the directory and rename
1423dcopy $dset "${outdir}/rm_rgm_00.nii.gz"
143# copy any atlas and label tables over
1443drefit -copytables $dset "${outdir}/rm_rgm_00.nii.gz"
145
146if ("$maskset" != "") then
147   3dcopy $maskset "${outdir}/maskset.nii.gz"
148   set maskset = maskset.nii.gz
149endif
150
151cd "${outdir}"
152set inset = rm_rgm_00.nii.gz
153
154foreach iter (`count -digits 2 1 $niters`)
155    3dLocalstat -nbhd "SPHERE(${NN})" -stat nzmode -overwrite \
156       -prefix rm_rgm_${iter}.nii.gz $inset
157    set inset = rm_rgm_${iter}.nii.gz
158end
159
160if ("$maskset" != "") then
161   3dcalc -datum short -nscale -a $inset -b $maskset  \
162      -expr 'a*step(b)' -prefix ${outprefix}.nii.gz -overwrite
163else
164   3dcalc -datum short -nscale -a $inset -expr a  \
165      -prefix ${outprefix}.nii.gz -overwrite
166endif
167
168# copy any atlas and label tables over to final output
1693drefit -cmap INT_CMAP -copytables rm_rgm_00.nii.gz ${outprefix}.nii.gz
170
171
172if ($keep_rm_files != "") then
173   \rm rm_rgm_${iter}.nii.gz
174endif
175
176exit 0
177
178
179# ===========================================================================
180HELP:
181
182cat << SCRIPT_HELP_STRING
183
184Overview ~1~
185
186Script to grow a set of regions in a volumetric dataset using modal
187smoothing.
188
189Usage Example ~1~
190
191    @ROI_modal_grow                                      \
192      -input  mydset.nii.gz                              \
193      -outdir myoutdir                                   \
194      -niters 5                                          \
195      -mask mymask.nii.gz                                \
196      -prefix rmg_name                                  
197
198   Note only the input dataset and the number of iteration levels
199   are required.
200
201Options ~1~
202
203    -input input_dset   :required input dataset. This dataset should be
204                         set of integer values. The program mostly assumes
205                         approximate isotropic voxels.
206
207    -outdir outdirname  :directory name for output. All output goes to
208                         this directory. Default is rmgrow.
209
210    -niters nn          :number of iterations for modal growth - 
211                         something like dilation level here - generally
212                         this will make sense for values from about 1-10
213
214    -mask maskset       :mask dataset at same grid as the input dataset.
215                         this could be a dilated version of the original mask
216                         or a larger region like a cortical ribbon mask.
217                         Not required but often desirable.
218
219    -prefix baseprefix  :base name of final output dataset, i.e. baseprefix.nii.gz
220                         Default is rmg, so output would be rmg.nii.gz
221
222    -NN [1,2,3]         :neighborhood type using in finding mode, 
223                         1 - facing neighbors, 2-edges, 3-corners
224
225Also see these programs with these similar functions: 
226   3dROImaker - grows regions using regular dilation iteratively
227   3dmask_tool - dilates and erodes sets of input files or a single output mask
228   3dmerge, 3dAutomask, 3dcalc - all can do dilation
229   3dLocalstat - modal and nonzero modal smoothing
230   ROIgrow - dilates surface ROI (patches) within the mesh of the surface
231
232SCRIPT_HELP_STRING
233
234   exit 0
235