1 #include "SUMA_suma.h"
2
usage_SurfExtrema(SUMA_GENERIC_ARGV_PARSE * ps)3 void usage_SurfExtrema (SUMA_GENERIC_ARGV_PARSE *ps)
4 {
5 static char FuncName[]={"usage_SurfExtrema"};
6 char * s = NULL, *sio=NULL, *st = NULL, *sts = NULL;
7 int i;
8 s = SUMA_help_basics();
9 sio = SUMA_help_IO_Args(ps);
10 printf ( "\n"
11 "Usage: A program finding the local extrema in a dataset.\n"
12 " The program finds the nodes with the highest value within an Rmm\n"
13 " radius, and at which the gradient of the signal meets a preset\n"
14 " threshold.\n"
15 " By default, the program searches for maxima.\n"
16 "\n"
17 " -input DSET = Dset in which Extrema are to be identified.\n"
18 " If you do not specify one, the program use the surface's\n"
19 " convexity dataset.\n"
20 " -hood R = Neighborhood of node n consists of nodes within R \n"
21 " -nbhd_rad R = distance from n as measured by the shortest \n"
22 " distance along the mesh.\n"
23 " Default is 8 mm\n"
24 " -thresh TH = Do not consider nodes with value less than TH\n"
25 " Default is 0\n"
26 " -gthresh GTH = Do not consider nodes with gradient less than GTH.\n"
27 " Default is 0.01\n"
28 " -gscale SCL = What scaling to apply to gradient computation.\n"
29 " Choose from:\n"
30 " NONE: g[n] = sum(v[n]-v[k])/Nk with k the neighbs. of n\n"
31 " LMEAN : Divide g[n] by mean of n and its neighbors * 100\n"
32 " GMEAN : Divide g[n] by mean of all nodes in mask * 100\n"
33 " Default is LMEAN\n"
34 " -extype TYP = Find maxima, minima, or extrema.\n"
35 " TYP is one of: MAX (default)\n"
36 " MIN \n"
37 " ABS \n"
38 " -prefix PREFIX = Prefix of two output data sets.\n"
39 " First dset is called PREFIX.grd and contains the \n"
40 " scaled average gradient values.\n"
41 " Second dset is called PREFIX.ext and contains the \n"
42 " nodes with maximum values. The value of a non-zero\n"
43 " node is its rank.\n"
44 " -table TABLE = Name of file in which to store a record of the extrema\n"
45 " found. The header part of TABLE contains examples\n"
46 " for easily extracting certain values from it.\n"
47 "\n"
48 " Examples:\n"
49 " ---------\n"
50 " 1- SurfExtrema -i SUMA/std141.rh.smoothwm.asc \\\n"
51 " -input pb05.rh.niml.dset'[1]' \\\n"
52 " -gscale LMEAN \\\n"
53 " -prefix ex1.rh \\\n"
54 " -table ex1.log \n"
55 " \n"
56 "%s"
57 "%s"
58 "\n", sio, s);
59 SUMA_free(s); s = NULL; SUMA_free(st); st = NULL;
60 SUMA_free(sio); sio = NULL;
61 printf(" Ziad S. Saad SSCC/NIMH/NIH saadz@mail.nih.gov \n");
62 exit(0);
63 }
64
65
66
SUMA_SurfExtrema_ParseInput(char * argv[],int argc,SUMA_GENERIC_ARGV_PARSE * ps)67 SUMA_GENERIC_PROG_OPTIONS_STRUCT *SUMA_SurfExtrema_ParseInput(
68 char *argv[], int argc, SUMA_GENERIC_ARGV_PARSE *ps)
69 {
70 static char FuncName[]={"SUMA_SurfExtrema_ParseInput"};
71 SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt=NULL;
72 int kar;
73 SUMA_Boolean brk;
74 SUMA_Boolean LocalHead = NOPE;
75
76 SUMA_ENTRY;
77
78 Opt = SUMA_Alloc_Generic_Prog_Options_Struct();
79 Opt->ps = ps; /* for convenience */
80 Opt->NodeDbg = -1;
81 Opt->out_prefix = NULL;
82 Opt->r = 8.0;
83 Opt->t = 0.0;
84 Opt->t2 = 0.01;
85 Opt->b1 = SUMA_MEAN_GRAD_SCALE;
86 Opt->b2 = SUMA_MAXIMUS;
87 Opt->s = NULL;
88 kar = 1;
89 brk = NOPE;
90 while (kar < argc) { /* loop accross command ine options */
91 /*fprintf(stdout, "%s verbose: Parsing command line...\n", FuncName);*/
92 if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
93 usage_SurfExtrema(ps);
94 exit (0);
95 }
96
97 SUMA_SKIP_COMMON_OPTIONS(brk, kar);
98
99 if (!brk && (strcmp(argv[kar], "-debug") == 0))
100 {
101 if (kar+1 >= argc)
102 {
103 fprintf (SUMA_STDERR, "need a number after -debug \n");
104 exit (1);
105 }
106
107 Opt->debug = atoi(argv[++kar]);
108 brk = YUP;
109 }
110
111 if (!brk && (strcmp(argv[kar], "-node_debug") == 0))
112 {
113 if (kar+1 >= argc)
114 {
115 fprintf (SUMA_STDERR, "need a node index after -node_debug \n");
116 exit (1);
117 }
118
119 Opt->NodeDbg = atoi(argv[++kar]);
120 brk = YUP;
121 }
122
123 if (!brk && (strcmp(argv[kar], "-table") == 0))
124 {
125 if (kar+1 >= argc)
126 {
127 fprintf (SUMA_STDERR, "need a string after -table \n");
128 exit (1);
129 }
130
131 Opt->s = SUMA_copy_string(argv[++kar]);
132 brk = YUP;
133 }
134
135 if (!brk && (strcmp(argv[kar], "-prefix") == 0))
136 {
137 if (kar+1 >= argc)
138 {
139 fprintf (SUMA_STDERR, "need a dset prefix after -prefix \n");
140 exit (1);
141 }
142
143 Opt->out_prefix = SUMA_RemoveDsetExtension_s(argv[++kar],
144 SUMA_NO_DSET_FORMAT);
145 brk = YUP;
146 }
147
148 if (!brk && (strcmp(argv[kar], "-thresh") == 0))
149 {
150 if (kar+1 >= argc)
151 {
152 fprintf (SUMA_STDERR, "need a value after -thresh \n");
153 exit (1);
154 }
155
156 Opt->t = (float)strtod(argv[++kar],NULL);
157 brk = YUP;
158 }
159
160 if (!brk && (strcmp(argv[kar], "-gthresh") == 0))
161 {
162 if (kar+1 >= argc)
163 {
164 fprintf (SUMA_STDERR, "need a value after -gthresh \n");
165 exit (1);
166 }
167
168 Opt->t2 = (float)strtod(argv[++kar],NULL);
169 brk = YUP;
170 }
171
172 if (!brk && (strcmp(argv[kar], "-gscale") == 0))
173 {
174 if (kar+1 >= argc)
175 {
176 fprintf (SUMA_STDERR, "need a value after -gscale \n");
177 exit (1);
178 }
179 ++kar;
180 if (!strcmp(argv[kar],"NONE")) Opt->b1=SUMA_NO_GRAD_SCALE;
181 else if (!strcmp(argv[kar],"LMEAN")) Opt->b1=SUMA_MEAN_GRAD_SCALE;
182 else if (!strcmp(argv[kar],"GMEAN"))
183 Opt->b1=SUMA_GMEAN_GRAD_SCALE;
184 else {
185 SUMA_S_Errv("Bad value (%s) for option -gscale\n", argv[kar]);
186 exit (1);
187 }
188 brk = YUP;
189 }
190
191 if (!brk && (strcmp(argv[kar], "-extype") == 0))
192 {
193 if (kar+1 >= argc)
194 {
195 fprintf (SUMA_STDERR, "need a string after -extype \n");
196 exit (1);
197 }
198 ++kar;
199 if (!strcmp(argv[kar],"MIN")) Opt->b2=SUMA_MINIMUS;
200 else if (!strcmp(argv[kar],"MAX")) Opt->b2=SUMA_MAXIMUS;
201 else if (!strcmp(argv[kar],"ABS")) Opt->b2=SUMA_EXTREMUS;
202 else {
203 SUMA_S_Errv("Bad value (%s) for option -gscale\n", argv[kar]);
204 exit (1);
205 }
206 brk = YUP;
207 }
208
209
210
211 if (!brk && (strcmp(argv[kar], "-hood") == 0 ||
212 strcmp(argv[kar], "-nbhd_rad") == 0))
213 {
214 if (kar+1 >= argc)
215 {
216 fprintf (SUMA_STDERR, "need a value after -nbhd_rad \n");
217 exit (1);
218 }
219
220 Opt->r = atof(argv[++kar]);
221 if (Opt->r <= 0.0) {
222 SUMA_S_Errv("neighborhood radius is not valid (have %f from %s).\n",
223 Opt->r, argv[kar]);
224 exit (1);
225 }
226 brk = YUP;
227 }
228
229 if (!brk && !ps->arg_checked[kar]) {
230 SUMA_S_Errv("Option %s not understood. Try -help for usage\n",
231 argv[kar]);
232 suggest_best_prog_option(argv[0], argv[kar]);
233 exit(1);
234 } else {
235 brk = NOPE;
236 kar ++;
237 }
238 }
239
240 if (!Opt->out_prefix) {
241 Opt->out_prefix = SUMA_copy_string("SurfExtrema");
242 }
243 if (Opt->r <= 0.0) {
244 SUMA_S_Errv("neighborhood radius is not set (have %f).\n", Opt->r);
245 exit (1);
246 }
247
248 SUMA_RETURN(Opt);
249 }
250
251
main(int argc,char * argv[])252 int main (int argc,char *argv[])
253 {/* Main */
254 static char FuncName[]={"SurfExtrema"};
255 SUMA_GENERIC_PROG_OPTIONS_STRUCT *Opt;
256 SUMA_GENERIC_ARGV_PARSE *ps=NULL;
257 SUMA_DSET_FORMAT iform = SUMA_NO_DSET_FORMAT, oform = SUMA_NO_DSET_FORMAT;
258 SUMA_DSET *din=NULL, *dout=NULL, *doute=NULL;
259 SUMA_SurfSpecFile *Spec = NULL;
260 int i, N_Spec, N_inmask = -1, exists=0;
261 SUMA_SurfaceObject *SO=NULL, *SOf=NULL;
262 char *ooo=NULL, *s=NULL;
263 SUMA_Boolean LocalHead = NOPE;
264
265 SUMA_STANDALONE_INIT;
266 SUMA_mainENTRY;
267
268 /* Allocate space for DO structure */
269 SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct (SUMA_MAX_DISPLAYABLE_OBJECTS);
270 ps = SUMA_Parse_IO_Args(argc, argv, "-i;-t;-spec;-sv;-m;-dset;-talk;");
271
272 if (argc < 2) {
273 usage_SurfExtrema(ps);
274 exit (1);
275 }
276
277 Opt = SUMA_SurfExtrema_ParseInput (argv, argc, ps);
278
279 SUMA_CHECK_OUTPUT_SDSET_STATUS(Opt->out_prefix,
280 Opt->ps->dsetname[0], oform,
281 NULL, ".grd", exists);
282 SUMA_CHECK_OUTPUT_SDSET_STATUS(Opt->out_prefix,
283 Opt->ps->dsetname[0], oform,
284 NULL, ".ext", exists);
285
286 if (Opt->debug > 2) LocalHead = YUP;
287 Spec = SUMA_IO_args_2_spec(ps, &N_Spec);
288 if (N_Spec == 0) {
289 SUMA_S_Err("No surfaces found.");
290 exit(1);
291 }
292 if (N_Spec != 1) {
293 SUMA_S_Err("Multiple spec at input.");
294 exit(1);
295 }
296
297 SUMA_LH("Loading surface...");
298 SO = SUMA_Load_Spec_Surf(Spec, 0, ps->sv[0], Opt->debug);
299 if (!SO) {
300 fprintf (SUMA_STDERR,"Error %s:\n"
301 "Failed to find surface\n"
302 "in spec file. \n",
303 FuncName );
304 exit(1);
305
306 }
307 if (Spec->N_Surfs == 2) {
308 SOf = SUMA_Load_Spec_Surf(Spec, 1, ps->sv[0], Opt->debug);
309 if (!SOf) {
310 fprintf (SUMA_STDERR,"Error %s:\n"
311 "Failed to find surface\n"
312 "in spec file. \n",
313 FuncName );
314 exit(1);
315 }
316 } else { SOf = NULL; }
317
318 if (!(Opt->nmask = SUMA_load_all_command_masks(Opt->ps->bmaskname,
319 Opt->ps->nmaskname, Opt->ps->cmask, SO->N_Node,
320 &N_inmask)) && N_inmask < 0) {
321 SUMA_S_Err("Failed loading mask");
322 exit(1);
323 }
324
325 if (Opt->ps->N_dsetname != 1 && Opt->ps->N_dsetname != 0) {
326 SUMA_S_Errv("Need one and only one dset please."
327 "Have %d on command line.\n", Opt->ps->N_dsetname);
328 exit(1);
329 }
330
331 if (Opt->ps->N_dsetname) {
332 if (!(din = SUMA_LoadDset_s (Opt->ps->dsetname[0], &iform, 0))) {
333 SUMA_S_Errv("Failed to load dset named %s\n", Opt->ps->dsetname[0]);
334 exit(1);
335 }
336 } else {
337 SUMA_S_Note("No input datasets, will use surface convexity as input.");
338 /* Use convexity as input dataset*/
339 if (!SUMA_SurfaceMetrics(SO, "Convexity", NULL)) {
340 SUMA_S_Errv("Failed to compute convexity for %s\n",
341 SO->Label);
342 exit(1);
343 }
344 if (!(din = (SUMA_DSET *)SUMA_GetCx(SO->idcode_str,
345 SUMAg_CF->DsetList, 1))) {
346 SUMA_S_Errv("Could not get convexity dset for surface %s\n",
347 SO->Label);
348 exit(1);
349 }
350 }
351
352 if (!(dout = SUMA_DsetAvgGradient(SO, NULL, din, Opt->nmask,
353 0, Opt->b1))) {
354 SUMA_S_Err("Failed in SUMA_DsetAvgGradient");
355 exit(1);
356 }
357 /* write it out */
358 s = SUMA_OutputDsetFileStatus(Opt->out_prefix, NULL, &oform,
359 NULL, ".grd", &exists);
360 SUMA_AddNgrHist(dout->ngr, FuncName, argc, argv);
361 ooo = SUMA_WriteDset_s(s, dout, oform,
362 THD_ok_overwrite(), 0);
363 SUMA_free(ooo); ooo=NULL; SUMA_free(s); s = NULL;
364
365 if (!(doute = SUMA_DsetExtrema(SO, NULL, din, dout,
366 Opt->r, Opt->t, Opt->t2,
367 Opt->nmask, 0, Opt->b2, Opt->s))) {
368 SUMA_S_Err("Failed in SUMA_DsetAvgGradient");
369 exit(1);
370 }
371 /* write it out */
372 s = SUMA_OutputDsetFileStatus(Opt->out_prefix, NULL, &oform,
373 NULL, ".ext", &exists);
374 SUMA_AddNgrHist(doute->ngr, FuncName, argc, argv);
375 ooo = SUMA_WriteDset_s(s, doute, oform,
376 THD_ok_overwrite(), 0);
377 SUMA_free(ooo); ooo=NULL; SUMA_free(s); s = NULL;
378
379
380 if (ps) SUMA_FreeGenericArgParse(ps); ps = NULL;
381 if (Opt) Opt = SUMA_Free_Generic_Prog_Options_Struct(Opt);
382 if (!SUMA_Free_CommonFields(SUMAg_CF))
383 SUMA_error_message(FuncName,"SUMAg_CF Cleanup Failed!",1);
384
385 exit(0);
386
387 }
388