1 /*----------------------------------------------------------------------
2 * main functions for afni (and others):
3 *
4 * v2s_results * afni_vol2surf - create surface data
5 * int free_v2s_results - free surface data
6 *
7 * main interface routing:
8 *
9 * v2s_results * vol2surf - create surface data
10 * (assumes valid parameters)
11 *
12 * display functions:
13 *
14 * int disp_mri_imarr ( char * info, MRI_IMARR * dp );
15 * int disp_v2s_opts_t ( char * info, v2s_opts_t * sopt );
16 * int disp_v2s_param_t ( char * info, v2s_param_t * p );
17 * int disp_v2s_plugin_opts ( char * info, v2s_plugin_opts * d );
18 * int disp_v2s_results ( char * mesg, v2s_results * d );
19 *
20 * Author: R Reynolds
21 *----------------------------------------------------------------------
22 */
23
24 #define _VOL2SURF_C_ /* so the header files know */
25
26 char gv2s_history[] =
27 "----------------------------------------------------------------------\n"
28 "vol2surf library history:\n"
29 "\n"
30 "September 01, 2004 [rickr]\n"
31 " - initial install into afni\n"
32 "\n"
33 "September 02, 2004 [rickr]\n"
34 " - moved gv2s_map_names here (from SUMA_3dVol2Surf.c)\n"
35 " - moved v2s_map_type (and test) here (from SUMA_3dVol2Surf.c)\n"
36 " - define _VOL2SURF_C_, to allow extern defines in vol2surf.h\n"
37 "\n"
38 "September 09, 2004 [rickr]\n"
39 " - in afni_vol2surf(), print v2s options when debug > 1\n"
40 " - allow (first_node > last_node) if (last == 0), then change to n-1\n"
41 "\n"
42 "September 16, 2004 [rickr]\n"
43 " - added support for -gp_index, computing over a single sub-brick\n"
44 " - altered subs in dump_surf_3dt(), max_index in set_surf_results(),\n"
45 " set brick_index in v2s_apply_filter(), mem in alloc_output_mem(),\n"
46 " and added an index check in validate_v2s_inputs()\n"
47 " - add gp_index as a parameter of afni_vol2surf()\n"
48 " - changed keep_norm_dir to norm_dir, allowing default/keep/reverse\n"
49 "\n"
50 "September 29, 2004 [rickr]\n"
51 " - added thd_mask_from_brick() (to make byte mask from sub-brick)\n"
52 " - added compact_results(), in case nalloc > nused\n"
53 " - added realloc_ints() and realloc_vals_list() (for compact_results())\n"
54 " - in afni_vol2surf(), if 1 surf and no norms, set steps to 1\n"
55 " - in set_surf_results(), pass gp_index to v2s_apply_filter\n"
56 " - segment oob only if both nodes are\n"
57 " - move dset bounds to range_3dmm struct\n"
58 " - in segment_imarr()\n"
59 " - changed THD_3dmm_to_3dind() to new THD_3dmm_to_3dind_no_wod()\n"
60 " - verify success of THD_extract_series()\n"
61 " - keep track of repeated and oob nodes\n"
62 " - in init_seg_endpoints(), nuke p1, pn; save dicom_to_mm until end\n"
63 " - changed THD_3dmm_to_3dind() to new THD_3dmm_to_3dind_no_wod()\n"
64 " - added function v2s_is_good_map_index()\n"
65 "\n"
66 "October 08, 2004 [rickr]\n"
67 " - added disp_v2s_plugin_opts()\n"
68 " - dealt with default v2s mapping of surface pairs\n"
69 " - added v2s_fill_sopt_default()\n"
70 " - moved v2s_write_outfile_*() here, with print_header()\n"
71 " - in afni_vol2surf(), actually write output files\n"
72 "\n"
73 "October 25, 2004 [rickr]\n"
74 " - apply debug and dnode, even for defaults\n"
75 " - if the user sets dnode, then skip any (debug > 0) tests for it\n"
76 " - check for out of bounds, even if an endpoint is in (e.g. midpoint)\n"
77 " - in thd_mask_from_brick(), allow for absolute thresholding\n"
78 "\n"
79 "March 22, 2005 [rickr] - removed tabs\n"
80 "March 28, 2006 [rickr] - fixed mode computation\n"
81 "June 30, 2006 [rickr] - segc_file functionality (-save_seg_coords)\n"
82 "\n"
83 "August 9, 2006 [rickr]\n"
84 " - create argc, argv from options in v2s_make_command()\n"
85 " - added loc_add_2_list() and v2s_free_cmd() for v2s_make_command()\n"
86 " - added labels, thres index/value and surf vol dset to gv2s_plug_opts\n"
87 "\n"
88 "August 23, 2006 [rickr]\n"
89 " - in v2s_make_command(), change -skip_col_NSD to -outcols_afni_NSD\n"
90 " - in v2s_write_outfile_NSD(), only output node list if it exists\n"
91 " - do not let set_sparse_data_attribs() set nodes_from_dset attrib\n"
92 "September 6, 2006 [rickr]\n"
93 " - use NI_free() with NI_search_group_shallow()\n"
94 "November 10, 2006 [rickr]\n"
95 " - added thd_multi_mask_from_brick()\n"
96 "---------------------------------------------------------------------\n";
97
98 #include "mrilib.h"
99 #include "vol2surf.h"
100 #include "suma_objs.h" /* 21 Apr 2020 */
101
102 /*----------------------------------------------------------------------*/
103 /* local typedefs */
104 typedef struct
105 {
106 int nused;
107 int nalloc;
108 float * list;
109 } float_list_t;
110
111 typedef struct
112 {
113 THD_3dim_dataset * dset; /* for data and geometry */
114 THD_fvec3 p1; /* segment endpoints */
115 THD_fvec3 pn;
116 THD_fvec3 dset_min; /* bounds on the dataset */
117 THD_fvec3 dset_max;
118 int oob_check; /* should we check for oob? */
119 int debug; /* for local control */
120 } range_3dmm;
121
122 typedef struct
123 {
124 MRI_IMARR ims; /* the image array struct */
125 int repeats; /* number of repeated nodes */
126 int masked; /* number of masked points */
127 int oob; /* number of oob points */
128 int ifirst; /* 1D index of first point */
129 THD_ivec3 i3first; /* i3ind index of first point */
130 THD_ivec3 * i3arr; /* i3ind index array */
131 } range_3dmm_res;
132
133 /*----------------------------------------------------------------------*/
134 /* local prototypes */
135 static v2s_results * alloc_output_mem (v2s_opts_t *sopt, v2s_param_t *p);
136
137 static int alloc_ints(int ** ptr, int length, char * dstr, int debug);
138 static int alloc_vals_list(float *** ptr, int length, int width, int debug);
139 static int check_SUMA_surface( SUMA_surface * s );
140 static int compact_results(v2s_results * sd, int debug);
141 static float directed_dist(float * pnew, float * pold, float *dir, float dist);
142 static float dist_f3mm( THD_fvec3 * p1, THD_fvec3 * p2 );
143 static int disp_range_3dmm( char * info, range_3dmm * dp );
144 static int disp_range_3dmm_res( char * info, range_3dmm_res * dp );
145 static int disp_surf_vals( char * mesg, v2s_results * sd, int node );
146 static int dump_surf_3dt(v2s_opts_t *sopt, v2s_param_t *p, v2s_results *sd);
147 static int f3mm_out_of_bounds(THD_fvec3 *cp, THD_fvec3 *min, THD_fvec3 *max);
148 static int float_list_alloc(float_list_t *f,int **ilist,int size,int trunc);
149 static int float_list_comp_mode(float_list_t *f, float *mode, int *nvals,
150 int *index);
151 static int float_list_comp_nzmode(float_list_t *f, float *mode, int *nvals,
152 int *index);
153 static int float_list_slow_sort(float_list_t * f, int * ilist);
154 static int init_seg_endpoints(v2s_opts_t * sopt, v2s_param_t * p,
155 range_3dmm * R, int node );
156 static int init_range_structs( range_3dmm * r3, range_3dmm_res * res3 );
157 static int loc_add_2_list( char *** list, int * nall, int * len, char * str);
158 static double magnitude_f( float * p, int length );
159 static int print_header(FILE *outfp, char *surf, char *map, v2s_results *sd);
160 static int realloc_ints( int ** ptr, int length, char * dstr, int debug );
161 static int realloc_vals_list(float ** ptr, int length, int width, int debug);
162 static int set_3dmm_bounds(THD_3dim_dataset *dset, THD_fvec3 *min,
163 THD_fvec3 *max);
164 static int set_all_surf_vals (v2s_results * sd, int node_ind, int vind,
165 int i, int j, int k, float fval);
166 static int set_output_labels(v2s_results * sd, v2s_param_t * p,
167 v2s_opts_t * sopt);
168 static int set_surf_results(v2s_param_t *p, v2s_opts_t *sopt,v2s_results *sd,
169 range_3dmm_res * r3res, int node, int findex);
170 static int segment_imarr(range_3dmm_res *res,range_3dmm *R,v2s_opts_t *sopt,
171 byte * cmask, FILE * cfp, int nindex);
172 static int v2s_adjust_endpts(v2s_opts_t *sopt, THD_fvec3 *p1, THD_fvec3 *pn);
173 static float v2s_apply_filter(range_3dmm_res *rr, v2s_opts_t *sopt, int index,
174 int * findex);
175 static int v2s_free_cmd(v2s_opts_t * sopt);
176 static int v2s_map_needs_sort(int map);
177 static int validate_v2s_inputs(v2s_opts_t * sopt, v2s_param_t * p);
178
179 /*----------------------------------------------------------------------*/
180 /* globals to be accessed by plugin and in afni_suma.c */
181 v2s_plugin_opts gv2s_plug_opts = {
182 0,0,0,0, /* ready, map_all, use0, use1 */
183 -1,-1,-1,-1, /* s0A, s0B, s1A, s1B */
184 -1,0.0, /* threshold index/value */
185 {NULL, NULL, NULL, NULL}, /* surface labels [4] */
186 NULL /* surf vol dataset */
187 };
188 /* this must match v2s_map_nums enum */
189 char * gv2s_map_names[] = { "none", "mask", "midpoint", "mask2", "ave",
190 "count", "min", "max", "max_abs", "seg_vals",
191 "median", "mode", "nzmode", "nzave",
192 "nzmin", "nzmax" };
193 char gv2s_no_label[] = "undefined";
194
195 /*----------------------------------------------------------------------
196 * afni_vol2surf - create v2s_results from gv2s_* afni globals
197 *
198 * input: gpar : AFNI dataset to be used as the grid parent
199 * gp_index : sub-brick selector
200 * sA : surface A structure
201 * sB : surface B structure
202 * mask : thresholding mask
203 * use_defaults : use default sopt structure
204 *
205 * output: sd : allocated v2s_results struct, with requested data
206 *
207 * This function is used to map data from an AFNI volume to a surface.
208 * These structures are expected to be complete.
209 *
210 * The function now relies on opt_vol2surf. 29 Apr 2009 [rickr]
211 *----------------------------------------------------------------------
212 */
afni_vol2surf(THD_3dim_dataset * gpar,int gp_index,SUMA_surface * sA,SUMA_surface * sB,byte * mask,int use_defaults)213 v2s_results * afni_vol2surf ( THD_3dim_dataset * gpar, int gp_index,
214 SUMA_surface * sA, SUMA_surface * sB,
215 byte * mask, int use_defaults )
216 {
217 v2s_opts_t * sopt, sopt_def;
218
219 ENTRY("afni_vol2surf");
220
221 if ( use_defaults )
222 {
223 sopt = &sopt_def;
224 v2s_fill_sopt_default(sopt, sB ? 2 : 1); /* 1 or 2 surfaces */
225
226 /* but apply any debug options */
227 sopt->debug = gv2s_plug_opts.sopt.debug;
228 sopt->dnode = gv2s_plug_opts.sopt.dnode;
229 }
230 else
231 sopt = &gv2s_plug_opts.sopt;
232
233 sopt->gp_index = gp_index;
234
235 RETURN(opt_vol2surf(gpar, sopt, sA, sB, mask));
236 }
237
238
239 /*----------------------------------------------------------------------
240 * opt_vol2surf - create v2s_results from v2s_plugin_opts struct
241 *
242 * Fill in v2s_param_t struct, call vol2surf, write any output files.
243 *
244 * input: gpar : AFNI dataset to be used as the grid parent
245 * sopt : vol2surf options struct
246 * sA : surface A structure
247 * sB : surface B structure
248 * mask : volume mask
249 *
250 * output: sd : allocated v2s_results struct, with requested data
251 *
252 * This function is used to map data from an AFNI volume to a surface.
253 * These structures are expected to be complete.
254 *----------------------------------------------------------------------
255 */
opt_vol2surf(THD_3dim_dataset * gpar,v2s_opts_t * sopt,SUMA_surface * sA,SUMA_surface * sB,byte * mask)256 v2s_results * opt_vol2surf (THD_3dim_dataset * gpar, v2s_opts_t * sopt,
257 SUMA_surface * sA, SUMA_surface * sB, byte * mask)
258 {
259 v2s_param_t P;
260 v2s_results * res;
261
262 ENTRY("opt_vol2surf");
263
264 if ( !gpar ) RETURN(NULL);
265 if ( check_SUMA_surface(sA) ) RETURN(NULL);
266 if ( sB && check_SUMA_surface(sB) ) RETURN(NULL);
267
268 /* now fill the param struct based on the inputs */
269 memset(&P, 0, sizeof(P));
270 P.gpar = gpar;
271 P.cmask = mask;
272 P.nvox = DSET_NVOX(gpar);
273 P.over_steps = v2s_vals_over_steps(sopt->map);
274 P.nsurf = sB ? 2 : 1;
275 P.surf[0] = *sA;
276 if ( sB ) P.surf[1] = *sB;
277
278 /* verify steps, in case the user has not selected 2 surfaces */
279 if ( P.nsurf == 1 && ! sopt->use_norms )
280 sopt->f_steps = 1;
281
282 if ( sopt->debug > 2 ) disp_v2s_opts_t("-- v2s options: ", sopt);
283
284 /* fire it up */
285
286 res = vol2surf(sopt, &P);
287
288 v2s_make_command(sopt, &P);
289 if( gv2s_plug_opts.sopt.debug > 1 ) disp_v2s_command(sopt);
290
291 /* if the user wants output files, here they are (don't error check) */
292 if( res && sopt->outfile_1D ) {
293 if( THD_is_file(sopt->outfile_1D) )
294 fprintf(stderr,"** over-writing 1D output file '%s'\n",
295 sopt->outfile_1D);
296 v2s_write_outfile_1D(sopt, res, P.surf[0].label);
297 }
298
299 if( res && sopt->outfile_niml )
300 {
301 if ( THD_is_file(sopt->outfile_niml) )
302 fprintf(stderr,"** over-writing niml output file '%s'\n",
303 sopt->outfile_niml);
304 v2s_write_outfile_NSD(res, sopt, &P, 0);
305 }
306
307 if( sopt->cmd.fake ) v2s_free_cmd( sopt );
308
309 RETURN(res);
310 }
311
312
313 /*----------------------------------------------------------------------
314 * vol2surf - produce a v2s_results surface dataset
315 *
316 * input: sopt : volume to surface options struct
317 * p : volume to surface parameter struct
318 *
319 * output: sd : allocated v2s_results struct, with requested data
320 *
321 * This function is used to map data from an AFNI volume to a surface.
322 * These structures are expected to be complete.
323 *----------------------------------------------------------------------
324 */
vol2surf(v2s_opts_t * sopt,v2s_param_t * p)325 v2s_results * vol2surf ( v2s_opts_t * sopt, v2s_param_t * p )
326 {
327 v2s_results * sd;
328 int rv;
329 ENTRY("vol2surf");
330
331 if ( sopt == NULL || p == NULL )
332 {
333 fprintf( stderr, "** smd_wo - bad params (%p,%p)\n", sopt, p );
334 RETURN(NULL);
335 }
336
337 if ( validate_v2s_inputs(sopt, p) )
338 RETURN(NULL);
339
340 if ( sopt->map == E_SMAP_INVALID )
341 {
342 fprintf(stderr,"** v2s wo: invalid map %d\n", sopt->map);
343 RETURN(NULL);
344 }
345
346 sd = alloc_output_mem( sopt, p );
347 if ( !sd ) RETURN(NULL);
348
349 if ( sopt->debug > 1 ) disp_v2s_param_t( "-d post alloc_output_mem : ", p );
350
351 rv = dump_surf_3dt( sopt, p, sd );
352
353 if ( compact_results(sd, sopt->debug) )
354 {
355 free_v2s_results(sd); /* free whatever didn't get burned */
356 RETURN(NULL);
357 }
358
359 if ( sopt->debug > 1 ) disp_v2s_results( "-d post surf creation : ", sd);
360
361 RETURN(sd);
362 }
363
364
365 /* compact_results - if nused < nalloc, realloc all pointers */
compact_results(v2s_results * sd,int debug)366 static int compact_results(v2s_results * sd, int debug)
367 {
368 int rv = 0, mem = 0;
369 ENTRY("compact_results");
370
371 if ( sd->nused > sd->nalloc ) /* should not happen, of course */
372 {
373 fprintf(stderr,"** cr: nused (%d) > nalloc (%d) !!\n",
374 sd->nused, sd->nalloc);
375 RETURN(-1);
376 }
377
378 if ( sd->nused == sd->nalloc ) RETURN(0); /* we're good */
379
380 /* otherwise, realloc everything */
381
382 sd->nalloc = sd->nused;
383
384 if ( sd->nodes )
385 {
386 mem += sizeof(int);
387 rv = realloc_ints(&sd->nodes, sd->nalloc, "nodes", debug);
388 }
389
390 if ( ! rv && sd->volind )
391 {
392 mem += sizeof(int);
393 rv = realloc_ints(&sd->volind, sd->nalloc, "volind", debug);
394 }
395
396 if ( ! rv && sd->i )
397 {
398 mem += sizeof(int);
399 rv = realloc_ints(&sd->i, sd->nalloc, "i", debug);
400 }
401
402 if ( ! rv && sd->j )
403 {
404 mem += sizeof(int);
405 rv = realloc_ints(&sd->j, sd->nalloc, "j", debug);
406 }
407
408 if ( ! rv && sd->k )
409 {
410 mem += sizeof(int);
411 rv = realloc_ints(&sd->k, sd->nalloc, "k", debug);
412 }
413
414 if ( ! rv && sd->nvals )
415 {
416 mem += sizeof(int);
417 rv = realloc_ints(&sd->nvals, sd->nalloc, "nvals", debug);
418 }
419
420 if ( ! rv )
421 {
422 mem += (sizeof(float) * sd->max_vals);
423 rv = realloc_vals_list(sd->vals, sd->nalloc, sd->max_vals, debug);
424 }
425
426 if ( rv ) RETURN(rv); /* if there was a failure, just leave */
427
428 mem *= sd->nalloc; /* now multiply be the array length */
429
430 if ( debug > 1 )
431 fprintf(stderr,"+d compact results: reallocated %d bytes down to %d\n",
432 sd->memory, mem);
433
434 sd->memory = mem;
435
436 RETURN(rv);
437 }
438
439
440 /*----------------------------------------------------------------------
441 * dump_surf_3dt - for each node index, get an appropriate node sampling,
442 * and compute and output results across sub-bricks
443 *----------------------------------------------------------------------
444 */
dump_surf_3dt(v2s_opts_t * sopt,v2s_param_t * p,v2s_results * sd)445 static int dump_surf_3dt( v2s_opts_t * sopt, v2s_param_t * p, v2s_results * sd )
446 {
447 range_3dmm_res r3mm_res;
448 range_3dmm r3mm;
449 float dist, min_dist, max_dist;
450 FILE * coord_fp = NULL;
451 int sub, nindex, findex = 0;
452 int oobc, oomc;
453 int oob1, oob2;
454
455 ENTRY("dump_surf_3dt");
456
457 if ( ! sopt || ! p || ! sd )
458 {
459 fprintf(stderr, "** ds3 : bad params (%p,%p,%p)\n", sopt, p, sd );
460 RETURN(-1);
461 }
462
463 /* possibly write to a coordinate output file 30 Jun 2006 [rickr] */
464 if ( sopt->segc_file )
465 {
466 coord_fp = fopen(sopt->segc_file, "w");
467 if ( !coord_fp ) /* complain, but continue */
468 fprintf(stderr,"** failed to open coord file '%s', continuing...\n", sopt->segc_file);
469 }
470
471 /* note the number of sub-bricks, unless the user has given just one */
472 init_range_structs( &r3mm, &r3mm_res ); /* to empty */
473 r3mm.dset = p->gpar;
474 set_3dmm_bounds( p->gpar, &r3mm.dset_min, &r3mm.dset_max );
475
476 if ( sopt->debug > 1 )
477 fprintf(stderr, "-d dset bounding box: (%f, %f, %f)\n"
478 " (%f, %f, %f)\n",
479 r3mm.dset_min.xyz[0],r3mm.dset_min.xyz[1],r3mm.dset_min.xyz[2],
480 r3mm.dset_max.xyz[0],r3mm.dset_max.xyz[1],r3mm.dset_max.xyz[2]);
481
482 min_dist = 9999.9; /* v2.3 */
483 max_dist = -1.0;
484 oobc = 0; /* init out-of-bounds counter */
485 oomc = 0; /* init out-of-mask counter */
486
487 /* note, NodeList elements are in dicomm mm orientation */
488
489 for ( nindex = sopt->first_node; nindex <= sopt->last_node; nindex++ )
490 {
491 init_seg_endpoints(sopt, p, &r3mm, nindex); /* segment endpoints */
492 v2s_adjust_endpts( sopt, &r3mm.p1, &r3mm.pn );
493
494 if ( r3mm.debug ) {
495 fprintf(stderr,"===== finished with debug node =====\n");
496 r3mm.debug = 0;
497 }
498
499 if ( nindex == sopt->dnode ) /* if we have dnode, forget debug */
500 r3mm.debug = sopt->debug > 0 ? sopt->debug : 1;
501
502 if ( r3mm.debug )
503 fprintf(stderr,"===== processing debug node %d =====\n", nindex);
504
505 /* if both points are outside our dataset, skip the pair v2.3 */
506 oob1 = f3mm_out_of_bounds( &r3mm.p1, &r3mm.dset_min, &r3mm.dset_max );
507 oob2 = f3mm_out_of_bounds( &r3mm.pn, &r3mm.dset_min, &r3mm.dset_max );
508 if ( oob1 && oob2 )
509 {
510 oobc++;
511 if ( sopt->oob.show )
512 if ( set_all_surf_vals( sd, nindex, sopt->oob.index,
513 sopt->oob.index, sopt->oob.index,
514 sopt->oob.index, sopt->oob.value) )
515 RETURN(1);
516 if ( nindex == sopt->dnode )
517 {
518 disp_surf_vals("-d debug node, out-of-bounds : ", sd, -1);
519 fprintf(stderr,"-d dnode coords: (%f, %f, %f)\n",
520 r3mm.p1.xyz[0], r3mm.p1.xyz[1], r3mm.p1.xyz[2]);
521 }
522 continue;
523 }
524 else
525 r3mm.oob_check = oob1 || oob2;
526
527 dist = dist_f3mm( &r3mm.p1, &r3mm.pn );
528 if ( dist < min_dist ) min_dist = dist;
529 if ( dist > max_dist ) max_dist = dist;
530
531 if ( segment_imarr(&r3mm_res,&r3mm,sopt,p->cmask,coord_fp,nindex) != 0 )
532 continue;
533
534 if ( r3mm_res.ims.num == 0 ) /* any good voxels in the bunch? */
535 {
536 /* oob or oom? */
537 if ( r3mm_res.oob == sopt->f_steps ) /* out of bounds */
538 {
539 oobc++;
540 if ( sopt->oob.show )
541 if ( set_all_surf_vals( sd, nindex, sopt->oob.index,
542 sopt->oob.index, sopt->oob.index,
543 sopt->oob.index, sopt->oob.value) )
544 RETURN(1);
545 if ( nindex == sopt->dnode )
546 disp_surf_vals("-d debug node, all out-of-bounds : ",sd,-1);
547 }
548 else /* then we consider it out of mask */
549 {
550 oomc++;
551 if ( sopt->oom.show )
552 if ( set_all_surf_vals( sd, nindex, r3mm_res.ifirst,
553 r3mm_res.i3first.ijk[0], r3mm_res.i3first.ijk[1],
554 r3mm_res.i3first.ijk[2], sopt->oom.value ) )
555 RETURN(1);
556 if ( nindex == sopt->dnode )
557 disp_surf_vals("-d debug node, out-of-mask : ", sd, -1);
558 }
559
560 if ( nindex == sopt->dnode )
561 fprintf(stderr,"-d dnode coords: (%f, %f, %f)\n",
562 r3mm.p1.xyz[0], r3mm.p1.xyz[1], r3mm.p1.xyz[2]);
563
564 continue; /* in any case */
565 }
566
567 /* get element 0, just for the findex */
568 (void)v2s_apply_filter(&r3mm_res, sopt, 0, &findex);
569
570 if ( set_surf_results(p, sopt, sd, &r3mm_res, nindex, findex) )
571 RETURN(-1);
572
573 /* clean up the MRI_IMARR struct, but don't free imarr */
574 for ( sub = 0; sub < r3mm_res.ims.num; sub++ )
575 {
576 mri_free(r3mm_res.ims.imarr[sub]);
577 r3mm_res.ims.imarr[sub] = NULL;
578 }
579 r3mm_res.ims.num = 0;
580 }
581
582 if ( r3mm.debug ) {
583 fprintf(stderr,"===== finished with debug node =====\n");
584 r3mm.debug = 0;
585 }
586
587 if ( sopt->debug > 0 ) /* v2.3 */
588 {
589 fprintf( stderr, "-d node pair dist (min,max) = (%f,%f)\n",
590 min_dist, max_dist );
591 fprintf( stderr, "-d out-of-bounds, o-o-mask counts : %d, %d (of %d)\n",
592 oobc, oomc, sopt->last_node - sopt->first_node + 1);
593 }
594
595 /* set up the labels before leaving */
596 set_output_labels(sd, p, sopt);
597
598 /* close any coord file */
599 if ( coord_fp ) fclose( coord_fp );
600
601 /* now we can free the imarr and voxel lists */
602 if ( r3mm_res.ims.nall > 0 )
603 {
604 free(r3mm_res.ims.imarr);
605 free(r3mm_res.i3arr);
606 r3mm_res.ims.nall = 0;
607 }
608
609 RETURN(0);
610 }
611
612
613 /***********************************************************************
614 * copy labels to v2s_results struct
615 ***********************************************************************
616 */
set_output_labels(v2s_results * sd,v2s_param_t * p,v2s_opts_t * sopt)617 static int set_output_labels(v2s_results * sd, v2s_param_t * p,
618 v2s_opts_t * sopt)
619 {
620 int c, nl;
621 char str[32], label[32];
622
623 ENTRY("set_output_labels");
624
625 if(!sd->labels){ fprintf(stderr,"** SOL: NULL labels!\n"); RETURN(1); }
626
627 if( sopt->gp_index >= 0 || p->over_steps ) /* get one label */
628 {
629 if(sd->nlab != sd->max_vals)
630 fprintf(stderr,"** SOL: nlabel mis-match: %d vs %d\n",
631 sd->nlab, sd->max_vals);
632
633 /* set label (prefix) from sub-brick, if possible */
634 nl = sopt->gp_index >= 0 ? sopt->gp_index : 0;
635 if( p->gpar->dblk->brick_lab && p->gpar->dblk->brick_lab[nl] )
636 MCW_strncpy(str, p->gpar->dblk->brick_lab[nl], 24);
637 else
638 sprintf(str, "vol %d\n", nl);
639
640 /* if one label use str, else index upto nlab */
641 if( sd->nlab == 1 ) sd->labels[0] = strdup(str);
642 else {
643 for( c = 0; c < sd->nlab; c++ ) {
644 sprintf(label,"%s #%d", str, c);
645 sd->labels[c] = strdup(label);
646 }
647 }
648 }
649 else /* use all sub-brick labels */
650 {
651 nl = DSET_NVALS(p->gpar);
652 if(sd->nlab != nl)
653 fprintf(stderr,"** SOL: %d labels != %d\n", sd->nlab, nl);
654 if(nl > sd->nlab) nl = sd->nlab;
655 for( c = 0; c < nl; c++ )
656 {
657 if( !p->gpar->dblk->brick_lab || !p->gpar->dblk->brick_lab[c] )
658 {
659 sprintf(str,"volume #%d\n",c);
660 sd->labels[c] = strdup(str);
661 }
662 else
663 sd->labels[c] = strdup(p->gpar->dblk->brick_lab[c]);
664 }
665 }
666
667 RETURN(0);
668 }
669
670
671 /***********************************************************************
672 * set_surf_results
673 ***********************************************************************
674 */
set_surf_results(v2s_param_t * p,v2s_opts_t * sopt,v2s_results * sd,range_3dmm_res * r3res,int node,int findex)675 static int set_surf_results(v2s_param_t *p, v2s_opts_t * sopt, v2s_results * sd,
676 range_3dmm_res * r3res, int node, int findex)
677 {
678 int i, j, k, volind;
679 int c, max_index;
680 ENTRY("set_surf_results");
681
682 if ( sd->nused >= sd->nalloc )
683 {
684 fprintf(stderr,"** ssr: nused (%d) >= nalloc (%d)!\n",
685 sd->nused, sd->nalloc);
686 RETURN(1);
687 }
688
689 i = r3res->i3arr[findex].ijk[0];
690 j = r3res->i3arr[findex].ijk[1];
691 k = r3res->i3arr[findex].ijk[2];
692
693 /* now get 3D and 1D coordinates */
694 volind = i + DSET_NX(p->gpar) * (j + DSET_NY(p->gpar) * k );
695
696 /* set everything but the values */
697 if (sd->nodes ) sd->nodes[sd->nused] = node;
698 if (sd->volind) sd->volind[sd->nused] = volind;
699 if (sd->i ) sd->i[sd->nused] = i;
700 if (sd->j ) sd->j[sd->nused] = j;
701 if (sd->k ) sd->k[sd->nused] = k;
702 if (sd->nvals ) sd->nvals[sd->nused] = r3res->ims.num;
703
704 /* set max_index, and adjust in case max_vals has been restricted */
705 max_index = p->over_steps ? r3res->ims.num : DSET_NVALS(p->gpar);
706 if ( max_index > sd->max_vals ) max_index = sd->max_vals;
707
708 if ( sopt->gp_index >= 0 )
709 sd->vals[0][sd->nused] =
710 v2s_apply_filter(r3res, sopt, sopt->gp_index, NULL);
711 else
712 for ( c = 0; c < max_index; c++ )
713 sd->vals[c][sd->nused] = v2s_apply_filter(r3res, sopt, c, NULL);
714
715 /* possibly fill line with default if by steps and short */
716 if ( max_index < sd->max_vals )
717 for ( c = max_index; c < sd->max_vals; c++ )
718 sd->vals[c][sd->nused] = 0.0;
719
720 /* maybe the user wants us to be verbose about this node */
721 if ( node == sopt->dnode )
722 {
723 fprintf(stderr,
724 "--------------------------------------------------\n");
725 if ( ! p->over_steps && sopt->gp_index >= 0 )
726 fprintf(stderr,"+d dnode %d gets %f from gp_index %d\n",
727 node, sd->vals[0][sd->nused], sopt->gp_index);
728 if ( sopt->debug > 1 )
729 fprintf(stderr, "-d debug: node %d, findex %d, vol_index %d\n",
730 node, findex, volind );
731 if ( sopt->use_norms )
732 {
733 float * fp = p->surf[0].norm[node].xyz;
734 fprintf(stderr,"-d normal %f, %f, %f\n", fp[0], fp[1], fp[2]);
735 }
736 disp_mri_imarr( "+d raw data: ", &r3res->ims );
737 }
738
739 sd->nused++;
740
741 RETURN(0);
742 }
743
744
745 /***********************************************************************
746 * set_all_surf_vals
747 * return 0 on success
748 ***********************************************************************
749 */
set_all_surf_vals(v2s_results * sd,int node_ind,int vind,int i,int j,int k,float fval)750 static int set_all_surf_vals( v2s_results * sd, int node_ind, int vind,
751 int i, int j, int k, float fval )
752 {
753 int c, nused;
754
755 ENTRY("set_all_surf_vals");
756
757 nused = sd->nused;
758
759 if ( nused >= sd->nalloc )
760 {
761 fprintf(stderr,"** sasv: nused=%d >= nalloc=%d!\n", nused, sd->nalloc);
762 RETURN(1);
763 }
764
765 if ( sd->nodes ) sd->nodes[nused] = node_ind;
766 if ( sd->volind ) sd->volind[nused] = vind;
767 if ( sd->i ) sd->i[nused] = i;
768 if ( sd->j ) sd->j[nused] = j;
769 if ( sd->k ) sd->k[nused] = k;
770 if ( sd->nvals ) sd->nvals[nused] = sd->max_vals;
771
772 for ( c = 0; c < sd->max_vals; c++ )
773 sd->vals[c][nused] = fval;
774
775 sd->nused++;
776
777 RETURN(0);
778 }
779
780
781 /***********************************************************************
782 * segment_imarr - get MRI_IMARR for steps over a segment
783 *
784 * The res->ims structure should be empty, except that it may
785 * optionally contain space for pointers in imarr. Note that nall
786 * should be accurate.
787 *
788 * return 0 on success
789 ***********************************************************************
790 */
segment_imarr(range_3dmm_res * res,range_3dmm * R,v2s_opts_t * sopt,byte * cmask,FILE * cfp,int nindex)791 static int segment_imarr( range_3dmm_res *res, range_3dmm *R, v2s_opts_t *sopt,
792 byte * cmask, FILE * cfp, int nindex )
793 {
794 THD_fvec3 f3mm;
795 THD_ivec3 i3ind;
796 float rat1, ratn;
797 int nx, ny, nvox;
798 int step, vindex, prev_ind;
799
800 ENTRY("segment_imarr");
801
802 /* check params for validity */
803 if ( !R || !sopt || !res || !R->dset )
804 {
805 fprintf(stderr, "** seg_imarr: invalid params (%p,%p,%p)\n",R,sopt,res);
806 if ( R ) disp_range_3dmm("segment_imarr: bad inputs:", R );
807 RETURN(-1);
808 }
809
810 if ( R->debug > 1 )
811 disp_range_3dmm("-d segment_imarr: ", R );
812
813 /* handle this as an acceptable, trivial case */
814 if ( sopt->f_steps < 1 )
815 {
816 res->ims.num = 0;
817 RETURN(0);
818 }
819
820 nx = DSET_NX(R->dset);
821 ny = DSET_NY(R->dset);
822 nvox = DSET_NVOX(R->dset);
823
824 /* if we don't have enough memory for results, (re)allocate some */
825 if ( res->ims.nall < sopt->f_steps )
826 {
827 if ( R->debug > 1 )
828 fprintf(stderr,"+d realloc of imarr (from %d to %d pointers)\n",
829 res->ims.nall,sopt->f_steps);
830
831 res->ims.nall = sopt->f_steps;
832 res->ims.imarr = realloc(res->ims.imarr,
833 sopt->f_steps*sizeof(MRI_IMAGE *));
834 res->i3arr = realloc(res->i3arr, sopt->f_steps*sizeof(THD_ivec3));
835 if ( !res->ims.imarr || !res->i3arr )
836 {
837 fprintf(stderr,"** seg_imarr: no memory for %d MRI_IMAGE ptrs\n",
838 sopt->f_steps);
839 res->ims.nall = res->ims.num = 0;
840 /* one might be good */
841 if ( res->ims.imarr ) free(res->ims.imarr);
842 if ( res->i3arr ) free(res->i3arr);
843 RETURN(-1);
844 }
845 }
846
847 /* init return structure */
848 res->ims.num = 0;
849 res->repeats = 0;
850 res->masked = 0;
851 res->oob = 0;
852 res->i3first = THD_3dmm_to_3dind_no_wod( R->dset, R->p1 );
853 res->ifirst = res->i3first.ijk[0] +
854 nx * (res->i3first.ijk[1] + ny * res->i3first.ijk[2] );
855
856 prev_ind = -1; /* in case we want unique voxels */
857
858 for ( step = 0; step < sopt->f_steps; step++ )
859 {
860 /* set our endpoint ratios */
861 if ( sopt->f_steps < 2 ) /* if this is based on a single point */
862 ratn = 0.0;
863 else
864 ratn = (float)step / (sopt->f_steps - 1);
865 rat1 = 1.0 - ratn;
866
867 f3mm.xyz[0] = rat1 * R->p1.xyz[0] + ratn * R->pn.xyz[0];
868 f3mm.xyz[1] = rat1 * R->p1.xyz[1] + ratn * R->pn.xyz[1];
869 f3mm.xyz[2] = rat1 * R->p1.xyz[2] + ratn * R->pn.xyz[2];
870
871 /* accept part being oob 30 Sep 2004 [rickr] */
872 if ( R->oob_check &&
873 f3mm_out_of_bounds( &f3mm, &R->dset_min, &R->dset_max ) )
874 {
875 res->oob++;
876 continue;
877 }
878
879 i3ind = THD_3dmm_to_3dind_no_wod( R->dset, f3mm );
880 vindex = i3ind.ijk[0] + nx * (i3ind.ijk[1] + ny * i3ind.ijk[2] );
881
882 /* is this voxel masked out? */
883 if ( cmask && !cmask[vindex] )
884 {
885 res->masked++;
886 continue;
887 }
888
889 /* is this voxel repeated, and if so, do we skip it? */
890 if ( sopt->f_index == V2S_INDEX_VOXEL && vindex == prev_ind )
891 {
892 res->repeats++;
893 continue;
894 }
895
896 /* Huston, we have a good voxel... */
897
898 /* now for the big finish, get and insert the actual data */
899
900 res->i3arr [res->ims.num] = i3ind; /* store the 3-D indices */
901 res->ims.imarr[res->ims.num] = THD_extract_series( vindex, R->dset, 0 );
902 res->ims.num++;
903
904 /* print 10 warnings, then wait for a long time
905 (done for PT - failing to process RGB dsets) 21 Sep 2015 [rickr] */
906 if ( ! res->ims.imarr[res->ims.num-1] )
907 {
908 static int nwarn = 0;
909 if( nwarn < 10 )
910 fprintf(stderr,"** failed THD_extract_series, vox %d\n", vindex);
911 if( nwarn == 9 )
912 fprintf(stderr," (no more 'failed' messages for a while)\n");
913
914 nwarn ++;
915 if( nwarn >= nvox ) nwarn = 0;
916
917 res->ims.num--; /* maybe we do not want to return a NULL? */
918 RETURN(-1);
919 }
920
921 if ( R->debug > 2 )
922 fprintf(stderr, "-d seg step %2d, vindex %d, coords %f %f %f\n",
923 step,vindex,f3mm.xyz[0],f3mm.xyz[1],f3mm.xyz[2]);
924 if ( cfp ) /* then write voxel coordinates to output file */
925 {
926 f3mm = THD_3dmm_to_dicomm(R->dset, f3mm); /* output in DICOM */
927 if ( prev_ind == -1 ) fprintf(cfp, "%8d", nindex);
928 fprintf(cfp, " %9.3f %9.3f %9.3f",
929 f3mm.xyz[0],f3mm.xyz[1],f3mm.xyz[2]);
930 }
931
932 prev_ind = vindex; /* note this for next time */
933
934 }
935
936 /* maybe write eol to coord file */
937 if ( cfp && prev_ind != -1 ) fputc('\n', cfp);
938
939 if ( R->debug > 1 )
940 disp_range_3dmm_res( "+d i3mm_seg_imarr results: ", res );
941
942 RETURN(0);
943 }
944
945
946 /*----------------------------------------------------------------------
947 * f3mm_out_of_bounds - check wether cp is between min and max
948 *
949 * - v2.3 [rickr]
950 *----------------------------------------------------------------------
951 */
f3mm_out_of_bounds(THD_fvec3 * cp,THD_fvec3 * min,THD_fvec3 * max)952 static int f3mm_out_of_bounds( THD_fvec3 * cp, THD_fvec3 * min, THD_fvec3 * max)
953 {
954 int c;
955
956 if ( !cp || !min || !max )
957 return(-1);
958
959 for ( c = 0; c < 3; c++ )
960 {
961 if ( ( cp->xyz[c] < min->xyz[c] ) ||
962 ( cp->xyz[c] > max->xyz[c] ) )
963 return(-1);
964 }
965
966 return(0);
967 }
968
969
970 /*----------------------------------------------------------------------
971 * v2s_adjust_endpoints - adjust endpoints for map and options
972 *
973 * return 0 on success
974 * < 0 on error
975 *----------------------------------------------------------------------
976 */
v2s_adjust_endpts(v2s_opts_t * sopt,THD_fvec3 * p1,THD_fvec3 * pn)977 static int v2s_adjust_endpts( v2s_opts_t * sopt, THD_fvec3 * p1, THD_fvec3 * pn)
978 {
979 THD_fvec3 f3_diff;
980 float dist, factor;
981
982 ENTRY("v2s_adjust_endpts");
983
984 if ( !sopt || !p1 || !pn )
985 {
986 fprintf(stderr,"** v2s_ae: invalid params (%p,%p,%p)\n", sopt, p1, pn);
987 RETURN(-1);
988 }
989
990 /* first, get the difference, and distance */
991 f3_diff.xyz[0] = pn->xyz[0] - p1->xyz[0];
992 f3_diff.xyz[1] = pn->xyz[1] - p1->xyz[1];
993 f3_diff.xyz[2] = pn->xyz[2] - p1->xyz[2];
994
995 dist = dist_f3mm( p1, pn );
996
997 if ( (sopt->f_p1_fr != 0.0) || (sopt->f_p1_mm != 0.0) )
998 {
999 if ( sopt->f_p1_fr != 0.0 ) /* what the heck, choose fr if both */
1000 factor = sopt->f_p1_fr;
1001 else
1002 factor = (dist == 0.0) ? 0.0 : sopt->f_p1_mm / dist;
1003
1004 p1->xyz[0] += factor * f3_diff.xyz[0];
1005 p1->xyz[1] += factor * f3_diff.xyz[1];
1006 p1->xyz[2] += factor * f3_diff.xyz[2];
1007 }
1008
1009 if ( (sopt->f_pn_fr != 0.0) || (sopt->f_pn_mm != 0.0) )
1010 {
1011 if ( sopt->f_pn_fr != 0.0 )
1012 factor = sopt->f_pn_fr;
1013 else
1014 factor = (dist == 0.0) ? 0.0 : sopt->f_pn_mm / dist;
1015
1016 pn->xyz[0] += factor * f3_diff.xyz[0];
1017 pn->xyz[1] += factor * f3_diff.xyz[1];
1018 pn->xyz[2] += factor * f3_diff.xyz[2];
1019 }
1020
1021 switch ( sopt->map )
1022 {
1023 default:
1024 fprintf(stderr,"** v2s_ae: mapping %d not ready\n", sopt->map );
1025 RETURN(-1);
1026
1027 case E_SMAP_AVE:
1028 case E_SMAP_MAX:
1029 case E_SMAP_MAX_ABS:
1030 case E_SMAP_MIN:
1031 case E_SMAP_MASK:
1032 case E_SMAP_SEG_VALS:
1033 case E_SMAP_MEDIAN:
1034 case E_SMAP_MODE:
1035 case E_SMAP_NZMODE:
1036 case E_SMAP_NZAVE:
1037 case E_SMAP_NZMIN:
1038 case E_SMAP_NZMAX:
1039 break;
1040
1041 case E_SMAP_MIDPT:
1042
1043 /* set the first point to be the average of the two */
1044 p1->xyz[0] = (p1->xyz[0] + pn->xyz[0]) / 2.0;
1045 p1->xyz[1] = (p1->xyz[1] + pn->xyz[1]) / 2.0;
1046 p1->xyz[2] = (p1->xyz[2] + pn->xyz[2]) / 2.0;
1047
1048 break;
1049 }
1050
1051 RETURN(0);
1052 }
1053
1054
1055 /*---------------------------------------------------------------------------
1056 * v2s_apply_filter - compute results for the given function and index
1057 *
1058 * As a side step, return any filter result index.
1059 *---------------------------------------------------------------------------
1060 */
v2s_apply_filter(range_3dmm_res * rr,v2s_opts_t * sopt,int index,int * findex)1061 static float v2s_apply_filter( range_3dmm_res * rr, v2s_opts_t * sopt,
1062 int index, int * findex )
1063 {
1064 static float_list_t flist = { 0, 0, NULL }; /* for sorting results */
1065 static int * ind_list = NULL; /* track index sources */
1066 double tmp, comp = 0.0;
1067 float fval;
1068 int count, source;
1069 int brick_index = 0;
1070
1071 ENTRY("v2s_apply_filter");
1072
1073 if ( !rr || !sopt || index < 0 )
1074 {
1075 fprintf(stderr,"** v2s_cm2: invalid params (%p,%p,%d)\n",
1076 rr, sopt, index);
1077 RETURN(0.0);
1078 }
1079
1080 if ( rr->ims.num <= 0 )
1081 RETURN(0.0);
1082
1083 /* if sorting is required for results, do it now */
1084 if ( v2s_map_needs_sort( sopt->map ) )
1085 {
1086 if ( float_list_alloc( &flist, &ind_list, rr->ims.num, 0 ) != 0 )
1087 RETURN(0.0);
1088 for ( count = 0; count < rr->ims.num; count++ )
1089 {
1090 flist.list[count] = MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
1091 ind_list [count] = count; /* init index sources */
1092 }
1093 flist.nused = rr->ims.num;
1094 float_list_slow_sort( &flist, ind_list );
1095 }
1096
1097 switch ( sopt->map )
1098 {
1099 default:
1100 if ( findex ) *findex = 0;
1101 RETURN(0.0);
1102
1103 case E_SMAP_AVE:
1104 if ( findex ) *findex = 0;
1105 for ( count = 0; count < rr->ims.num; count++ )
1106 comp += MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
1107
1108 comp = comp / rr->ims.num;
1109 break;
1110
1111 case E_SMAP_NZAVE:
1112 {
1113 int nnum;
1114
1115 if ( findex ) *findex = 0;
1116 nnum = 0;
1117 comp = 0.0;
1118 for ( count = 0; count < rr->ims.num; count++ ){
1119 tmp = MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
1120 if(tmp!=0){
1121 comp += tmp;
1122 nnum++;
1123 }
1124 }
1125 if (nnum != 0)
1126 comp = comp / nnum;
1127 break;
1128 }
1129
1130 case E_SMAP_MASK:
1131 case E_SMAP_MIDPT:
1132 if ( findex ) *findex = 0;
1133 /* we have only the one point */
1134 comp = MRI_FLOAT_PTR(rr->ims.imarr[0])[index];
1135 break;
1136
1137 case E_SMAP_MAX:
1138 comp = MRI_FLOAT_PTR(rr->ims.imarr[0])[index];
1139 if ( findex ) *findex = 0;
1140
1141 for ( count = 1; count < rr->ims.num; count++ )
1142 {
1143 tmp = MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
1144 if ( tmp > comp )
1145 {
1146 if ( findex ) *findex = count;
1147 comp = tmp;
1148 }
1149 }
1150 break;
1151
1152 case E_SMAP_NZMAX:
1153 comp = MRI_FLOAT_PTR(rr->ims.imarr[0])[index];
1154 if ( findex ) *findex = 0;
1155
1156 for ( count = 1; count < rr->ims.num; count++ )
1157 {
1158 tmp = MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
1159 if (( tmp > comp ) && (tmp != 0.0))
1160 {
1161 if ( findex ) *findex = count;
1162 comp = tmp;
1163 }
1164 }
1165 break;
1166
1167 case E_SMAP_MAX_ABS:
1168 comp = MRI_FLOAT_PTR(rr->ims.imarr[0])[index];
1169 if ( findex ) *findex = 0;
1170
1171 for ( count = 1; count < rr->ims.num; count++ )
1172 {
1173 tmp = MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
1174 if ( fabs(tmp) > fabs(comp) )
1175 {
1176 if ( findex ) *findex = count;
1177 comp = tmp;
1178 }
1179 }
1180 break;
1181
1182 case E_SMAP_MIN:
1183 comp = MRI_FLOAT_PTR(rr->ims.imarr[0])[index];
1184 if ( findex ) *findex = 0;
1185
1186 for ( count = 1; count < rr->ims.num; count++ )
1187 {
1188 tmp = MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
1189 if ( tmp < comp )
1190 {
1191 if ( findex ) *findex = count;
1192 comp = tmp;
1193 }
1194 }
1195 break;
1196
1197 case E_SMAP_NZMIN:
1198 comp = MRI_FLOAT_PTR(rr->ims.imarr[0])[index];
1199 if ( findex ) *findex = 0;
1200
1201 for ( count = 1; count < rr->ims.num; count++ )
1202 {
1203 tmp = MRI_FLOAT_PTR(rr->ims.imarr[count])[index];
1204 if (( tmp < comp ) && (tmp != 0.0))
1205 {
1206 if ( findex ) *findex = count;
1207 comp = tmp;
1208 }
1209 }
1210 break;
1211
1212 case E_SMAP_SEG_VALS:
1213 /* if the user has specified a brick index, use it */
1214 if ( sopt->gp_index >= 0 ) brick_index = sopt->gp_index;
1215 if ( findex ) *findex = 0;
1216 comp = MRI_FLOAT_PTR(rr->ims.imarr[index])[brick_index];
1217 break;
1218
1219 case E_SMAP_MEDIAN:
1220 count = flist.nused >> 1;
1221 if ( (flist.nused & 1) || (count == 0) )
1222 {
1223 comp = flist.list[count];
1224 if ( findex ) *findex = ind_list[count];
1225 }
1226 else
1227 {
1228 comp = (flist.list[count-1] + flist.list[count]) / 2;
1229 if ( findex ) *findex = ind_list[count-1]; /* take first */
1230 }
1231 break;
1232
1233 case E_SMAP_MODE:
1234 float_list_comp_mode(&flist, &fval, &count, &source);
1235 comp = fval;
1236 if ( findex ) *findex = ind_list[source];
1237 break;
1238
1239 case E_SMAP_NZMODE:
1240 float_list_comp_nzmode(&flist, &fval, &count, &source);
1241 comp = fval;
1242 if ( findex ) *findex = ind_list[source];
1243 break;
1244 }
1245
1246 RETURN((float)comp);
1247 }
1248
1249
1250 /*----------------------------------------------------------------------
1251 * v2s_map_needs_sort - does this map function require sorting?
1252 *
1253 * return 1 on true
1254 * 0 on false
1255 *----------------------------------------------------------------------
1256 */
v2s_map_needs_sort(int map)1257 static int v2s_map_needs_sort( int map )
1258 {
1259 if ( (map == E_SMAP_MEDIAN) || (map == E_SMAP_MODE) || (map == E_SMAP_NZMODE))
1260 return 1;
1261
1262 return 0;
1263 }
1264
1265
1266 /*----------------------------------------------------------------------
1267 * float_list_comp_mode - compute the mode of the sorted list
1268 *
1269 * return 0 : on success
1270 * -1 : on error
1271 *----------------------------------------------------------------------
1272 */
float_list_comp_mode(float_list_t * f,float * mode,int * nvals,int * index)1273 static int float_list_comp_mode( float_list_t *f, float *mode, int *nvals,
1274 int *index )
1275 {
1276 float fcur;
1277 int ncur, c;
1278
1279 ENTRY("float_list_comp_mode");
1280
1281 /* init default results */
1282 *nvals = ncur = 1;
1283 *mode = fcur = f->list[0];
1284 *index = 0;
1285
1286 for ( c = 1; c < f->nused; c++ )
1287 {
1288 if ( f->list[c] == fcur )
1289 ncur ++;
1290 else /* found a new entry to count */
1291 {
1292 if ( ncur > *nvals ) /* keep track of any new winner */
1293 {
1294 *mode = fcur;
1295 *nvals = ncur;
1296 *index = c - ncur; /* note first occurrence */
1297 }
1298
1299 fcur = f->list[c];
1300 ncur = 1;
1301 }
1302 }
1303
1304 if ( ncur > *nvals ) /* keep track of any new winner */
1305 {
1306 *mode = fcur;
1307 *nvals = ncur;
1308 *index = c - ncur; /* note first occurrence */
1309 }
1310
1311 RETURN(0);
1312 }
1313
1314
1315
1316 /*----------------------------------------------------------------------
1317 * float_list_comp_nzmode - compute the non-zero mode of the sorted list
1318 *
1319 * return 0 : on success
1320 * -1 : on error
1321 *----------------------------------------------------------------------
1322 */
float_list_comp_nzmode(float_list_t * f,float * mode,int * nvals,int * index)1323 static int float_list_comp_nzmode( float_list_t *f, float *mode, int *nvals,
1324 int *index )
1325 {
1326 float fcur;
1327 int ncur, c;
1328
1329 ENTRY("float_list_comp_nzmode");
1330
1331 /* init default results */
1332 if(f->list[0] != 0.0) {
1333 *nvals = ncur = 1;
1334 }
1335 else {
1336 *nvals = ncur = 0;
1337 }
1338
1339 *mode = fcur = f->list[0];
1340 *index = 0;
1341
1342 for ( c = 1; c < f->nused; c++ )
1343 {
1344 if ( f->list[c] == fcur )
1345 ncur ++;
1346 else /* found a new entry to count */
1347 {
1348 if (( ncur > *nvals ) && (fcur != 0.0)) /* keep track of any new winner */
1349 {
1350 *mode = fcur;
1351 *nvals = ncur;
1352 *index = c - ncur; /* note first occurrence */
1353 }
1354
1355 fcur = f->list[c];
1356 ncur = 1;
1357 }
1358 }
1359
1360
1361 if (( ncur > *nvals ) && (fcur != 0.0)) /* keep track of any new winner */
1362 {
1363 *mode = fcur;
1364 *nvals = ncur;
1365 *index = c - ncur; /* note first occurrence */
1366 }
1367 RETURN(0);
1368 }
1369
1370
1371 /*----------------------------------------------------------------------
1372 * float_list_slow_sort - sort (small) float list
1373 *
1374 * If ilist, track index sources.
1375 *
1376 * return 0 on success
1377 * < 0 on error
1378 *----------------------------------------------------------------------
1379 */
float_list_slow_sort(float_list_t * f,int * ilist)1380 static int float_list_slow_sort( float_list_t * f, int * ilist )
1381 {
1382 float * list, save;
1383 int c0, c1, sindex;
1384
1385 ENTRY("float_list_slow_sort");
1386
1387 list = f->list; /* for any little speed gain */
1388
1389 for ( c0 = 0; c0 < f->nused-1; c0++ )
1390 {
1391 sindex = c0;
1392 save = list[c0];
1393
1394 /* find smallest remaining */
1395 for ( c1 = c0+1; c1 < f->nused; c1++ )
1396 if ( list[c1] < save )
1397 {
1398 sindex = c1;
1399 save = list[sindex];
1400 }
1401
1402 /* swap if smaller was found */
1403 if ( sindex > c0 )
1404 {
1405 list[sindex] = list[c0];
1406 list[c0] = save;
1407
1408 if ( ilist ) /* make same swap of indices */
1409 {
1410 c1 = ilist[sindex];
1411 ilist[sindex] = ilist[c0];
1412 ilist[c0] = c1;
1413 }
1414 }
1415 }
1416
1417 RETURN(0);
1418 }
1419
1420
1421 /*----------------------------------------------------------------------
1422 * float_list_alloc - verify float list memory
1423 *
1424 * If trunc(ate), realloc down to necessary size.
1425 * If ilist, make space for accompanying int list.
1426 *
1427 * return 0 on success
1428 * < 0 on error
1429 *----------------------------------------------------------------------
1430 */
float_list_alloc(float_list_t * f,int ** ilist,int size,int trunc)1431 static int float_list_alloc(float_list_t *f, int **ilist, int size, int trunc)
1432 {
1433 ENTRY("float_list_alloc");
1434
1435 if ( (f->nalloc < size) ||
1436 (trunc && (f->nalloc > size)) )
1437 {
1438 f->list = (float *)realloc(f->list, size * sizeof(float));
1439 if ( ! f->list )
1440 {
1441 fprintf(stderr,"** float_list_alloc: failed for %d floats\n", size);
1442 RETURN(-2);
1443 }
1444 f->nalloc = size;
1445
1446 if ( ilist ) /* then allocate accompanying ilist */
1447 {
1448 *ilist = (int *)realloc(*ilist, size * sizeof(int));
1449 if ( ! *ilist )
1450 {
1451 fprintf(stderr,"** float_list_alloc: failed for %d ints\n",
1452 size);
1453 RETURN(-2);
1454 }
1455 }
1456
1457 if ( trunc && (f->nused > f->nalloc) )
1458 f->nused = f->nalloc;
1459 }
1460
1461 RETURN(0);
1462 }
1463
1464
1465 /*---------------------------------------------------------------------------
1466 * directed_dist - travel from pold, along dir, a distance of dist
1467 *---------------------------------------------------------------------------
1468 */
directed_dist(float * pnew,float * pold,float * dir,float dist)1469 static float directed_dist(float * pnew, float * pold, float * dir, float dist)
1470 {
1471 double mag, ratio;
1472 int c;
1473
1474 ENTRY("directed_dist");
1475
1476 mag = magnitude_f(dir, 3);
1477
1478 if ( mag < V2S_EPSILON ) /* can't be negative */
1479 ratio = 0.0;
1480 else
1481 ratio = dist/mag;
1482
1483 for ( c = 0; c < 3; c++ )
1484 pnew[c] = pold[c] + dir[c]*ratio;
1485
1486 RETURN(dist);
1487 }
1488
1489
1490 /*----------------------------------------------------------------------
1491 * init_seg_endpoints - initialize segment endpoints
1492 *----------------------------------------------------------------------
1493 */
init_seg_endpoints(v2s_opts_t * sopt,v2s_param_t * p,range_3dmm * R,int node)1494 static int init_seg_endpoints ( v2s_opts_t * sopt, v2s_param_t * p,
1495 range_3dmm * R, int node )
1496 {
1497 SUMA_surface * sp;
1498 float * norm;
1499
1500 ENTRY("init_seg_endpoints");
1501
1502 /* get node from first surface */
1503 sp = p->surf;
1504 R->p1.xyz[0] = sp->ixyz[node].x;
1505 R->p1.xyz[1] = sp->ixyz[node].y;
1506 R->p1.xyz[2] = sp->ixyz[node].z;
1507
1508 /* set pn based on other parameters */
1509 if ( sopt->use_norms ) {
1510 /* actually reverse norm direction 4 Feb 2009 [rickr] */
1511 norm = sp->norm[node].xyz;
1512 if ( sopt->norm_dir == V2S_NORM_REVERSE )
1513 {
1514 if( node == sopt->dnode )
1515 fprintf(stderr,"-d reversing norm dir at dnode %d\n", node);
1516 directed_dist(R->pn.xyz, R->p1.xyz, norm, -sopt->norm_len);
1517 } else
1518 directed_dist(R->pn.xyz, R->p1.xyz, norm, sopt->norm_len);
1519 }
1520 else if ( p->nsurf > 1 )
1521 {
1522 /* get node from second surface */
1523 sp = p->surf + 1;
1524 R->pn.xyz[0] = sp->ixyz[node].x;
1525 R->pn.xyz[1] = sp->ixyz[node].y;
1526 R->pn.xyz[2] = sp->ixyz[node].z;
1527 }
1528 else
1529 R->pn = R->p1;
1530
1531 R->p1 = THD_dicomm_to_3dmm(R->dset, R->p1);
1532 R->pn = THD_dicomm_to_3dmm(R->dset, R->pn);
1533
1534 RETURN(0);
1535 }
1536
1537
1538 /*----------------------------------------------------------------------
1539 * init_range_structs
1540 *----------------------------------------------------------------------
1541 */
init_range_structs(range_3dmm * r3,range_3dmm_res * res3)1542 static int init_range_structs( range_3dmm * r3, range_3dmm_res * res3 )
1543 {
1544 ENTRY("init_range_structs");
1545
1546 r3->dset = NULL;
1547 r3->debug = 0;
1548
1549 res3->ims.num = 0;
1550 res3->ims.nall = 0;
1551 res3->ims.imarr = NULL;
1552 res3->i3arr = NULL;
1553
1554 RETURN(0);
1555 }
1556
1557
1558 /*----------------------------------------------------------------------
1559 * set_3dmm_bounds - note 3dmm bounding values
1560 *
1561 * This is an outer bounding box, like FOV, not SLAB.
1562 *----------------------------------------------------------------------
1563 */
set_3dmm_bounds(THD_3dim_dataset * dset,THD_fvec3 * min,THD_fvec3 * max)1564 int set_3dmm_bounds ( THD_3dim_dataset *dset, THD_fvec3 *min, THD_fvec3 *max)
1565 {
1566 float tmp;
1567 int c;
1568
1569 ENTRY("set_3dmm_bounds");
1570
1571 if ( !dset || !min || !max )
1572 {
1573 fprintf(stderr, "** invalid params to set_3dmm_bounds: (%p,%p,%p)\n",
1574 dset, min, max );
1575 RETURN(-1);
1576 }
1577
1578 /* get undirected bounds */
1579 min->xyz[0] = DSET_XORG(dset) - 0.5 * DSET_DX(dset);
1580 max->xyz[0] = min->xyz[0] + DSET_NX(dset) * DSET_DX(dset);
1581
1582 min->xyz[1] = DSET_YORG(dset) - 0.5 * DSET_DY(dset);
1583 max->xyz[1] = min->xyz[1] + DSET_NY(dset) * DSET_DY(dset);
1584
1585 min->xyz[2] = DSET_ZORG(dset) - 0.5 * DSET_DZ(dset);
1586 max->xyz[2] = min->xyz[2] + DSET_NZ(dset) * DSET_DZ(dset);
1587
1588 for ( c = 0; c < 3; c++ )
1589 if ( min->xyz[c] > max->xyz[c] )
1590 {
1591 tmp = min->xyz[c];
1592 min->xyz[c] = max->xyz[c];
1593 max->xyz[c] = tmp;
1594 }
1595
1596 RETURN(0);
1597 }
1598
1599
1600 /*----------------------------------------------------------------------
1601 * alloc_output_mem - for output surface dataset
1602 *----------------------------------------------------------------------
1603 */
alloc_output_mem(v2s_opts_t * sopt,v2s_param_t * p)1604 static v2s_results * alloc_output_mem(v2s_opts_t *sopt, v2s_param_t *p)
1605 {
1606 v2s_results * sd;
1607 int rv = 0, mem, nnodes;
1608
1609 ENTRY("allocate_output_mem");
1610
1611 /* if last <= 0, it will be set to nodes-1 */
1612 if ( sopt->first_node > sopt->last_node && sopt->last_node > 0 )
1613 {
1614 fprintf(stderr,"** error : first_node (%d) > last_node (%d)\n",
1615 sopt->first_node, sopt->last_node);
1616 RETURN(NULL);
1617 }
1618
1619 nnodes = p->surf[0].num_ixyz; /* just for typing ease */
1620
1621 sd = calloc(1, sizeof(v2s_results));
1622 if ( ! sd )
1623 {
1624 fprintf(stderr,"** aom: failed to allocate v2s_results struct\n");
1625 RETURN(NULL);
1626 }
1627
1628 /* explicitly initialize all pointers with NULL */
1629 sd->nodes = sd->volind = sd->i = sd->j = sd->k = sd->nvals = NULL;
1630 sd->vals = NULL;
1631 sd->labels = NULL;
1632
1633 /* verify first and last node indices */
1634 if ( sopt->first_node < 0 ) sopt->first_node = 0;
1635 if ( sopt->first_node >= nnodes ) sopt->first_node = nnodes-1;
1636 if ( sopt->last_node <= 0 ) sopt->last_node = nnodes-1;
1637 if ( sopt->last_node >= nnodes ) sopt->last_node = nnodes-1;
1638
1639 sd->nalloc = sopt->last_node - sopt->first_node + 1;
1640 sd->nused = 0;
1641
1642 /* decide the maximum number of entries per row */
1643 if ( p->over_steps ) sd->max_vals = sopt->f_steps;
1644 else if ( sopt->gp_index >= 0 ) sd->max_vals = 1;
1645 else sd->max_vals = DSET_NVALS(p->gpar);
1646
1647 if ( sopt->skip_cols & V2S_SKIP_VALS ) sd->max_vals = 1;
1648
1649 if ( p->over_steps && (sopt->skip_cols & V2S_SKIP_VALS) )
1650 {
1651 fprintf(stderr,"** if SKIP_VALS, cannot compute results over steps\n");
1652 free(sd);
1653 RETURN(NULL);
1654 }
1655
1656 /* allocate labels array - always max_vals 13 Jan 2010 */
1657 sd->nlab = sd->max_vals;
1658 sd->labels = (char **)malloc(sd->nlab * sizeof(char *));
1659 if(!sd->labels){ fprintf(stderr,"** AOM malloc failure\n"); rv = 1; }
1660
1661 sd->memory = 0;
1662
1663 /* first, compute the memory needed for one row */
1664 mem = 0;
1665 if ( !(sopt->skip_cols & V2S_SKIP_NODES ) ) mem += sizeof(int);
1666 if ( !(sopt->skip_cols & V2S_SKIP_VOLIND) ) mem += sizeof(int);
1667 if ( !(sopt->skip_cols & V2S_SKIP_I ) ) mem += sizeof(int);
1668 if ( !(sopt->skip_cols & V2S_SKIP_J ) ) mem += sizeof(int);
1669 if ( !(sopt->skip_cols & V2S_SKIP_K ) ) mem += sizeof(int);
1670 if ( !(sopt->skip_cols & V2S_SKIP_NVALS ) ) mem += sizeof(int);
1671
1672 /* now note the actual output values */
1673 if ( !(sopt->skip_cols & V2S_SKIP_VALS ) ) mem += sizeof(float);
1674 else mem += sd->max_vals * sizeof(float);
1675
1676 mem *= sd->nalloc; /* and multiply by the height of each column */
1677 sd->memory = mem;
1678
1679 if ( sopt->debug > 1 )
1680 fprintf(stderr,"+d alloc result mem: %d bytes, max_vals = %d\n",
1681 mem, sd->max_vals);
1682
1683 /* okay, this time let's allocate something... */
1684
1685 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_NODES) )
1686 rv = alloc_ints(&sd->nodes, sd->nalloc, "nodes", sopt->debug);
1687
1688 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_VOLIND) )
1689 rv = alloc_ints(&sd->volind, sd->nalloc, "volind", sopt->debug);
1690
1691 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_I) )
1692 rv = alloc_ints(&sd->i, sd->nalloc, "i", sopt->debug);
1693
1694 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_J) )
1695 rv = alloc_ints(&sd->j, sd->nalloc, "j", sopt->debug);
1696
1697 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_K) )
1698 rv = alloc_ints(&sd->k, sd->nalloc, "k", sopt->debug);
1699
1700 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_NVALS) )
1701 rv = alloc_ints(&sd->nvals, sd->nalloc, "nvals", sopt->debug);
1702
1703 if ( ! rv && ! (sopt->skip_cols & V2S_SKIP_VALS) )
1704 rv = alloc_vals_list(&sd->vals, sd->nalloc, sd->max_vals, sopt->debug);
1705 else if ( ! rv )
1706 rv = alloc_vals_list(&sd->vals, sd->nalloc, 1, sopt->debug);
1707
1708 if ( rv )
1709 {
1710 fprintf(stderr,"** failed v2s_results allocation\n");
1711 free_v2s_results(sd);
1712 sd = NULL;
1713 }
1714
1715 RETURN(sd);
1716 }
1717
1718
1719 /*----------------------------------------------------------------------
1720 * alloc_vals_list - allocate 2D array for surface data values
1721 *----------------------------------------------------------------------
1722 */
alloc_vals_list(float *** ptr,int length,int width,int debug)1723 static int alloc_vals_list(float *** ptr, int length, int width, int debug)
1724 {
1725 int c;
1726
1727 ENTRY("alloc_vals_list");
1728
1729 *ptr = (float **)malloc(width * sizeof(float *));
1730 if ( !*ptr )
1731 fprintf(stderr,"** avl: failed to alloc %d floats pointers\n", width);
1732
1733 for ( c = 0; c < width; c++ )
1734 {
1735 (*ptr)[c] = (float *)malloc(length * sizeof(float));
1736 if ( (*ptr)[c] == NULL )
1737 fprintf(stderr,"** avl: failed to alloc %d floats (# %d of %d)\n",
1738 length, c, width);
1739 }
1740
1741 if ( debug > 1 )
1742 fprintf(stderr,"-d alloc'd %d x %d floats for surf data\n",
1743 width, length);
1744
1745 RETURN(0);
1746 }
1747
1748
1749 /*----------------------------------------------------------------------
1750 * realloc_vals_list - reallocate 2D arrays for surface data values
1751 *----------------------------------------------------------------------
1752 */
realloc_vals_list(float ** ptr,int length,int width,int debug)1753 static int realloc_vals_list(float ** ptr, int length, int width, int debug)
1754 {
1755 int c, count;
1756
1757 ENTRY("realloc_vals_list");
1758
1759 count = 0;
1760 for ( c = 0; c < width; c++ )
1761 {
1762 if ( ptr[c] )
1763 {
1764 ptr[c] = (float *)realloc(ptr[c], length * sizeof(float));
1765 if ( ptr[c] == NULL )
1766 fprintf(stderr,"** rvl: fail to realloc %d floats (%d of %d)\n",
1767 length, c, width);
1768 count++;
1769 }
1770 }
1771
1772 if ( debug > 1 )
1773 fprintf(stderr,"-d realloc'd %d x %d floats for surf data\n",
1774 count, length);
1775
1776 RETURN(0);
1777 }
1778
1779
1780 /*----------------------------------------------------------------------
1781 * alloc_ints - allocate 1D array of ints
1782 *----------------------------------------------------------------------
1783 */
alloc_ints(int ** ptr,int length,char * dstr,int debug)1784 static int alloc_ints( int ** ptr, int length, char * dstr, int debug )
1785 {
1786 ENTRY("alloc_ints");
1787
1788 *ptr = (int *)malloc(length * sizeof(int));
1789 if ( ! *ptr )
1790 {
1791 fprintf(stderr,"** ai: failed to alloc %d ints for '%s'\n",length,dstr);
1792 RETURN(1);
1793 }
1794 if ( debug > 1 )
1795 fprintf(stderr,"-d ai: alloc'd %d ints for '%s'\n", length, dstr);
1796
1797 RETURN(0);
1798 }
1799
1800
1801 /*----------------------------------------------------------------------
1802 * realloc_ints - reallocate 1D array of ints (could replace alloc_ints)
1803 *----------------------------------------------------------------------
1804 */
realloc_ints(int ** ptr,int length,char * dstr,int debug)1805 static int realloc_ints( int ** ptr, int length, char * dstr, int debug )
1806 {
1807 ENTRY("realloc_ints");
1808
1809 *ptr = (int *)realloc(*ptr, length * sizeof(int));
1810 if ( ! *ptr )
1811 {
1812 fprintf(stderr,"** ri: failed to alloc %d ints for '%s'\n",length,dstr);
1813 RETURN(1);
1814 }
1815 if ( debug > 1 )
1816 fprintf(stderr,"-d ri: alloc'd %d ints for '%s'\n", length, dstr);
1817
1818 RETURN(0);
1819 }
1820
1821
1822 /*----------------------------------------------------------------------
1823 * disp_v2s_param_t - display the contents of the v2s_param_t struct
1824 *----------------------------------------------------------------------
1825 */
disp_v2s_param_t(char * info,v2s_param_t * p)1826 int disp_v2s_param_t ( char * info, v2s_param_t * p )
1827 {
1828 ENTRY("disp_v2s_param_t");
1829
1830 if ( info )
1831 fputs( info, stderr );
1832
1833 if ( p == NULL )
1834 {
1835 fprintf(stderr, "disp_v2s_param_t: p == NULL\n" );
1836 RETURN(-1);
1837 }
1838
1839 fprintf(stderr,
1840 "v2s_param_t struct at %p :\n"
1841 " gpar : vcheck = %p : %s\n"
1842 " cmask = %p\n"
1843 " nvox, over_steps = %d, %d\n"
1844 " nsurf = %d\n"
1845 , p,
1846 p->gpar, ISVALID_DSET(p->gpar) ? "valid" : "invalid",
1847 p->cmask, p->nvox, p->over_steps, p->nsurf
1848 );
1849
1850 RETURN(0);
1851 }
1852
1853
1854 /*----------------------------------------------------------------------
1855 * disp_v2s_opts_t - display the contents of the v2s_opts_t struct
1856 *----------------------------------------------------------------------
1857 */
disp_v2s_opts_t(char * info,v2s_opts_t * sopt)1858 int disp_v2s_opts_t ( char * info, v2s_opts_t * sopt )
1859 {
1860 ENTRY("disp_v2s_opts_t");
1861
1862 if ( info )
1863 fputs( info, stderr );
1864
1865 if ( sopt == NULL )
1866 {
1867 fprintf(stderr, "disp_v2s_opts_t: sopt == NULL\n");
1868 RETURN(-1);
1869 }
1870
1871 fprintf(stderr,
1872 "v2s_opts_t struct at %p :\n"
1873 " map, gp_index = %d, %d\n"
1874 " debug, dnode = %d, %d\n"
1875 " no_head, skip_cols = %d, %d\n"
1876 " first_node, last_node = %d, %d\n"
1877 " use_norms, norm_len = %d, %f\n"
1878 " norm_dir = %d\n"
1879 " f_index, f_steps = %d, %d\n"
1880 " f_p1_fr, f_pn_fr = %f, %f\n"
1881 " f_p1_mm, f_pn_mm = %f, %f\n"
1882 " outfile_1D = %s\n"
1883 " outfile_niml = %s\n"
1884 " segc_file = %s\n"
1885 " fake, argc, argv = %d, %d, %p\n"
1886 , sopt,
1887 sopt->map, sopt->gp_index, sopt->debug, sopt->dnode,
1888 sopt->no_head, sopt->skip_cols,
1889 sopt->first_node, sopt->last_node,
1890 sopt->use_norms, sopt->norm_len, sopt->norm_dir,
1891 sopt->f_index, sopt->f_steps,
1892 sopt->f_p1_fr, sopt->f_pn_fr, sopt->f_p1_mm, sopt->f_pn_mm,
1893 CHECK_NULL_STR(sopt->outfile_1D),
1894 CHECK_NULL_STR(sopt->outfile_niml),
1895 CHECK_NULL_STR(sopt->segc_file),
1896 sopt->cmd.fake, sopt->cmd.argc, sopt->cmd.argv
1897 );
1898
1899 RETURN(0);
1900 }
1901
1902
1903 /*----------------------------------------------------------------------
1904 * magnitude_f - return magnitude of float vector
1905 *----------------------------------------------------------------------
1906 */
magnitude_f(float * p,int length)1907 static double magnitude_f( float * p, int length )
1908 {
1909 int c;
1910 double sums = 0.0;
1911
1912 for ( c = 0; c < length; c++, p++ )
1913 sums += (*p) * (*p);
1914
1915 return(sqrt(sums));
1916 }
1917
1918
1919 /*----------------------------------------------------------------------
1920 * dist_f3mm - return Euclidean distance between the points
1921 *----------------------------------------------------------------------
1922 */
dist_f3mm(THD_fvec3 * p1,THD_fvec3 * p2)1923 static float dist_f3mm( THD_fvec3 * p1, THD_fvec3 * p2 )
1924 {
1925 double d0, d1, d2;
1926
1927 if ( p1 == NULL || p2 == NULL )
1928 {
1929 fprintf( stderr, "** dist_f3mm: invalid params (%p,%p)\n", p1, p2 );
1930 return(0.0);
1931 }
1932
1933 d0 = p1->xyz[0] - p2->xyz[0];
1934 d1 = p1->xyz[1] - p2->xyz[1];
1935 d2 = p1->xyz[2] - p2->xyz[2];
1936
1937 return(sqrt(d0*d0 + d1*d1 + d2*d2));
1938 }
1939
1940
1941 /*---------------------------------------------------------------------------
1942 * disp_range_3dmm - display the contents of the range_3dmm struct
1943 *---------------------------------------------------------------------------
1944 */
disp_range_3dmm(char * info,range_3dmm * dp)1945 static int disp_range_3dmm ( char * info, range_3dmm * dp )
1946 {
1947 ENTRY("disp_range_3dmm");
1948
1949 if ( info )
1950 fputs( info, stderr );
1951
1952 if ( dp == NULL )
1953 {
1954 fprintf(stderr, "disp_range_3dmm: dp == NULL\n");
1955 RETURN(-1);
1956 }
1957
1958 fprintf(stderr,
1959 "range_3dmm struct at %p :\n"
1960 " dset = %p : %s\n"
1961 " p1 = (%f, %f, %f)\n"
1962 " pn = (%f, %f, %f)\n"
1963 " dset_min = (%f, %f, %f)\n"
1964 " dset_max = (%f, %f, %f)\n"
1965 " oob_check, debug = (%d, %d)\n",
1966 dp, dp->dset, ISVALID_DSET(dp->dset) ? "valid" : "invalid",
1967 dp->p1.xyz[0], dp->p1.xyz[1], dp->p1.xyz[2],
1968 dp->pn.xyz[0], dp->pn.xyz[1], dp->pn.xyz[2],
1969 dp->dset_min.xyz[0], dp->dset_min.xyz[1], dp->dset_min.xyz[2],
1970 dp->dset_max.xyz[0], dp->dset_max.xyz[1], dp->dset_max.xyz[2],
1971 dp->oob_check, dp->debug
1972 );
1973
1974 RETURN(0);
1975 }
1976
1977
1978 /*---------------------------------------------------------------------------
1979 * disp_range_3dmm_res - display the contents of the range_3dmm_res struct
1980 *---------------------------------------------------------------------------
1981 */
disp_range_3dmm_res(char * info,range_3dmm_res * dp)1982 static int disp_range_3dmm_res ( char * info, range_3dmm_res * dp )
1983 {
1984 ENTRY("disp_range_3dmm_res");
1985
1986 if ( info )
1987 fputs( info, stderr );
1988
1989 if ( dp == NULL )
1990 {
1991 fprintf(stderr, "disp_range_3dmm: dp == NULL\n");
1992 RETURN(-1);
1993 }
1994
1995 fprintf(stderr,
1996 "range_3dmm_res struct at %p :\n"
1997 " ims.imarr = %p\n"
1998 " ims.num, ims.nall = %d, %d\n"
1999 " repeats, masked = %d, %d\n"
2000 " oob, ifirst = %d, %d\n"
2001 " i3first[0,1,2] = %d, %d, %d\n"
2002 " i3arr = %p\n"
2003 , dp,
2004 dp->ims.imarr, dp->ims.num, dp->ims.nall,
2005 dp->repeats, dp->masked, dp->oob, dp->ifirst,
2006 dp->i3first.ijk[0], dp->i3first.ijk[1], dp->i3first.ijk[2],
2007 dp->i3arr );
2008
2009 if ( dp->i3arr )
2010 fprintf(stderr,
2011 " i3arr[0].ijk = %d, %d, %d\n",
2012 dp->i3arr[0].ijk[0], dp->i3arr[0].ijk[1], dp->i3arr[0].ijk[2] );
2013
2014 RETURN(0);
2015 }
2016
2017
2018 /*---------------------------------------------------------------------------
2019 * disp_mri_imarr - display the contents of the MRI_IMARR struct
2020 *---------------------------------------------------------------------------
2021 */
disp_mri_imarr(char * info,MRI_IMARR * dp)2022 int disp_mri_imarr ( char * info, MRI_IMARR * dp )
2023 {
2024 float * fp;
2025 int cr, cc;
2026
2027 ENTRY("disp_mri_imarr");
2028
2029 if ( info )
2030 fputs( info, stderr );
2031
2032 if ( dp == NULL )
2033 {
2034 fprintf(stderr, "disp_mri_imarr: dp == NULL\n");
2035 RETURN(-1);
2036 }
2037
2038 fprintf(stderr,
2039 "mri_imarr struct at %p :\n"
2040 " num, nall = %d, %d\n",
2041 dp, dp->num, dp->nall );
2042
2043 for ( cr = 0; cr < dp->num; cr++ )
2044 {
2045 fp = MRI_FLOAT_PTR(dp->imarr[cr]);
2046 fprintf(stderr, " %3d: ", cr);
2047 for ( cc = 0; cc < dp->imarr[cr]->nx; cc++, fp++ )
2048 fprintf(stderr, "%f ", *fp );
2049 fputc( '\n', stderr );
2050 }
2051
2052 RETURN(0);
2053 }
2054
2055
2056 /***********************************************************************
2057 * disp_surf_vals
2058 ***********************************************************************
2059 */
disp_surf_vals(char * mesg,v2s_results * sd,int node)2060 static int disp_surf_vals( char * mesg, v2s_results * sd, int node )
2061 {
2062 int index, c;
2063
2064 ENTRY("disp_surf_vals");
2065
2066 fprintf(stderr, "-------------------------------------------------\n");
2067 if ( mesg ) fputs( mesg, stderr );
2068 if ( sd->nused < 1 )
2069 {
2070 fprintf(stderr,"** no surf nodes defined\n");
2071 RETURN(-1);
2072 }
2073
2074 index = (node >= 0) ? node : sd->nused - 1;
2075
2076 fprintf(stderr, "v2s_results vals for sd_index %d, node %d :\n"
2077 " volind, (i, j, k) = %d, (%d, %d, %d)\n"
2078 " nvals: values... = %d: ", index,
2079 sd->nodes ? sd->nodes[index] : 0,
2080 sd->volind ? sd->volind[index] : 0,
2081 sd->i ? sd->i[index] : 0,
2082 sd->j ? sd->j[index] : 0,
2083 sd->k ? sd->k[index] : 0,
2084 sd->nvals ? sd->nvals[index] : 0 );
2085
2086 for ( c = 0; c < sd->max_vals; c++ )
2087 fprintf(stderr,"%s ", MV_format_fval(sd->vals[c][index]));
2088 fputc('\n', stderr);
2089
2090 RETURN(0);
2091 }
2092
2093
2094 /***********************************************************************
2095 * disp_v2s_results
2096 ***********************************************************************
2097 */
disp_v2s_results(char * mesg,v2s_results * d)2098 int disp_v2s_results( char * mesg, v2s_results * d )
2099 {
2100 ENTRY("disp_v2s_results");
2101
2102 if ( mesg ) fputs( mesg, stderr );
2103
2104 fprintf(stderr, "v2s_results @ %p\n"
2105 " nalloc, nused = %d, %d\n"
2106 " max_vals, memory = %d, %d\n"
2107 " nodes, volind = %p, %p\n"
2108 " i, j, k = %p, %p, %p\n"
2109 " nvals, vals = %p, %p\n"
2110 " labels, nlab = %p, %d\n",
2111 d, d->nalloc, d->nused, d->max_vals, d->memory,
2112 d->nodes, d->volind, d->i, d->j, d->k, d->nvals, d->vals,
2113 d->labels, d->nlab);
2114
2115 RETURN(0);
2116 }
2117
2118
2119 /***********************************************************************
2120 * disp_v2s_plugin_opts
2121 ***********************************************************************
2122 */
disp_v2s_plugin_opts(char * mesg,v2s_plugin_opts * d)2123 int disp_v2s_plugin_opts( char * mesg, v2s_plugin_opts * d )
2124 {
2125 ENTRY("disp_v2s_plugin_opts");
2126
2127 if ( mesg ) fputs( mesg, stderr );
2128
2129 fprintf(stderr, "v2s_plugin_opts @ %p\n"
2130 " ready = %d\n"
2131 " map_all = %d\n"
2132 " use0, use1 = %d, %d\n"
2133 " s0A, s0B = %d, %d\n"
2134 " s1A, s1B = %d, %d\n"
2135 " gpt_index = %d\n"
2136 " gpt_thresh = %f\n"
2137 " label[0,1] = %s, %s\n"
2138 " label[2,3] = %s, %s\n"
2139 " surf_vol = %s\n",
2140 d, d->ready, d->map_all, d->use0, d->use1,
2141 d->s0A, d->s0B, d->s1A, d->s1B,
2142 d->gpt_index, d->gpt_thresh,
2143 CHECK_NULL_STR(d->label[0]),
2144 CHECK_NULL_STR(d->label[1]),
2145 CHECK_NULL_STR(d->label[2]),
2146 CHECK_NULL_STR(d->label[3]),
2147 d->sv_dset ? DSET_FILECODE(d->sv_dset) : "NULL"
2148 );
2149 RETURN(0);
2150 }
2151
2152
2153 /*----------------------------------------------------------------------
2154 * free_v2s_results - free contents of v2s_results struct
2155 *----------------------------------------------------------------------
2156 */
free_v2s_results(v2s_results * sd)2157 int free_v2s_results( v2s_results * sd )
2158 {
2159 int c;
2160
2161 ENTRY("free_v2s_results");
2162
2163 if( ! sd ) RETURN(0);
2164
2165 if (sd->nodes) { free(sd->nodes); sd->nodes = NULL; }
2166 if (sd->volind) { free(sd->volind); sd->volind = NULL; }
2167 if (sd->i) { free(sd->i); sd->i = NULL; }
2168 if (sd->j) { free(sd->j); sd->j = NULL; }
2169 if (sd->k) { free(sd->k); sd->k = NULL; }
2170 if (sd->nvals) { free(sd->nvals); sd->nvals = NULL; }
2171
2172 if (sd->vals)
2173 {
2174 for ( c = 0; c < sd->max_vals; c++ )
2175 if ( sd->vals[c] ) { free(sd->vals[c]); sd->vals[c] = NULL; }
2176
2177 free(sd->vals);
2178 sd->vals = NULL;
2179 }
2180
2181 if (sd->labels && sd->nlab > 0)
2182 {
2183 for ( c = 0; c < sd->nlab; c++ )
2184 if ( sd->labels[c] ) { free(sd->labels[c]); sd->labels[c] = NULL; }
2185
2186 free(sd->labels);
2187 sd->labels = NULL;
2188 }
2189
2190 free(sd);
2191
2192 RETURN(0);
2193 }
2194
2195
2196 /*----------------------------------------------------------------------
2197 * validate structure contents
2198 *----------------------------------------------------------------------
2199 */
validate_v2s_inputs(v2s_opts_t * sopt,v2s_param_t * p)2200 static int validate_v2s_inputs ( v2s_opts_t * sopt, v2s_param_t * p )
2201 {
2202 int c;
2203
2204 ENTRY("validate_v2s_inputs");
2205
2206 if ( !sopt || !p )
2207 {
2208 fprintf(stderr,"** validate_v2s_inputs: bad params (%p,%p)\n", sopt, p);
2209 RETURN(-1);
2210 }
2211
2212 /* validate surface options structure */
2213 if ( sopt->map <= E_SMAP_INVALID || sopt->map >= E_SMAP_FINAL )
2214 {
2215 fprintf(stderr,"** invalid mapping index %d\n", sopt->map);
2216 RETURN(1);
2217 }
2218 else if ( v2s_map_type(gv2s_map_names[sopt->map]) != sopt->map )
2219 {
2220 fprintf(stderr,"** mapping index failure for %d\n", sopt->map);
2221 RETURN(1);
2222 }
2223
2224 if ( sopt->f_steps <= 0 || sopt->f_steps >= V2S_STEPS_TOOOOO_BIG )
2225 {
2226 fprintf(stderr,"** invalid f_steps = %d\n", sopt->f_steps);
2227 RETURN(1);
2228 }
2229
2230 /* validate the contents of p, proper */
2231 if ( !p->gpar ) {fprintf(stderr,"** vv2si: no dset?\n"); RETURN(2);}
2232
2233 if ( p->nvox != DSET_NVOX(p->gpar) )
2234 {
2235 fprintf(stderr,"** invalid voxel count (%d) for grid_parent\n",p->nvox);
2236 RETURN(2);
2237 }
2238
2239 if ( sopt->gp_index >= DSET_NVALS(p->gpar) )
2240 {
2241 fprintf(stderr,"** gp_index (%d) > max grid_parent index (%d)\n",
2242 sopt->gp_index, DSET_NVALS(p->gpar) - 1);
2243 RETURN(2);
2244 }
2245
2246 if ( p->nsurf < 1 || p->nsurf > 2 ) /* see V2S_MAX_SURFS */
2247 {
2248 fprintf(stderr,"** invalid: nsurf = %d must be in [%d,%d]\n",
2249 p->nsurf, 1, 2);
2250 RETURN(2);
2251 }
2252
2253 /* validate individual SUMA_surface structs */
2254 for ( c = 0; c < p->nsurf; c++ )
2255 if ( check_SUMA_surface(p->surf + c) )
2256 RETURN(3);
2257
2258 /* now look for consistency */
2259 if ( p->nsurf > 1 )
2260 {
2261 if ( p->surf[0].num_ixyz != p->surf[1].num_ixyz )
2262 {
2263 fprintf(stderr,"** invalid surfaces, different # nodes (%d,%d)\n",
2264 p->surf[0].num_ixyz, p->surf[1].num_ixyz);
2265 RETURN(4);
2266 }
2267 }
2268 else if ( sopt->use_norms && !p->surf[0].norm )
2269 {
2270 fprintf(stderr,"** missing surface normals for surface '%s'\n",
2271 CHECK_EMPTY_STR(p->surf[0].label));
2272 RETURN(4);
2273 }
2274
2275 RETURN(0);
2276 }
2277
2278 /* return 0 if valid, > 0 otherwise */
check_SUMA_surface(SUMA_surface * s)2279 static int check_SUMA_surface( SUMA_surface * s )
2280 {
2281 int rv = 0;
2282 ENTRY("is_valid_SUMA_surface");
2283
2284 if ( !s ) { fprintf(stderr,"** ivSs: no surface\n"); RETURN(0); }
2285
2286 if ( s->type != SUMA_SURFACE_TYPE )
2287 {
2288 fprintf(stderr,"** surf '%s': invalid type %d\n",
2289 CHECK_EMPTY_STR(s->label), s->type);
2290 rv++;
2291 }
2292
2293 if ( s->num_ixyz < 0 || s->nall_ixyz < 0 || s->num_ixyz < s->nall_ixyz )
2294 {
2295 fprintf(stderr,"** surf '%s': invalid num_ixyz = %d, nall_ixyz = %d\n",
2296 CHECK_EMPTY_STR(s->label), s->num_ixyz, s->nall_ixyz);
2297 rv++;
2298 }
2299
2300 if ( s->seq == 0 || s->sorted == 0 || s->seqbase != 0 )
2301 {
2302 fprintf(stderr,"** surf '%s': invalid seq, sort or base (%d, %d, %d)\n",
2303 CHECK_EMPTY_STR(s->label), s->seq, s->sorted, s->seqbase);
2304 rv++;
2305 }
2306
2307 if ( !s->ixyz )
2308 {
2309 fprintf(stderr,"** surf '%s': invalid, missing nodes\n",
2310 CHECK_EMPTY_STR(s->label));
2311 rv++;
2312 }
2313
2314 RETURN(rv);
2315 }
2316
2317
2318 /*---------------------------------------------------------------------------
2319 * v2s_vals_over_steps - return whether a function is displayed over steps
2320 *
2321 * Most function results are output per sub-brick. These functions will
2322 * have results displayed over the segment steps.
2323 *---------------------------------------------------------------------------
2324 */
v2s_vals_over_steps(int map)2325 int v2s_vals_over_steps( int map )
2326 {
2327 if ( map == E_SMAP_SEG_VALS )
2328 return 1;
2329
2330 return 0;
2331 }
2332
2333
2334 /*----------------------------------------------------------------------
2335 * v2s_map_type - return an E_SMAP_XXX code
2336 *
2337 * on failure, return -1 (E_SMAP_INVALID)
2338 * else return >0 (a valid map code)
2339 *----------------------------------------------------------------------
2340 */
v2s_map_type(char * map_str)2341 int v2s_map_type ( char * map_str )
2342 {
2343 v2s_map_nums map;
2344
2345 ENTRY("v2s_map_type");
2346
2347 if ( map_str == NULL )
2348 {
2349 fprintf( stderr, "** v2s_map_type: missing map_str parameter\n" );
2350 RETURN((int)E_SMAP_INVALID);
2351 }
2352
2353 if ( sizeof(gv2s_map_names) / sizeof(char *) != (int)E_SMAP_FINAL )
2354 {
2355 fprintf( stderr, "** error: gv2s_map_names/v2s_map_num mismatch\n");
2356 RETURN((int)E_SMAP_INVALID);
2357 }
2358
2359 /* not ready for E_SMAP_COUNT yet (until someone wants it) */
2360 if ( !strcmp( map_str, gv2s_map_names[E_SMAP_COUNT] ) )
2361 RETURN((int)E_SMAP_INVALID);
2362
2363 for ( map = E_SMAP_INVALID; map < E_SMAP_FINAL; map++ )
2364 if ( !strcmp( map_str, gv2s_map_names[map] ) )
2365 RETURN((int)map);
2366
2367 RETURN((int)E_SMAP_INVALID);
2368 }
2369
2370
2371 /*----------------------------------------------------------------------
2372 * check for a map index that we consider valid
2373 *
2374 * from anywhere, E_SMAP_MASK2 and E_SMAP_COUNT are not yet implemented
2375 * from afni, E_SMAP_SEG_VALS is not acceptable (only allow 1 output)
2376 *----------------------------------------------------------------------
2377 */
v2s_is_good_map(int map,int from_afni)2378 int v2s_is_good_map( int map, int from_afni )
2379 {
2380 ENTRY("v2s_good_map_index");
2381
2382 if ( map <= E_SMAP_INVALID || map >= E_SMAP_FINAL )
2383 RETURN(0);
2384
2385 /* these have not (yet? do we care?) been implemented */
2386 if ( map == E_SMAP_MASK2 || map == E_SMAP_COUNT )
2387 RETURN(0);
2388
2389 if ( from_afni && map == E_SMAP_SEG_VALS )
2390 RETURN(0);
2391
2392 RETURN(1);
2393 }
2394
2395 /*----------------------------------------------------------------------
2396 * check for a map index that we consider valid
2397 *
2398 * from anywhere, E_SMAP_MASK2 and E_SMAP_COUNT are not yet implemented
2399 * from afni, E_SMAP_SEG_VALS is not acceptable (only allow 1 output)
2400 *----------------------------------------------------------------------
2401 */
v2s_fill_sopt_default(v2s_opts_t * sopt,int nsurf)2402 int v2s_fill_sopt_default(v2s_opts_t * sopt, int nsurf )
2403 {
2404
2405 ENTRY("v2s_fill_sopt_default");
2406
2407 if ( !sopt || nsurf < 1 || nsurf > 2 )
2408 {
2409 fprintf(stderr,"** FSAD: bad params (%p,%d)\n",sopt,nsurf);
2410 RETURN(1);
2411 }
2412
2413 /* first set any zeros */
2414 memset(sopt, 0, sizeof(*sopt));
2415
2416 if ( nsurf == 2 )
2417 sopt->map = E_SMAP_MIDPT;
2418 else
2419 sopt->map = E_SMAP_MASK;
2420
2421 sopt->gp_index = -1;
2422 sopt->no_head = 1;
2423 sopt->skip_cols = V2S_SKIP_ALL ^ V2S_SKIP_NODES; /* nodes and 1 val */
2424 sopt->f_steps = 1;
2425 sopt->outfile_1D = NULL;
2426 sopt->outfile_niml = NULL;
2427 sopt->segc_file = NULL;
2428
2429 sopt->cmd.fake = 0;
2430 sopt->cmd.argc = 0;
2431 sopt->cmd.argv = NULL;
2432
2433 RETURN(0);
2434 }
2435
2436
2437 /*----------------------------------------------------------------------
2438 * v2s_write_outfile_1D - write results to 1D file
2439 *----------------------------------------------------------------------
2440 */
v2s_write_outfile_1D(v2s_opts_t * sopt,v2s_results * sd,char * label)2441 int v2s_write_outfile_1D ( v2s_opts_t * sopt, v2s_results * sd, char * label )
2442 {
2443 FILE * fp;
2444 int c, c2;
2445 ENTRY("v2s_write_outfile_1D");
2446
2447 fp = fopen( sopt->outfile_1D, "w" );
2448 if ( fp == NULL )
2449 {
2450 fprintf( stderr, "** failure to open '%s' for writing\n",
2451 sopt->outfile_1D );
2452 RETURN(-1);
2453 }
2454
2455 if ( ! sopt->no_head )
2456 print_header(fp, label, gv2s_map_names[sopt->map], sd);
2457
2458 for ( c = 0; c < sd->nused; c++ )
2459 {
2460 /* keep old spacing */
2461 fputc(' ', fp);
2462 if ( sd->nodes ) fprintf(fp, " %8d", sd->nodes[c]);
2463 if ( sd->volind ) fprintf(fp, " %8d ", sd->volind[c]);
2464 if ( sd->i ) fprintf(fp, " %3d", sd->i[c]);
2465 if ( sd->j ) fprintf(fp, " %3d", sd->j[c]);
2466 if ( sd->k ) fprintf(fp, " %3d", sd->k[c]);
2467 if ( sd->nvals ) fprintf(fp, " %3d", sd->nvals[c]);
2468
2469 for ( c2 = 0; c2 < sd->max_vals; c2++ )
2470 fprintf(fp, " %10s", MV_format_fval(sd->vals[c2][c]));
2471 fputc('\n', fp);
2472 }
2473
2474 fclose(fp);
2475
2476 RETURN(0);
2477 }
2478
2479
2480 /*----------------------------------------------------------------------
2481 * v2s_write_outfile_niml - write results to niml file
2482 * - free data pointers as we go
2483 *----------------------------------------------------------------------
2484 */
v2s_write_outfile_niml(v2s_opts_t * sopt,v2s_results * sd,int free_vals)2485 int v2s_write_outfile_niml ( v2s_opts_t * sopt, v2s_results * sd, int free_vals)
2486 {
2487 static char v2s_name[] = "3dVol2Surf_dataset";
2488 NI_element * nel = NULL;
2489 NI_stream ns;
2490 char * ni_name;
2491 int c;
2492
2493 ENTRY("v2s_write_outfile_niml");
2494
2495 if ( !sopt->outfile_niml ) RETURN(0);
2496
2497 nel = NI_new_data_element( v2s_name, sd->nused );
2498 if ( !nel )
2499 {
2500 fprintf(stderr,"** file NI_new_data_element, n = '%s', len = %d\n",
2501 v2s_name, sd->nused);
2502 RETURN(1);
2503 }
2504
2505 ni_name = (char *)calloc(strlen(sopt->outfile_niml)+6, sizeof(char));
2506 if ( !ni_name ) { fprintf(stderr,"** ni_name failed\n"); RETURN(1); }
2507 sprintf(ni_name, "file:%s", sopt->outfile_niml);
2508
2509 ns = NI_stream_open(ni_name, "w");
2510
2511 NI_add_column(nel,NI_INT,sd->nodes);
2512
2513 for ( c = 0; c < sd->max_vals; c++ )
2514 {
2515 NI_add_column(nel, NI_FLOAT, sd->vals[c]);
2516 if ( free_vals ) { free(sd->vals[c]); sd->vals[c] = NULL; }
2517 }
2518 if ( free_vals ) { free(sd->vals); sd->vals = NULL; }
2519
2520 if ( NI_write_element(ns, nel, NI_BINARY_MODE) < 0 )
2521 {
2522 fprintf(stderr,"** NI_write_element failed for: '%s'\n", ni_name);
2523 RETURN(1);
2524 }
2525
2526 NI_free_element( nel );
2527 NI_stream_close( ns );
2528 free(ni_name);
2529
2530 RETURN(0);
2531 }
2532
2533
2534 /*----------------------------------------------------------------------
2535 * v2s_write_outfile_NSD - write results to NI_SURF_DSET file
2536 *----------------------------------------------------------------------
2537 */
v2s_write_outfile_NSD(v2s_results * sd,v2s_opts_t * sopt,v2s_param_t * p,int free_vals)2538 int v2s_write_outfile_NSD(v2s_results *sd, v2s_opts_t * sopt,
2539 v2s_param_t * p, int free_vals)
2540 {
2541 SUMA_DSET * sdset;
2542 NI_element * nel;
2543 void ** elist = NULL;
2544 char * oname, * statsym;
2545 int c, rv, ind;
2546
2547 ENTRY("v2s_write_outfile_NSD");
2548
2549 if ( !sd || !p || !sopt || !sopt->outfile_niml ) RETURN(1);
2550
2551 /* create an empty dataset without an idcode or domain string */
2552 sdset = SUMA_CreateDsetPointer(sopt->outfile_niml, SUMA_NODE_BUCKET,
2553 NULL, NULL, sd->nused);
2554 /* add node indices, if they exist */
2555 if( sd->nodes )
2556 {
2557 rv = SUMA_AddDsetNelCol(sdset, "Node Indices", SUMA_NODE_INDEX,
2558 (void *)sd->nodes, NULL, 1);
2559 if( !rv ){ fprintf(stderr,"** WO_NSD add nodes failure\n"); RETURN(1); }
2560
2561 if( free_vals ){ free(sd->nodes); sd->nodes = NULL; }
2562 if( sopt->debug>1 ) fprintf(stderr,"+d adding node indices to NSD\n");
2563 }
2564
2565 for( c = 0; c < sd->max_vals; c++ )
2566 {
2567 rv = SUMA_AddDsetNelCol(sdset, sd->labels[c],
2568 SUMA_NODE_FLOAT, sd->vals[c],NULL,1);
2569 if( !rv ){ fprintf(stderr,"** WO_NSD add col failure\n"); RETURN(1); }
2570 if( free_vals ){ free(sd->vals[c]); sd->vals[c] = NULL; }
2571 }
2572
2573 /* add history (may need to fake it) */
2574 SUMA_AddNgrHist(sdset->ngr,
2575 sopt->cmd.fake ? "Vol2Surf_plugin" : "3dVol2Surf",
2576 sopt->cmd.argc, sopt->cmd.argv);
2577
2578 set_ni_globs_from_env(); /* init niml globals from environment */
2579 set_gni_debug(sopt->debug);
2580
2581 /*--- add COLMS_STATSYM (probably already there) ---*/
2582 if( sopt->gp_index >= 0 ) ind = sopt->gp_index; /* pick one sub-brick */
2583 else if( p->over_steps ) ind = 0; /* must be sub-brick 0 */
2584 else ind = -1; /* get all sub-bricks */
2585 statsym = THD_make_statsym_string(p->gpar, ind);
2586 if( statsym ) {
2587 nel = SUMA_FindDsetAttributeElement(sdset, "COLMS_STATSYM");
2588 if( nel ) { /* if it exists, replace the old, else make one */
2589 SUMA_NEL_REPLACE_STRING(nel, 0, 0, statsym);
2590 if( sopt->debug > 2 )
2591 fprintf(stderr,"++ replaced COLMS_STATSYM with '%s'\n",statsym);
2592 } else {
2593 nel = NI_new_data_element("AFNI_atr", 1);
2594 NI_set_attribute(nel,"atr_name", "COLMS_STATSYM");
2595 NI_add_column(nel, NI_STRING, &statsym);
2596 NI_add_to_group(sdset->ngr, nel);
2597 if( sopt->debug > 2 )
2598 fprintf(stderr,"++ added COLMS_STATSYM as '%s'\n", statsym);
2599 }
2600 free(statsym);
2601 } else if ( sopt->debug > 1 )
2602 fprintf(stderr,"** failed to make statsym_string...\n");
2603
2604 /* find the data element and set the output format */
2605 c = NI_search_group_shallow(sdset->ngr, "SPARSE_DATA", &elist);
2606 if( c == 1 && (nel = (NI_element *)elist[0]) )
2607 { NI_free(elist); set_sparse_data_attribs(nel, p->gpar, 0); }
2608 else
2609 fprintf(stderr, "** WO_NSD: missing SPARSE_DATA?\n");
2610
2611 oname = SUMA_WriteDset_ns(sopt->outfile_niml, sdset,
2612 SUMA_NO_DSET_FORMAT, 1,1);
2613 if(sopt->debug && oname) fprintf(stderr,"+d wrote NI_SURF_DSET %s\n",oname);
2614
2615 SUMA_FreeDset(sdset);
2616
2617 RETURN(0);
2618 }
2619
2620
2621 /*----------------------------------------------------------------------
2622 * print_header - dump standard header for node output - v2.4
2623 *----------------------------------------------------------------------
2624 */
print_header(FILE * outfp,char * surf,char * map,v2s_results * sd)2625 static int print_header(FILE * outfp, char * surf, char * map, v2s_results * sd)
2626 {
2627 int val;
2628
2629 ENTRY("print_header");
2630
2631 fprintf( outfp, "# --------------------------------------------------\n" );
2632 fprintf( outfp, "# surface '%s', '%s' :\n", surf, map );
2633 fprintf( outfp, "#\n" );
2634
2635 /* keep old style, but don't presume all columns get used (v 6.0) :
2636 * fprintf( outfp, "# node 1dindex i j k vals" );
2637 * fprintf( outfp, "# ------ ------- --- --- --- ----" );
2638 */
2639
2640 /* output column headers */
2641 fputc( '#', outfp ); /* still comment line */
2642
2643 if ( sd->nodes ) fprintf(outfp, " node ");
2644 if ( sd->volind ) fprintf(outfp, " 1dindex ");
2645 if ( sd->i ) fprintf(outfp, " i ");
2646 if ( sd->j ) fprintf(outfp, " j ");
2647 if ( sd->k ) fprintf(outfp, " k ");
2648 if ( sd->nvals ) fprintf(outfp, " vals");
2649
2650 for ( val = 0; val < sd->max_vals; val++ )
2651 fprintf( outfp, " v%-2d ", val );
2652 fputc( '\n', outfp );
2653
2654 fputc( '#', outfp );
2655 /* underline the column headers */
2656 if ( sd->nodes ) fprintf(outfp, " ------");
2657 if ( sd->volind ) fprintf(outfp, " ------- ");
2658 if ( sd->i ) fprintf(outfp, " ---");
2659 if ( sd->j ) fprintf(outfp, " ---");
2660 if ( sd->k ) fprintf(outfp, " ---");
2661 if ( sd->nvals ) fprintf(outfp, " ----");
2662
2663 fputs( " ", outfp );
2664 for ( val = 0; val < sd->max_vals; val++ )
2665 fprintf( outfp, " -------- " );
2666 fputc( '\n', outfp );
2667
2668 RETURN(0);
2669 }
2670
v2s_free_cmd(v2s_opts_t * sopt)2671 static int v2s_free_cmd(v2s_opts_t * sopt)
2672 {
2673 int c;
2674
2675 ENTRY("v2s_free_cmd");
2676
2677 if( ! sopt->cmd.fake ) RETURN(0);
2678 if( sopt->cmd.argc <= 0 || !sopt->cmd.argv ) RETURN(0);
2679
2680 for( c = 0; c < sopt->cmd.argc; c++ )
2681 if( sopt->cmd.argv[c] ) free(sopt->cmd.argv[c]);
2682 free(sopt->cmd.argv);
2683
2684 sopt->cmd.fake = 0;
2685 sopt->cmd.argc = 0;
2686 sopt->cmd.argv = NULL;
2687
2688 RETURN(0);
2689 }
2690
2691 /*----------------------------------------------------------------------
2692 * create a command string from the passed structures 9 Aug 2006 [rickr]
2693 *
2694 * allocate and fill a string that might work as a command
2695 *----------------------------------------------------------------------
2696 */
v2s_make_command(v2s_opts_t * opt,v2s_param_t * p)2697 int v2s_make_command( v2s_opts_t * opt, v2s_param_t * p )
2698 {
2699 char ** argv = NULL, str[2048];
2700 char * dset_file = NULL, * sv_file = NULL;
2701 int argc = 0, acnall = 0;
2702
2703 ENTRY("v2s_make_command");
2704
2705 /* set the sv and grid_parent filenames */
2706 if(gv2s_plug_opts.sv_dset) sv_file = DSET_FILECODE(gv2s_plug_opts.sv_dset);
2707 else sv_file = "UNKNOWN_SURF_VOL";
2708 dset_file = DSET_FILECODE(p->gpar);
2709
2710 /* start setting options (3dVol2Surf may get replaced) */
2711 loc_add_2_list(&argv, &acnall, &argc, "3dVol2Surf");
2712
2713 loc_add_2_list(&argv, &acnall, &argc, "-spec");
2714 if( p->surf[0].spec_file[0] )
2715 loc_add_2_list(&argv, &acnall, &argc, p->surf[0].spec_file);
2716 else
2717 loc_add_2_list(&argv, &acnall, &argc, "NO_SPEC_FILE");
2718
2719 loc_add_2_list(&argv, &acnall, &argc, "-surf_A");
2720 loc_add_2_list(&argv, &acnall, &argc, p->surf[0].label);
2721 if( p->nsurf == 2 ){
2722 loc_add_2_list(&argv, &acnall, &argc, "-surf_B");
2723 loc_add_2_list(&argv, &acnall, &argc, p->surf[1].label);
2724 }
2725
2726 loc_add_2_list(&argv, &acnall, &argc, "-sv");
2727 loc_add_2_list(&argv, &acnall, &argc, sv_file);
2728
2729 loc_add_2_list(&argv, &acnall, &argc, "-grid_parent");
2730 loc_add_2_list(&argv, &acnall, &argc, dset_file);
2731 loc_add_2_list(&argv, &acnall, &argc, "-gp_index");
2732 sprintf(str,"%d",opt->gp_index);
2733 loc_add_2_list(&argv, &acnall, &argc, str);
2734
2735 loc_add_2_list(&argv, &acnall, &argc, "-map_func");
2736 loc_add_2_list(&argv, &acnall, &argc, gv2s_map_names[opt->map]);
2737
2738 sprintf(str, "%d", opt->f_steps);
2739 loc_add_2_list(&argv, &acnall, &argc, "-f_steps");
2740 loc_add_2_list(&argv, &acnall, &argc, str);
2741
2742 loc_add_2_list(&argv, &acnall, &argc, "-f_index");
2743 if( opt->f_index == V2S_INDEX_VOXEL )
2744 loc_add_2_list(&argv, &acnall, &argc, "voxels");
2745 else
2746 loc_add_2_list(&argv, &acnall, &argc, "nodes");
2747
2748 if( gv2s_plug_opts.gpt_index >= 0 ){
2749 loc_add_2_list(&argv, &acnall, &argc, "-cmask");
2750 if( gv2s_plug_opts.gpt_thresh > 0 )
2751 sprintf(str,"-a %s[%d] -expr astep(a,%f)+equals(a,%f)",
2752 dset_file, gv2s_plug_opts.gpt_index,
2753 gv2s_plug_opts.gpt_thresh,
2754 gv2s_plug_opts.gpt_thresh);
2755 else
2756 sprintf(str,"-a %s[%d] -expr astep(a,%f)",
2757 dset_file, gv2s_plug_opts.gpt_index,
2758 gv2s_plug_opts.gpt_thresh);
2759 loc_add_2_list(&argv, &acnall, &argc, str);
2760 }
2761
2762 if( opt->first_node > 0 ){
2763 loc_add_2_list(&argv, &acnall, &argc, "-first_node");
2764 sprintf(str,"%d",opt->first_node);
2765 loc_add_2_list(&argv, &acnall, &argc, str);
2766 }
2767
2768 if( opt->last_node > 0 && opt->last_node < p->surf[0].num_ixyz-1 ){
2769 loc_add_2_list(&argv, &acnall, &argc, "-last_node");
2770 sprintf(str,"%d",opt->last_node);
2771 loc_add_2_list(&argv, &acnall, &argc, str);
2772 }
2773
2774 if( opt->use_norms ){
2775 loc_add_2_list(&argv, &acnall, &argc, "-use_norms");
2776 if( opt->norm_len != 0.0 ){
2777 loc_add_2_list(&argv, &acnall, &argc, "-norm_len");
2778 loc_add_2_list(&argv, &acnall, &argc,MV_format_fval(opt->norm_len));
2779 }
2780 }
2781
2782 if( opt->f_p1_fr != 0 ){
2783 loc_add_2_list(&argv, &acnall, &argc, "-f_p1_fr");
2784 loc_add_2_list(&argv, &acnall, &argc,MV_format_fval(opt->f_p1_fr));
2785 }
2786
2787 if( opt->f_pn_fr != 0 ){
2788 loc_add_2_list(&argv, &acnall, &argc, "-f_pn_fr");
2789 loc_add_2_list(&argv, &acnall, &argc,MV_format_fval(opt->f_pn_fr));
2790 }
2791
2792 if( opt->f_p1_mm != 0 ){
2793 loc_add_2_list(&argv, &acnall, &argc, "-f_p1_mm");
2794 loc_add_2_list(&argv, &acnall, &argc,MV_format_fval(opt->f_p1_mm));
2795 }
2796
2797 if( opt->f_pn_mm != 0 ){
2798 loc_add_2_list(&argv, &acnall, &argc, "-f_pn_mm");
2799 loc_add_2_list(&argv, &acnall, &argc,MV_format_fval(opt->f_pn_mm));
2800 }
2801
2802 if( opt->oob.show && opt->oob.index > 0 ){
2803 loc_add_2_list(&argv, &acnall, &argc, "-oob_index");
2804 sprintf(str,"%d",opt->oob.index);
2805 loc_add_2_list(&argv, &acnall, &argc, str);
2806 }
2807
2808 if( opt->oob.show && opt->oob.value != 0.0 ){
2809 loc_add_2_list(&argv, &acnall, &argc, "-oob_value");
2810 loc_add_2_list(&argv, &acnall, &argc,MV_format_fval(opt->oob.value));
2811 }
2812
2813 if( DSET_NVALS(p->gpar) > 1 )
2814 loc_add_2_list(&argv, &acnall, &argc, "-outcols_afni_NSD");
2815
2816 if( opt->debug ){
2817 loc_add_2_list(&argv, &acnall, &argc, "-debug");
2818 sprintf(str,"%d",opt->debug);
2819 loc_add_2_list(&argv, &acnall, &argc, str);
2820 }
2821
2822 if( opt->dnode ){
2823 loc_add_2_list(&argv, &acnall, &argc, "-dnode");
2824 sprintf(str,"%d",opt->dnode);
2825 loc_add_2_list(&argv, &acnall, &argc, str);
2826 }
2827
2828 if( opt->segc_file ){
2829 loc_add_2_list(&argv, &acnall, &argc, "-save_seg_coords");
2830 loc_add_2_list(&argv, &acnall, &argc, opt->segc_file);
2831 }
2832
2833 if( opt->outfile_1D ){
2834 loc_add_2_list(&argv, &acnall, &argc, "-out_1D");
2835 loc_add_2_list(&argv, &acnall, &argc, opt->outfile_1D);
2836 }
2837
2838 if( opt->outfile_niml ){
2839 loc_add_2_list(&argv, &acnall, &argc, "-out_niml");
2840 loc_add_2_list(&argv, &acnall, &argc, opt->outfile_niml);
2841 }
2842
2843 /* insert this into the cmd struct */
2844 opt->cmd.fake = 1;
2845 opt->cmd.argc = argc;
2846 opt->cmd.argv = argv;
2847
2848 RETURN(0);
2849 }
2850
2851 /*----------------------------------------------------------------------
2852 * display command
2853 *----------------------------------------------------------------------
2854 */
disp_v2s_command(v2s_opts_t * sopt)2855 int disp_v2s_command( v2s_opts_t * sopt )
2856 {
2857 char * arg;
2858 int ac, quote;
2859
2860 ENTRY("disp_v2s_command");
2861
2862 if( sopt->cmd.argc <= 1 || ! sopt->cmd.argv || ! sopt->cmd.argv[0] )
2863 return 1;
2864
2865 printf("------------------------------------------------------\n"
2866 "+d applying vol2surf similar to the following command:\n");
2867 for( ac = 0; ac < sopt->cmd.argc; ac++ )
2868 if( sopt->cmd.argv[ac] )
2869 {
2870 arg = sopt->cmd.argv[ac];
2871 /* if there are special char, quote the option */
2872 if( strchr(arg, '(') || strchr(arg, '[') ) quote = 1;
2873 else quote = 0;
2874
2875 if( quote ) putchar('\'');
2876 fputs(arg, stdout);
2877 if( quote ) putchar('\'');
2878 putchar(' ');
2879 }
2880 putchar('\n');
2881 fflush(stdout);
2882
2883 RETURN(0);
2884 }
2885
2886
2887 /*----------------------------------------------------------------------
2888 * add 'str' to 'list' via strcpy
2889 *
2890 * if list is not long enough, realloc more pointers
2891 *----------------------------------------------------------------------
2892 */
loc_add_2_list(char *** list,int * nall,int * len,char * str)2893 static int loc_add_2_list( char *** list, int * nall, int * len, char * str)
2894 {
2895 ENTRY("loc_add_2_list");
2896 if( !list || !nall || !len || !str ) RETURN(-1);
2897
2898 if( *nall < 0 ) *nall = 0;
2899 if( *nall == 0 ){ *list = NULL; *nall = 0; *len = 0; } /* init */
2900
2901 if( *nall <= *len ) /* then allocate more memory */
2902 {
2903 *nall += 32;
2904 *list = (char **)realloc(*list, *nall * sizeof(char *));
2905 if(!*list){ fprintf(stderr,"** LA2L: cannot alloc list\n"); RETURN(1); }
2906 }
2907
2908 /* add the string to the list */
2909 (*list)[*len] = strdup(str);
2910 (*len)++;
2911
2912 RETURN(0);
2913 }
2914
2915