1 /*----------------------------------------------------------------------
2 * afni_vol2surf.c - afni's interface to the vol2surf library
3 *
4 * AFNI_vol2surf_func_overlay() - similar to AFNI_vnlist_func_overlay()
5 *----------------------------------------------------------------------
6 */
7
8 /*----------------------------------------------------------------------
9 * R. Reynolds September, 2004
10 *
11 * history:
12 *
13 * 08 Oct 2004 [rickr]:
14 * - AFNI_vol2surf_func_overlay() has new params, surfA,surfB,use_defaults
15 * - pass use_defaults to afni_vol2surf()
16 * 25 Oct 2004 [rickr]:
17 * - accept Rdata and Rthr pointers, for optionally returning the data
18 * and global threshold
19 * - make threshold masks absolute
20 * 09 Aug 2006 [rickr]:
21 * - set the surface volume dataset in the global v2s_plugin_opts struct
22 * - also store the index and threshold value of the threshold sub-brick
23 *----------------------------------------------------------------------
24 */
25
26 #include "afni.h"
27 #include "vol2surf.h"
28
29 static int map_v2s_results(v2s_results *res, Three_D_View *im3d,
30 SUMA_irgba **map, int debug, int dnode );
31 /*-----------------------------------------------------------------------*/
32 /*! Create a nodal color overlay from a voxel map.
33 - Surface index(es) come from global v2s_plugin_opts struct
34 - Return value is number of nodes overlaid:
35 - 0 return ==> no overlay was computed
36 - -1 return ==> some error (e.g., no surface nodes on this dataset)
37 = positive ==> can use *map
38 - *map is set to a newly malloc()-ed array (if return > 0)
39 - if (Rdata) return data from results structure
40 - if (Rthr) return overlay threshold from afni
41
42 Sample usage:
43 - SUMA_irgba * map;
44 - int nmap;
45 - nmap = AFNI_vnlist_func_overlay( im3d, &map );
46 - if( nmap < 0 ){ ** error ** }
47 - else if( nmap == 0 ){ ** nothing to show ** }
48 - else { ** show map[0..nmap-1] ** }
49
50 * based on AFNI_vnlist_func_overlay()
51 -------------------------------------------------------------------------*/
AFNI_vol2surf_func_overlay(Three_D_View * im3d,SUMA_irgba ** map,int surfA,int surfB,int use_defaults,float ** Rdata,float * Rthr)52 int AFNI_vol2surf_func_overlay(Three_D_View *im3d, SUMA_irgba **map,
53 int surfA, int surfB, int use_defaults, float ** Rdata, float * Rthr )
54 {
55 THD_3dim_dataset * oset; /* overlay dataset */
56 THD_slist_find find;
57 THD_session * ss;
58 SUMA_surface * sA, * sB;
59 MRI_IMAGE * im_oim;
60 v2s_plugin_opts * go;
61 v2s_results * results;
62 byte * cmask;
63 int nout, oind, debug;
64
65 ENTRY("AFNI_vol2surf_func_overlay") ;
66
67 /* check inputs */
68 if ( map == NULL || !IM3D_VALID(im3d) ) RETURN(-1);
69
70 go = &gv2s_plug_opts;
71 if ( ! use_defaults && ! go->ready ) RETURN(-1);
72
73 debug = go->sopt.debug; /* because I'm lazy */
74
75 ss = im3d->ss_now; /* session must have needed surface(s) */
76 if( ss == NULL ||
77 ss->su_num <= 0 ||
78 surfA < 0 ||
79 ss->su_num <= surfA ||
80 ss->su_num <= surfB ||
81 ss->su_surf[surfA] == NULL ||
82 (surfB >= 0 && ss->su_surf[surfB] == NULL) )
83 {
84 if ( debug > 1 )
85 {
86 if ( !ss ) fprintf(stderr,"** v2s: NULL session\n");
87 else
88 fprintf(stderr,"** v2s: bad session data:\n"
89 " su_num, surfA, surfB = %d, %d, %d\n",
90 ss->su_num, surfA, surfB);
91 }
92 RETURN(-1);
93 }
94
95 /* init return values */
96 if ( Rdata ) *Rdata = NULL;
97 if ( Rthr ) *Rthr = 0.0;
98
99 /* set surface pointers */
100 sA = ss->su_surf[surfA];
101 sB = ( surfB >= 0 ) ? ss->su_surf[surfB] : NULL;
102
103 /* store the surface volume pointer 9 Aug 2006 [rickr] */
104 find = PLUTO_dset_finder(sA->idcode_dset);
105 go->sv_dset = find.dset;
106
107 if ( debug )
108 {
109 fprintf(stderr,"++ v2s overlay: sa,sb = %d,%d\n", surfA, surfB);
110 if ( debug > 1 )
111 fprintf(stderr," surfA is %s, surfB is %s\n",
112 sA->label[0] ? sA->label : "<no label>",
113 sB ? (sB->label[0] ? sB->label : "<no label>") : "<not used>");
114 }
115
116 /*-------------------- overlay image --------------------*/
117 oset = im3d->fim_now;
118 oind = im3d->vinfo->fim_index; /* overlay sub-brick index */
119
120 if( oset == NULL )
121 {
122 if ( debug > 1 ) fprintf(stderr,"** v2s: no overlay dset\n");
123 RETURN(-1);
124 }
125 im_oim = DSET_BRICK(oset,oind);
126
127 if( ! im_oim || ! AFNI_GOOD_FUNC_DTYPE(im_oim->kind) )
128 {
129 if ( debug > 1 ) fprintf(stderr,"** v2s: no overlay image\n");
130 RETURN(-1);
131 }
132 DSET_load(oset); /* to be sure */
133 if( !DSET_LOADED(oset) ) RETURN(-1);
134
135 /*-------------------- mask from threshold --------------------*/
136 go->gpt_index = -1; go->gpt_thresh = 0.0; /* init threshold options */
137 cmask = NULL;
138
139 /* if a new overlay has been created from clustering, get it as a mask
140 (for now - since it's an MRI_IMAGE, and we can't pass to v2s)
141 11 Jan 2008 [rickr] */
142
143 if( debug ) {
144 fprintf(stderr,"-- CLUST: ");
145 if( ! im3d->vinfo )
146 fprintf(stderr,"no vinfo to cluster on\n");
147 else if( ! oset->dblk->vedim )
148 fprintf(stderr,"no dblk->vedim to cluster on\n");
149 else {
150 fprintf(stderr,"fim_index %d, thr_index %d, VEDIT_IVAL %d\n",
151 im3d->vinfo->fim_index,
152 im3d->vinfo->thr_index,
153 DSET_VEDIT_IVAL(oset));
154 fprintf(stderr," (kind %d, brick type %d)\n",
155 oset->dblk->vedim->kind, DSET_BRICK_TYPE(oset,oind));
156 }
157 }
158
159 if( im3d->vinfo && oset->dblk->vedim &&
160 (im3d->vinfo->fim_index == DSET_VEDIT_IVAL(oset)) &&
161 (oset->dblk->vedim->kind == DSET_BRICK_TYPE(oset,oind)) ) {
162 cmask = mri_to_bytemask(oset->dblk->vedim, 1.0, 0.0);
163 if( debug > 1 ) fprintf(stderr,"++ applying mask from edited image\n");
164 }
165 else if( im3d->vinfo->func_threshold > 0.0 &&
166 im3d->vinfo->thr_onoff ) /* then want a threshold */
167 {
168 MRI_IMAGE * im_thr;
169 float thresh;
170 int tind = im3d->vinfo->thr_index;
171
172 im_thr = DSET_BRICK(oset,tind);
173
174 /* note real threshold */
175 thresh = get_3Dview_func_thresh(im3d,1) ;
176
177 /* maybe we want to return this */
178 if ( Rthr ) *Rthr = thresh;
179
180 if( im_thr && !AFNI_GOOD_FUNC_DTYPE(im_thr->kind) )
181 im_thr = NULL;
182
183 /* create the mask, if this fails, just continue... */
184 if( im_thr != NULL )
185 {
186 int nset;
187 if ( debug > 1 )
188 fprintf(stderr,"++ mask from index %d and thresh %f\n",
189 tind,thresh);
190 nset = thd_mask_from_brick(oset, tind, thresh, &cmask, 1);
191 if ( debug > 1 )
192 {
193 if ( ! cmask )
194 fprintf(stderr,"-- no mask created\n");
195 else
196 fprintf(stderr,"++ mask has %d set voxels\n", nset);
197 }
198 }
199 else if ( debug > 1 )
200 fprintf(stderr,"-- no threshold mask\n");
201
202 /* note mask options for v2s command output 9 Aug 2006 [rickr] */
203 if( cmask ) { go->gpt_index = tind; go->gpt_thresh = thresh; }
204 } else { /* then make the mask from non-zero overlay values */
205 cmask = THD_makemask(oset, oind, 1.0, 0.0);
206 /* we are applying a cmask, so make it clear in any printed command
207 8 Dec 2016 [rickr] */
208 go->gpt_index = oind; go->gpt_thresh = 0;
209 }
210
211 /*-------------------- vol2surf computation --------------------*/
212 results = afni_vol2surf(oset, oind, sA, sB, cmask, use_defaults);
213
214 if ( cmask ) free(cmask); /* we're done with the mask */
215 if ( ! results )
216 {
217 if ( debug ) fprintf(stderr,"-- vol2surf failure\n");
218 RETURN(-1); /* failure? */
219 }
220
221 /*-------------------- lose old vnlist --------------------*/
222 /* vol2surf has no method for counting voxels (right now) */
223 if( sA->vn )
224 {
225 if ( debug > 1 ) fprintf(stderr,"-- destroying vnlist\n");
226 SUMA_destroy_vnlist( sA->vn ) ;
227 sA->vn = NULL;
228 }
229
230 /*-------------------- set overlay colors --------------------*/
231 nout = map_v2s_results(results, im3d, map, debug, go->sopt.dnode);
232
233 /* before free()ing results, check whether we want to return the values */
234 if ( Rdata )
235 {
236 *Rdata = results->vals[0];
237 results->vals[0] = NULL; /* do not let free_v2s_results() free it */
238 }
239
240 free_v2s_results(results);
241
242 if ( debug > 1 ) fprintf(stderr,"++ map_v2s_results: nout = %d\n", nout);
243
244 /* check failure */
245 if ( ! *map ) RETURN(-1);
246
247 RETURN(nout);
248 }
249
250 /*! given vol2surf results, fill SUMA_irgba struct */
map_v2s_results(v2s_results * res,Three_D_View * im3d,SUMA_irgba ** map,int debug,int dnode)251 static int map_v2s_results(v2s_results *res, Three_D_View *im3d,
252 SUMA_irgba **map, int debug, int dnode )
253 {
254 SUMA_irgba * mptr;
255 MCW_pbar * pbar;
256 rgbyte * cmap;
257 float * result_vals;
258 float bbot, btop, fac, fval;
259 float pane_scale = 1.0;
260 byte r, g, b;
261 int nindex, node, npanes, ival;
262 int map_index = 0; /* counter for results */
263
264 ENTRY("map_v2s_results");
265
266 /*-------------------- create output node list --------------------*/
267 /* from malloc 12 Feb 2009 [lesstif patrol] */
268 mptr = (SUMA_irgba *) calloc(res->nused, sizeof(SUMA_irgba));
269 if ( ! mptr )
270 {
271 fprintf(stderr,"** malloc failure: %d irgba's\n", res->nused);
272 RETURN(-1);
273 }
274 *map = mptr; /* and store the address */
275
276 pbar = im3d->vwid->func->inten_pbar;
277 npanes = pbar->num_panes;
278
279 /* with AFNI_PBAR_FULLRANGE, pane_scale = 1.0 [20 May 2019 rickr] */
280 /*
281 pane_scale = im3d->vinfo->fim_range;
282 if ( pane_scale == 0.0 ) pane_scale = im3d->vinfo->fim_autorange;
283 if ( pane_scale == 0.0 ) pane_scale = 1.0;
284 */
285 if( PBAR_FULLRANGE ) {
286 pane_scale = 1.0;
287 if( debug > 2 )
288 fprintf(stderr,"+d mvr: have PBAR_FULLRANGE, resetting pane_scale\n");
289 } else {
290 pane_scale = FIM_RANGE(im3d);
291 if( debug > 2 )
292 fprintf(stderr,"+d mvr: no PBAR_FULLRANGE, pane_scale = FIM_RANGE\n");
293 }
294
295 if ( debug > 1 ) {
296 fprintf(stderr,"+d mvr: npanes = %d, pane_scale = %f\n",
297 pbar->bigmode ? NPANE_BIG : npanes, pane_scale);
298 fprintf(stderr," fim_range = %f, fim_autorange = %f\n",
299 im3d->vinfo->fim_range, im3d->vinfo->fim_autorange);
300 }
301
302 result_vals = res->vals[0]; /* for typing and potential speed */
303
304 if( pbar->bigmode ) /* colorscale */
305 {
306 int zbot;
307
308 bbot = pane_scale*pbar->bigbot;
309 btop = pane_scale*pbar->bigtop;
310 fac = NPANE_BIG / (btop-bbot);
311 cmap = pbar->bigcolor;
312 zbot = (bbot == 0.0);
313
314 if ( debug > 1 ) {
315 fprintf(stderr,"+d bigmode: bbot,btop,fac, zbot = %f,%f,%f, %d\n",
316 bbot, btop, fac, zbot);
317 fprintf(stderr," NPANE_BIG = %d, BIGGEST = %d\n",
318 NPANE_BIG, NPANE_BIGGEST);
319 }
320
321 for ( nindex = 0; nindex < res->nused; nindex++ )
322 {
323 if ( res->nodes ) node = res->nodes[nindex];
324 else node = nindex;
325
326 /* note value, and check to ignore */
327 fval = result_vals[nindex];
328
329 if ( debug > 1 && node == dnode )
330 fprintf(stderr, "+d dnode %d, fval %f\n", dnode, fval);
331
332 if ( zbot && fval < 0 ) continue;
333
334 /* note the color panel index, and bound it in [0,NPANE_BIG-1] */
335 if( fval >= btop ) ival = 0; /* guard against overflow */
336 else if ( fval <= bbot ) ival = NPANE_BIG - 1;
337 /* trunc: matching afni_func.c [20 May 2019 rickr] */
338 else ival = (int)(fac * (btop - fval));
339
340 if ( ival < 0 ) ival = 0;
341 if ( ival >= NPANE_BIG ) ival = NPANE_BIG - 1;
342
343 /* get color map, and check if non-zero */
344 r = cmap[ival].r; g = cmap[ival].g; b = cmap[ival].b;
345
346 if ( debug > 1 && node == dnode )
347 fprintf(stderr, "+d pane, r,g,b : %d, %d,%d,%d\n",ival,r,g,b);
348
349 if ( r == 0 && g == 0 && b == 0 ) continue;
350
351 /* we are set, fill the SUMA_irgba struct */
352 mptr->id = res->nodes[nindex];
353 mptr->r = r; mptr->g = g; mptr->b = b; mptr->a = 255;
354
355 /* increment pointer and counter (okay, so counter is unneeded) */
356 mptr++; map_index++;
357 }
358 }
359 else /* indexed colors */
360 {
361 float othr[NPANE_MAX+1]; /* threshold */
362 short ovc[NPANE_MAX+1]; /* color */
363 byte ovc_r[NPANE_MAX+1], ovc_g[NPANE_MAX+1], ovc_b[NPANE_MAX+1];
364
365 /* set the overlay color indices */
366 for( ival=0 ; ival < npanes ; ival++ )
367 ovc[ival] = pbar->ov_index[ival]; /* from top of pbar down */
368 ovc[npanes] = im3d->vinfo->use_posfunc ? 0 : ovc[npanes-1];
369
370 /* get the actual RGB colors of each pane on the pbar */
371 for( ival=0 ; ival <= npanes ; ival++ ) /* include npanes 27 Aug 2007 */
372 {
373 ovc_r[ival] = DCOV_REDBYTE (im3d->dc,ovc[ival]);
374 ovc_g[ival] = DCOV_GREENBYTE(im3d->dc,ovc[ival]);
375 ovc_b[ival] = DCOV_BLUEBYTE (im3d->dc,ovc[ival]);
376 }
377
378 /* compute the thresholds */
379 for( ival=0 ; ival < npanes ; ival++ )
380 othr[ival] = pane_scale * pbar->pval[ival+1];
381 othr[npanes] = othr[npanes-1] ;
382
383 if ( debug > 2 )
384 {
385 for( ival=0 ; ival <= npanes ; ival++ )
386 fprintf(stderr,"+d pane #%2d, ovc = %d, othr = %f\n",
387 ival, ovc[ival], othr[ival]);
388 }
389
390 for ( nindex = 0; nindex < res->nused; nindex++ )
391 {
392 if ( res->nodes ) node = res->nodes[nindex];
393 else node = nindex;
394
395 /* note value, and check to ignore */
396 fval = result_vals[nindex];
397
398 if ( debug > 1 && node == dnode )
399 fprintf(stderr, "+d dnode %d, fval %f\n", dnode, fval);
400
401 if ( fval == 0.0 ) continue;
402
403 for ( ival = 0; ival < npanes && fval < othr[ival]; ival++ )
404 ;
405 if ( ovc[ival] == 0 ) continue; /* no color in this pane */
406 r = ovc_r[ival]; g = ovc_g[ival]; b = ovc_b[ival];
407
408 /* we are set, fill the SUMA_irgba struct */
409 mptr->id = res->nodes[nindex];
410 mptr->r = r; mptr->g = g; mptr->b = b; mptr->a = 255;
411
412 if ( debug > 1 && node == dnode )
413 fprintf(stderr, "+d pane, r,g,b = %d, %d,%d,%d\n",ival,r,g,b);
414
415 if ( r == 0 && g == 0 && b == 0 ) continue;
416
417 /* increment pointer and counter (okay, so counter is unneeded) */
418 mptr++; map_index++;
419 }
420 }
421
422 if ( debug > 0 )
423 fprintf(stderr,"+d mvr v2s map size = %d\n", map_index);
424
425 /* forfeit unused memory */
426 *map = (SUMA_irgba *)realloc(*map, sizeof(SUMA_irgba)*map_index);
427 if ( ! *map )
428 {
429 fprintf(stderr,"** mvr: failed realloc of map, %d irgba's\n",map_index);
430 map_index = -1;
431 }
432
433 RETURN(map_index) ; /* number of entries in map */
434 }
435
436
437 /*-----------------------------------------------------------------------*/
438 /*! Return the vol2surf time series data at a single node.
439 *
440 * sess : session to get surfaces from (e.g. im3d->ss_now)
441 * dset : dataset to get values from (e.g. im3d->anat_now)
442 * surf : surface index (will use LDP of this surface)
443 * (e.g. in {0..from im3d->ss_now->su_num-1})
444 * node : node index in surface
445 *
446 * Note: surf/node might come from AFNI_get_xhair_node(), say.
447 *
448 * returns: a pointer to the float data or NULL on failure.
449 *
450 * Array will be of length DSET_NVALS(dset), and should be free'd by
451 * the calling function.
452 * ----------------------------------------------------------------------*/
AFNI_v2s_node_timeseries(THD_session * sess,THD_3dim_dataset * dset,int surfA,int surfB,int node,int use_v2s)453 float * AFNI_v2s_node_timeseries(THD_session * sess, THD_3dim_dataset * dset,
454 int surfA, int surfB, int node, int use_v2s)
455 {
456 v2s_results * results = NULL;
457 v2s_opts_t sopt;
458 SUMA_surface * sA, * sB;
459 float * ts = NULL; /* time series to return */
460 int ind, len, verb = gv2s_plug_opts.sopt.debug;
461
462 ENTRY("AFNI_v2s_node_timeseries");
463
464 if( !sess || !dset || surfA < 0 ){
465 fprintf(stderr,"** v2s_node_ts - bad inputs: %p,%p,%d\n",sess,dset,surfA);
466 RETURN(NULL);
467 }
468
469 sA = (surfA >= 0) ? sess->su_surf[surfA] : NULL;
470 sB = (surfB >= 0) ? sess->su_surf[surfB] : NULL;
471
472 /* check and apply the node index */
473 if( node < 0 || node >= sA->num_ixyz ) {
474 fprintf(stderr,"** v2s_NTS: node %d outside [0,%d)\n",node,sA->num_ixyz);
475 RETURN(NULL);
476 } else if( verb > 1 )
477 fprintf(stderr,"-- v2s_NTS: getting time series at node %d\n", node);
478
479 /* fill options struct (from defaults or plugin opts) */
480 if( use_v2s ) {
481 /* is using plugin, nuke any gp_index, filenames or structs */
482 sopt = gv2s_plug_opts.sopt;
483
484 sopt.gp_index = -1; /* get all sub-bricks */
485 sopt.outfile_1D = sopt.outfile_niml = sopt.segc_file = NULL;
486 memset(&sopt.cmd, 0, sizeof(v2s_cmd_t));
487 memset(&sopt.oob, 0, sizeof(v2s_oob_t));
488 memset(&sopt.oom, 0, sizeof(v2s_oob_t));
489 } else
490 v2s_fill_sopt_default(&sopt, sB ? 2 : 1);
491
492 sopt.skip_cols = V2S_SKIP_ALL ^ V2S_SKIP_VALS; /* get all data */
493 sopt.first_node = sopt.last_node = node; /* get only 1 node */
494 if( verb > 2 ) sopt.dnode = node; /* maybe make verbose */
495
496 /* get the results */
497 results = opt_vol2surf(dset, &sopt, sA, sB, NULL);
498 if( !results ) RETURN(NULL);
499
500 /* extract single time series */
501 len = DSET_NVALS(dset);
502 if( len != results->nlab ){
503 fprintf(stderr,"** v2s_NTS: nvals != nlab (%d, %d)\n",len,results->nlab);
504 free_v2s_results(results); RETURN(NULL);
505 }
506
507 ts = (float *)malloc(len*sizeof(float));
508 if( !ts ) {
509 fprintf(stderr,"** v2s_NTS: failed to alloc %d floats\n", len);
510 free_v2s_results(results); RETURN(NULL);
511 }
512
513 /* actually fill the array, woohoo! */
514 if(verb > 3) fprintf(stderr,"-- v2s_NTS: copying %d vals for return\n",len);
515 for( ind = 0; ind < len; ind++ )
516 ts[ind] = results->vals[ind][0];
517
518 if(verb > 3) fprintf(stderr,"-- v2s_NTS: freeing results... ");
519 free_v2s_results(results);
520
521 if(verb > 3) fprintf(stderr,"done\n");
522 RETURN(ts);
523 }
524
525