1 #include "SUMA_suma.h"
2
3 extern SUMA_CommonFields *SUMAg_CF;
4 extern SUMA_DO *SUMAg_DOv;
5 extern SUMA_SurfaceViewer *SUMAg_SVv;
6 extern int SUMAg_N_SVv;
7 extern int SUMAg_N_DOv;
8
SUMA_DotXform_MakeOrts(NI_element * dotopt,int ts_len,int polort,char * ortname)9 SUMA_Boolean SUMA_DotXform_MakeOrts( NI_element *dotopt, int ts_len,
10 int polort, char *ortname)
11 {
12 static char FuncName[]={"SUMA_DotXform_MakeOrts"};
13 char stmp[256];
14 float *fort=NULL, **refvec=NULL;
15 int nort=0, nts = 0 ;
16 int suc=0, i=0, nref=0;
17 SUMA_Boolean LocalHead = NOPE;
18
19 SUMA_ENTRY;
20
21
22 if (ortname) {
23 SUMA_LHv("Loading %s\n", ortname);
24 fort = SUMA_Load1D_s(ortname, &nort, &nts, 0, 0);
25 if (!fort) {
26 SUMA_S_Err("Could not load orts");
27 SUMA_RETURN(NOPE);
28 }
29 if (nts!= ts_len) {
30 SUMA_S_Err("mismatch between polort length and time series length");
31 SUMA_RETURN(NOPE);
32 }
33 } else {
34 fort = NULL;
35 nort = 0;
36 }
37
38 /* form the baseline regressors */
39 nref = polort+1;
40 if (nref >= ts_len-3) {
41 SUMA_S_Errv("Number of baseline regressors(%d) \n"
42 "is too high compared to the number of samples (%d)\n"
43 , nref, ts_len);
44 SUMA_RETURN(NOPE);
45 }
46 if (nref) {
47 SUMA_LHv("Building %d polyrefs\n", nref);
48 refvec = THD_build_polyref( nref ,
49 ts_len ) ;
50 } else {
51 SUMA_LH("No polyrefs");
52 refvec = NULL;
53 nref = 0;
54 }
55
56 /* now add the orts */
57 if (fort) {
58 SUMA_LHv("Adding %d from fort, length %d to previous %d orts\n",
59 nort, nts, nref);
60 refvec = (float **)SUMA_realloc(refvec, (nref+nort)*sizeof(float*));
61 for (i=0; i<nort; ++i) {
62 refvec[i+nref] = (float *)SUMA_calloc(nts, sizeof(float));
63 memcpy(refvec[i+nref], &(fort[i*nts]), sizeof(float)*nts);
64 }
65 free(fort); fort = NULL;
66 nref += nort;
67 } else {
68 SUMA_LH("No fort\n");
69 }
70
71
72 /* put a copy of the regressors in dotopt */
73 if (dotopt->vec_num) {
74 SUMA_LH("Cleaning up dotopt");
75 while (dotopt->vec_num) {
76 NI_remove_column(dotopt,-1);
77 }
78 }
79
80 SUMA_LHv("Adding %d columns of length %d\n", nref, ts_len);
81 /* add the columns here */
82 NI_alter_veclen(dotopt, ts_len);
83 for (i=0; i<nref; ++i) {
84 NI_add_column(dotopt, NI_FLOAT, refvec[i]);
85 }
86
87 if (LocalHead) {
88 sprintf(stmp,"file:%s.dotopt.1D", FuncName);
89 SUMA_LHv("Writing %s\n", stmp);
90 NEL_WRITE_1D(dotopt, stmp, suc);
91 }
92
93 /* clean regressors */
94 if (refvec) {
95 for (i=0; i<nref; ++i) {
96 if (refvec[i]) free(refvec[i]); refvec[i] = NULL;
97 }
98 free(refvec); refvec=NULL;
99 }
100
101 /* flag num_ort_parameters as not set */
102 NI_SET_INT(dotopt, "num_ort_parameters", -1);
103
104 SUMA_RETURN(YUP);
105 }
106
107 /*
108 Function to create dotopt NI_element, set some of its
109 parameters, AND recreate the vector of regressors.
110
111 doptop (NI_element *) if NULL, a new one is created
112 ts_len (int): length of regressor time series
113 ftop: if > 0 set the top pass frequency
114 else nothing is done
115 fbot: if > 0 set the bot pass frequency
116 else nothing is done
117 norm: if == 1 set normalize_dset = "y"
118 if == 0 set normalize_dset = "n"
119 else do nothing
120 prec: if > 0 set precision
121 else do nothing
122 polort: if > -2 set polort value and create regressors
123 else whatever is in doptopt is used.
124 ortname: if not NULL, load orts from file ortname and add
125 the other regressors
126 WARNING!: If polort, or ortname need changing, then both
127 parameters must be set.
128 */
SUMA_set_dotopts(NI_element * dotopt,int ts_len,float ftop,float fbot,int norm,int prec,int polort,char * ortname)129 NI_element *SUMA_set_dotopts(NI_element *dotopt, int ts_len,
130 float ftop, float fbot,
131 int norm, int prec,
132 int polort, char *ortname)
133 {
134 static char FuncName[]={"SUMA_set_dotopts"};
135 SUMA_Boolean LocalHead = NOPE;
136
137 SUMA_ENTRY;
138
139 if (ts_len < 2) {
140 SUMA_S_Errv("bad ts_len of %d\n", ts_len);
141 SUMA_RETURN(dotopt);
142 }
143
144 if (!dotopt) { /* new with defaults */
145 dotopt = NI_new_data_element("dotopts", 0);
146 NI_SET_FLOAT (dotopt,"filter_above", 99999999.9);
147 NI_SET_FLOAT (dotopt,"filter_below", 0.0);
148 NI_SET_INT(dotopt,"polort",-1);
149 NI_SET_INT(dotopt,"prec",1);
150 }
151
152 /* initialize by input */
153 if (ftop > 0) NI_SET_FLOAT(dotopt, "filter_above", ftop);
154 if (fbot >= 0) NI_SET_FLOAT(dotopt, "filter_below", fbot);
155 if (norm == 1) {
156 NI_set_attribute(dotopt, "normalize_dset", "y");
157 } else if (norm == 0) {
158 NI_set_attribute(dotopt, "normalize_dset", "n");
159 }
160 if (prec > 0) NI_SET_INT(dotopt,"numeric_precision", prec);
161 if (polort > -2) NI_SET_INT(dotopt,"polort", polort);
162 if (ortname) NI_set_attribute(dotopt,"ortname",ortname);
163
164 /* now enforce settings */
165
166 /* get the file-based ort functions */
167 ortname = NI_get_attribute(dotopt,"ortname");
168 NI_GET_INT(dotopt,"polort", polort);
169 if (!NI_GOT) polort = -1;
170
171 if (!SUMA_DotXform_MakeOrts( dotopt, ts_len, polort, ortname)){
172 SUMA_S_Err("Failed to make orts");
173 SUMA_RETURN(dotopt);
174 }
175
176 /* get the filtering option */
177 NI_GET_FLOAT(dotopt, "filter_above", ftop);
178 if (!NI_GOT) ftop = 99999999.9; /* Hz, what elese did you expect? */
179 NI_GET_FLOAT(dotopt, "filter_below", fbot);
180 if (!NI_GOT) fbot = 0;
181
182 /* initialize pending to nothing */
183 NI_set_attribute(dotopt, "pending", "");
184
185 if (LocalHead) {
186 SUMA_ShowNel(dotopt);
187 }
188 SUMA_RETURN(dotopt);
189 }
190
SUMA_DotXform_GetRecomputeForDset(NI_element * dotopts,char * id)191 int SUMA_DotXform_GetRecomputeForDset (NI_element *dotopts, char *id)
192 {
193 static char FuncName[]={"SUMA_DotXform_GetRecomputeForDset"};
194 int recompute = 0;
195 char *cs;
196
197 SUMA_ENTRY;
198
199 cs = NI_get_attribute(dotopts, "pending");
200 if (strstr(cs, id)) SUMA_RETURN(1);
201 else SUMA_RETURN(0);
202
203 }
204
SUMA_DotXform_SetPending(NI_element * dotopts,int pending,char * id)205 void SUMA_DotXform_SetPending (NI_element *dotopts, int pending, char *id)
206 {
207 static char FuncName[]={"SUMA_DotXform_SetPending"};
208 int ii;
209 char stmp[10*SUMA_IDCODE_LENGTH+11]={""};
210 char *cs=NULL;
211 SUMA_Boolean LocalHead = NOPE;
212
213 SUMA_ENTRY;
214
215
216 if (!dotopts) {
217 SUMA_S_Err("No dotopts");
218 SUMA_RETURNe;
219 }
220 if (!id && pending) {
221 SUMA_S_Err("Cannot set pending to 1 with no id");
222 SUMA_RETURNe;
223 }
224 if (!pending) {
225 if (!id) {
226 /* kill all */
227 NI_set_attribute(dotopts, "pending", "");
228 } else {
229 cs = NI_get_attribute(dotopts, "pending");
230 SUMA_LHv("Have pending of \n"
231 ">>>%s<<<\n", cs);
232 if (SUMA_Remove_Sub_String(cs, ";", id) == 1) {
233 NI_set_attribute(dotopts, "pending", cs);
234 }
235 }
236 } else {
237 {
238 cs = NI_get_attribute(dotopts,"pending");
239 if (cs) {
240 if (!strstr(cs, id)) {
241 /* add it */
242 strcat(stmp, cs);
243 strcat(stmp, id);
244 NI_set_attribute(dotopts, "pending", stmp);
245 } else {/* alread pending */
246 }
247 } else { /* nothing pending, add it */
248 sprintf(stmp, "%s;", id);
249 NI_set_attribute(dotopts, "pending", stmp);
250 }
251 }
252 }
253
254 SUMA_RETURNe;
255 }
256
SUMA_GetDotPreprocessedDset(SUMA_DSET * in_dset,NI_element * dotopt)257 SUMA_DSET *SUMA_GetDotPreprocessedDset(SUMA_DSET *in_dset,
258 NI_element *dotopt)
259 {
260 static char FuncName[]={"SUMA_GetDotPreprocessedDset"};
261 SUMA_DSET *pdset=NULL;
262 float **refvec=NULL, ftop=9999999.9, fbot = 0.0;
263 int recompute=0, i, N_ort_param=0;
264 char stmp[256], *ppid=NULL;
265 SUMA_Boolean LocalHead = NOPE;
266
267 SUMA_ENTRY;
268
269 /* look for pre-existing dset */
270 SUMA_LH("Looking for preprocessed dset");
271
272 /* first recreate the unique identifier of the preprocessed dset */
273 sprintf(stmp,"dot.preprocess.%s", SDSET_ID(in_dset));
274 ppid = UNIQ_hashcode(stmp);
275
276 recompute = SUMA_DotXform_GetRecomputeForDset(dotopt, SDSET_ID(in_dset));
277
278 /* search if that dset exists already */
279 if (!recompute &&
280 (pdset = SUMA_FindDset_s(ppid, SUMAg_CF->DsetList))) {
281 SUMA_LH("got it!");
282 /* got it, get out */
283 SUMA_free(ppid); ppid = NULL;
284 SUMA_RETURN(pdset);
285 } else {/* (re) create the bastard */
286 SUMA_LH("Don't got it!, or recompute");
287 /* make sure dotopt is ready (detrending vectors are loaded)*/
288 if (!dotopt->vec_num) {
289 dotopt = SUMA_set_dotopts(dotopt, SDSET_VECNUM(in_dset),
290 -1, -1,
291 -1, -1,
292 -2, NULL);
293 }
294
295 if (!dotopt->vec_num) {
296 SUMA_S_Err("Nothing to do here!");
297 goto BAIL;
298 }
299
300 /* build a refvec */
301 refvec = (float **)SUMA_calloc(dotopt->vec_num, sizeof(float *));
302 for (i=0; i<dotopt->vec_num; ++i) {
303 refvec[i] = (float *)dotopt->vec[i];
304 }
305 /* now detrend the whole thing */
306 NI_GET_FLOAT(dotopt,"filter_below",fbot);
307 NI_GET_FLOAT(dotopt,"filter_above",ftop);
308
309 SUMA_LHv("Detrending with polort of 2, %d extra orts and BP %f %f\n",
310 dotopt->vec_num, fbot, ftop);
311 pdset = SUMA_DotDetrendDset(in_dset, refvec, dotopt->vec_num,
312 fbot, ftop, 1, &N_ort_param);
313 /* the '1' flag forces a second order detrending,
314 in addition to orts and smoothing. This is done
315 in keeping with AFNI's implementation */
316 NI_SET_INT(dotopt,"num_ort_parameters", N_ort_param);
317 NI_set_attribute(pdset->ngr, "self_idcode", ppid);
318 NI_set_attribute(pdset->ngr,"domain_parent_idcode",
319 SDSET_IDMDOM(in_dset));
320 NI_set_attribute(pdset->ngr,"geometry_parent_idcode",
321 SDSET_IDGDOM(in_dset));
322 SUMA_free(ppid); ppid = NULL;
323
324 /* put the dset in the global list (allow for replace)*/
325 if (!SUMA_InsertDsetPointer(&pdset, SUMAg_CF->DsetList,1)){
326 SUMA_S_Err("Failed to insert pointer");
327 goto BAIL;
328 }
329
330
331 #if 1 /* make it displayable */
332 {
333 SUMA_OVERLAYS *child=NULL;
334 SUMA_SurfaceObject *SO=NULL;
335 char *stmp=NULL;
336 SUMA_LIST_WIDGET *LW = NULL;
337 int OverInd=0;
338
339 SUMA_LH("Create its overlay (child)");
340
341 if (!(SO =
342 SUMA_findSOp_inDOv(SDSET_IDMDOM(pdset), SUMAg_DOv, SUMAg_N_DOv))) {
343 SUMA_S_Errv("Failed to find domain surface '%s'\n",
344 SDSET_IDMDOM(pdset));
345 if (LocalHead) {
346 SUMA_ShowDset(pdset, 0, NULL);
347 SUMA_ShowDset(in_dset, 0, NULL);
348 }
349 goto BAIL;
350 }
351
352 if (!(child = SUMA_Fetch_OverlayPointerByDset (
353 (SUMA_ALL_DO *)SO, pdset, &OverInd))) {
354 /* need a new one */
355 if (!(SDSET_LABEL(pdset))) SUMA_LabelDset(pdset,NULL);
356 stmp = SUMA_append_string("dotprep.", SDSET_LABEL(pdset));
357 child = SUMA_CreateOverlayPointer (
358 stmp, pdset, SO->idcode_str, NULL);
359 SUMA_free(stmp); stmp=NULL;
360 if (!child) {
361 SUMA_S_Err("Failed in CreateOverlayPointer." );
362 goto BAIL;
363 }
364 SUMA_LH("Add overlay to SO");
365 /* Add this plane to SO->Overlays */
366 if (!SUMA_AddNewPlane ( (SUMA_ALL_DO *)SO, child,
367 SUMAg_DOv, SUMAg_N_DOv, 0)) {
368 SUMA_SL_Crit("Failed in SUMA_AddNewPlane");
369 SUMA_FreeOverlayPointer(child);
370 goto BAIL;
371 }
372
373 /* set the opacity, index column and the range */
374 child->GlobalOpacity = YUP;
375 child->ShowMode = SW_SurfCont_DsetViewCol;
376 child->OptScl->BrightFact = 0.6;
377
378 child->OptScl->find = 0;
379 child->OptScl->tind = 0;
380 child->OptScl->bind = 0;
381
382 SUMA_LH("Refreshing Dset list");
383 /*update the list widget if open */
384 LW = SO->SurfCont->SwitchDsetlst;
385 if (LW) {
386 if (!LW->isShaded) SUMA_RefreshDsetList ((SUMA_ALL_DO *)SO);
387 }
388
389 /* this chunk from SUMA_ColPlane_NewOrder */
390 if (!SUMA_MovePlaneDown((SUMA_ALL_DO *)SO, child->Name)) {
391 SUMA_L_Err("Error in SUMA_MovePlaneUp.");
392 goto BAIL;
393 }
394 }
395 }
396 #endif
397
398 SUMA_DotXform_SetPending (dotopt, 0, SDSET_ID(in_dset));
399
400 SUMA_RETURN(pdset);
401 }
402
403 BAIL:
404 if (ppid) SUMA_free(ppid); ppid = NULL;
405 if (refvec) {
406 free(refvec); refvec=NULL;
407 }
408 SUMA_RETURN(NULL);
409 }
410
SUMA_DotPreProcessTimeSeries(float * fv,int N_ts,float TR,NI_element * dotopts)411 double *SUMA_DotPreProcessTimeSeries(float *fv, int N_ts,
412 float TR, NI_element *dotopts)
413 {
414 static char FuncName[]={"SUMA_DotPreProcessTimeSeries"};
415 float **ort=NULL, ftop=9999999.9, fbot=0.0;
416 double *ts=NULL;
417 int ii =0;
418 SUMA_Boolean LocalHead = NOPE;
419
420 SUMA_ENTRY;
421
422 if (!fv || !dotopts) SUMA_RETURN(NULL);
423
424 /* detrend it ? */
425 if (NI_YES_ATTR(dotopts,"normalize_dset")) {
426 if (dotopts->vec_num) { /* something to ort */
427 if (dotopts->vec_len != N_ts) {
428 SUMA_S_Err("bad dotopts->vec_len");
429 SUMA_RETURN(NULL);
430 }
431 ort = (float **)SUMA_calloc(dotopts->vec_num, sizeof(float *));
432 for (ii=0; ii<dotopts->vec_num; ++ii)
433 ort[ii] = (float *)dotopts->vec[ii];
434 } else {
435 ort = NULL;
436 }
437 NI_GET_FLOAT(dotopts, "filter_above", ftop);
438 if (!NI_GOT) ftop = 99999999.9; /* Hz, what elese did you expect? */
439 NI_GET_FLOAT(dotopts, "filter_below", fbot);
440 if (!NI_GOT) fbot = 0.0;
441
442 SUMA_LHv("HERE: fbot: %f\tftop: %f\n", fbot, ftop);
443 /* SUMA_ShowNel(dotopts); */
444
445 if (!THD_bandpass_vectors( N_ts , 1 , &fv ,
446 (float)TR , fbot , ftop ,
447 0 , dotopts->vec_num , ort )) {
448 SUMA_S_Err("Bad bandpass call");
449 SUMA_RETURN(NULL);
450 }
451 SUMA_LH("Back from purgatorium");
452 if (ort) free(ort); ort = NULL;
453 }
454
455 THD_normalize( N_ts , fv ) ;
456
457
458 /* now stash fv in ts */
459 ts = (double *)SUMA_calloc(N_ts, sizeof(double));
460 for (ii=0; ii<N_ts; ++ii) {
461 ts[ii] = (double)fv[ii];
462 }
463 SUMA_RETURN(ts);
464 }
465
SUMA_dot_product_CB(void * params)466 void SUMA_dot_product_CB( void *params)
467 {
468 static char FuncName[]={"SUMA_dot_product_CB"};
469
470 char *SO_idcode=NULL, *ts_dset_id=NULL,
471 *dot_dset_id=NULL, stmp[300], ident[300], prefix[300],
472 *p1=NULL, *p2=NULL, *Cside;
473 SUMA_DSET *in_dset=NULL, *ts_src_dset=NULL;
474 double TR = 0;
475 double *ts=NULL;
476 int ts_node=0, N_ts, ii, jj;
477 float *fv=NULL;
478 SUMA_DSET *out_dset=NULL, *child=NULL, *dt_dset=NULL;
479 SUMA_CALLBACK *cb= (SUMA_CALLBACK *)params;
480 SUMA_XFORM *xf=NULL;
481 NI_element *nelts = NULL, *nelpars=NULL, *dotopts=NULL;
482 NI_group *ngr = NULL;
483 SUMA_SurfaceObject *SO=NULL;
484 SUMA_OVERLAYS *Sover=NULL;
485 SUMA_Boolean LocalHead = NOPE;
486
487 SUMA_ENTRY;
488
489 if (LocalHead) {
490 SUMA_LH("Callbacks upon entering");
491 SUMA_Show_Callbacks(SUMAg_CF->callbacks, SUMA_STDERR, 1);
492 }
493
494 if (!cb) {
495 SUMA_S_Err("NULL input");
496 SUMA_RETURNe;
497 }
498
499 if (!cb->pending) {
500 SUMA_S_Err("Why am I called?")
501 SUMA_RETURNe;
502 }
503
504 if (strcmp(cb->FunctionName, "SUMA_dot_product_CB")) {
505 SUMA_S_Errv("Baddilicious: %s\n", cb->FunctionName);
506 SUMA_RETURNe;
507 }
508
509 ngr = cb->FunctionInput;
510 if (!ngr) {
511 SUMA_S_Err("NULL input");
512 SUMA_RETURNe;
513 }
514
515 /* get some key elements */
516 nelts = SUMA_FindNgrDataElement(ngr, "callback.data", "ts_vec");
517
518 if (!(nelpars = SUMA_FindNgrNamedElement(ngr, "event_parameters"))) {
519 SUMA_S_Err("parameters element not found!");
520 SUMA_RETURNe;
521 }
522
523 /* Go find out, which xform produced this cb */
524 if (!(xf=SUMA_Find_XformByID(cb->creator_xform))) {
525 SUMA_S_Err("Have no way to reliably select time series dset"
526 "that produced cb");
527 SUMA_RETURNe;
528 }
529
530 if (!(dotopts = SUMA_FindNgrNamedElement(xf->XformOpts, "dotopts"))) {
531 SUMA_S_Err("dotopts element not found!");
532 SUMA_RETURNe;
533 }
534
535 if (LocalHead) {
536 int suc=0;
537 sprintf(stmp,"file:%s.nelpars.1D", FuncName);
538 SUMA_LH("Writing %s", stmp);
539 NEL_WRITE_1D(nelpars, stmp, suc);
540 }
541 if (LocalHead) {
542 int suc=0;
543 sprintf(stmp,"file:%s.dotopts.1D", FuncName);
544 SUMA_LH("Writing %s", stmp);
545 NEL_WRITE_1D(dotopts, stmp, suc);
546 }
547
548 /* get one of the ts_dsets */
549 ts_dset_id = SUMA_GetNgrColStringAttr(ngr, 0, "ts_dsets_idcode");
550 if (!(SUMA_is_ID_4_DSET(ts_dset_id, &in_dset))) {
551 SUMA_S_Err("Could not find ts dset");
552 SUMA_RETURNe;
553 }
554 if (ts_dset_id) SUMA_free(ts_dset_id); ts_dset_id=NULL;
555
556 /* See if ts_node is specified */
557 NI_GET_INT(nelpars, "event.new_node", ts_node);
558 if (!NI_GOT) ts_node = -1;
559
560 /* get yer time series kid */
561 N_ts = -1;
562 if (nelts) {
563 if (nelts->vec_len && nelts->vec_len != SDSET_VECNUM(in_dset)) {
564 SUMA_S_Errv( "ts_vec exists, but is of different length (%d) \n"
565 " than in_dset (%d)\n",
566 nelts->vec_len, SDSET_VECNUM(in_dset));
567 SUMA_RETURNe;
568 } else if (nelts->vec_len) {
569 if (!(SUMA_is_TimeSeries_dset(in_dset, &TR))) {
570 TR = 0.0;
571 }
572 if (ts_node >= 0) {
573 SUMA_LH("Have ts_node AND ts. ts takes precedence\n");
574 }
575 N_ts = nelts->vec_len;
576 /* OK, now get your ts */
577 fv = (float *)nelts->vec[0];
578 if (!(ts = SUMA_DotPreProcessTimeSeries(fv, N_ts,(float)TR, dotopts))) {
579 SUMA_S_Err("Failed processing time series");
580 SUMA_RETURNe;
581 }
582 }
583 }
584
585 if (N_ts < 0) { /* try to get ts from a dset, based on where click happened */
586 if (!(SO = SUMA_findSOp_inDOv(NI_get_attribute(nelpars,"event.DO_idcode"),
587 SUMAg_DOv, SUMAg_N_DOv))) {
588 SUMA_S_Notev("Could not find event's SO (%s)\n"
589 "Use Shift+Ctrl+Right Click on the surface "
590 "to calculate correlation with new changes\n",
591 NI_get_attribute(nelpars,"event.DO_idcode"));
592 /* SUMA_DUMP_TRACE("event.DO_idcode"); */
593 SUMA_RETURNe;
594 }
595 /* find out which overlay was clicked on */
596 if (!(Sover = SUMA_Fetch_OverlayPointer((SUMA_ALL_DO*)SO,
597 NI_get_attribute(nelpars,"event.overlay_name"),&jj))) {
598 SUMA_S_Err("Could not find event's overlay");
599 SUMA_RETURNe;
600 }
601
602 /* Get the dset that corresponds to this overlay */
603 child = Sover->dset_link;
604
605 /* find child dset */
606 if (!SUMA_is_XformChild(xf, SDSET_ID(child), &jj)) {
607 SUMA_S_Err("Failed to find child in xform");
608 SUMA_RETURNe;
609 }
610
611 /* here is the relevant parent dset */
612 if (!(SUMA_is_ID_4_DSET(xf->parents[jj],
613 &ts_src_dset))) {
614 SUMA_S_Err("Could not find ts source dset");
615 SUMA_RETURNe;
616 }
617 if (SDSET_VECNUM(in_dset) != SDSET_VECNUM(ts_src_dset)) {
618 SUMA_S_Errv("Mismatch in ts length (%d vs %d)\n",
619 SDSET_VECNUM(in_dset), SDSET_VECNUM(ts_src_dset));
620 }
621
622 if (NI_YES_ATTR(dotopts,"normalize_dset")) {
623 SUMA_LHv("Getting time series from node %d on detrended version"
624 "of dset %s\n",
625 ts_node, SDSET_FILENAME(ts_src_dset));
626 dt_dset = SUMA_GetDotPreprocessedDset(ts_src_dset, dotopts);
627 } else {
628 SUMA_LHv("Getting time series from node %d of dset %s\n"
629 " (normalize_dset=%s)\n",
630 ts_node, SDSET_FILENAME(ts_src_dset),
631 CHECK_NULL_STR(NI_get_attribute(dotopts,"normalize_dset")));
632 dt_dset = ts_src_dset;
633 }
634 /* get the time series */
635 if (!(ts = (double*)SUMA_GetDsetAllNodeValsInCols2(dt_dset,
636 NULL, 0,
637 ts_node, -1,
638 &N_ts,
639 SUMA_double))) {
640 SUMA_S_Err("Failed to extract time series.");
641 SUMA_RETURNe;
642 }
643 }
644 if (N_ts < 0) {
645 SUMA_S_Err("Nothing to work with here!");
646 SUMA_RETURNe;
647 }
648
649 /* call the function for each of the parents */
650 for (ii=0; ii<cb->N_parents; ++ii) {
651 if (!(out_dset = SUMA_FindDset_s(cb->parents[ii], SUMAg_CF->DsetList))) {
652 SUMA_S_Err("Failed to find out_dset");
653 SUMA_RETURNe;
654 }
655 ts_dset_id = SUMA_GetNgrColStringAttr(ngr, ii, "ts_dsets_idcode");
656 if (!(SUMA_is_ID_4_DSET(ts_dset_id, &in_dset))) {
657 SUMA_S_Err("Could not find ts dset");
658 SUMA_RETURNe;
659 }
660 if (ts_dset_id) SUMA_free(ts_dset_id); ts_dset_id=NULL;
661
662 if (LocalHead) {
663 char stmp[1000];
664 if (!SDSET_LABEL(in_dset)) SUMA_LabelDset(in_dset, NULL);
665 sprintf(stmp,"%s.ts.%d.1D", SDSET_LABEL(in_dset), ts_node);
666 SUMA_LH("Writing %s", stmp);
667 SUMA_WRITE_ARRAY_1D(ts, N_ts, 1, stmp);
668 }
669 SUMA_LHv(" Calculating dot for %d/%d: \n"
670 " Time series dset: %s\n"
671 " Dot product output in: %s\n",
672 ii, cb->N_parents,
673 SDSET_FILENAME(in_dset),
674 SDSET_FILENAME(out_dset) );
675 /* call dot product computer */
676 if (!SUMA_dot_product(in_dset,ts,&out_dset,dotopts)) {
677 SUMA_S_Err("Failed to compute dot product");
678 SUMA_RETURNe;
679 }
680
681 SUMA_LH("Now clearing overlay vecs");
682 /* Make sure color overlay datacopies get reset */
683 if (!SUMA_DSET_ClearOverlay_Vecs(out_dset)) {
684 SUMA_S_Err("Failed to clear overlay copies");
685 SUMA_RETURNe;
686 }
687
688 snprintf(ident,298*sizeof(char), "filename:%s", SDSET_FILENAME(out_dset));
689
690 SUMA_LH("Adding to save list");
691 {
692 SUMA_PARSED_NAME *p1p, *p2p;
693 p1p = SUMA_ParseFname_eng(SDSET_FILENAME(in_dset), SUMAg_CF->cwd, 0);
694 p1 = SUMA_RemoveDsetExtension_s(p1p->FileName_NoExt,
695 SUMA_NO_DSET_FORMAT);
696 p2p = SUMA_ParseFname_eng(SDSET_FILENAME(ts_src_dset), SUMAg_CF->cwd, 0);
697 p2 = SUMA_RemoveDsetExtension_s(p2p->FileName_NoExt,
698 SUMA_NO_DSET_FORMAT);
699 SUMA_Free_Parsed_Name(p1p); p1p = NULL;
700 SUMA_Free_Parsed_Name(p2p); p2p = NULL;
701 }
702
703 if (SO->Side == SUMA_LEFT) Cside = "L";
704 else if (SO->Side == SUMA_RIGHT) Cside = "R";
705 else Cside = "";
706 snprintf(prefix,298*sizeof(char),"seed_%d%s.DOT.%s",
707 ts_node, Cside, p1);
708 SUMA_free(p1); SUMA_free(p2);
709 if (!(SUMA_Add_to_SaveList(&SUMAg_CF->SaveList, "sdset", ident, prefix))) {
710 SUMA_S_Warnv("Failed to add to save list %s %s\n", ident, prefix);
711 }
712 }
713
714 if (ts) SUMA_free(ts); ts = NULL;
715
716 /* Now clear callback specific elements for caution
717 This way I'll know if this function
718 is being called repeatedly by mistake*/
719 SUMA_LH("Clearing callback specific elements");
720 if (nelts) {
721 NI_remove_from_group(ngr, nelts); NI_free_element(nelts); nelts = NULL;
722 }
723
724
725 SUMA_RETURNe;
726 }
727
SUMA_DotDetrendDset(SUMA_DSET * in_dset,float ** refvec,int nref,float fbot,float ftop,int qdet,int * num_ort)728 SUMA_DSET *SUMA_DotDetrendDset( SUMA_DSET *in_dset,
729 float **refvec, int nref,
730 float fbot, float ftop,
731 int qdet, int *num_ort)
732 {
733 static char FuncName[]={"SUMA_DotDetrendDset"};
734 float **fvec=NULL;
735 double TR=0;
736 int i, N_ret=0, nnort = 0;
737 SUMA_DSET *o_dset=NULL;
738 SUMA_Boolean LocalHead = NOPE;
739
740 SUMA_ENTRY;
741
742 if (!refvec || !nref || !in_dset) SUMA_RETURN(NULL);
743
744 if (!(SUMA_is_TimeSeries_dset(in_dset, &TR))) {
745 TR = 0.0;
746 }
747
748 /* turn the dset to a float vector array */
749 if (!(fvec = (float **)SUMA_Dset2VecArray( in_dset,
750 NULL, -1, /* all columns */
751 NULL, -1, /* all nodes */
752 -1, /* don't care to figure out max node */
753 &N_ret,
754 SUMA_float /* float me back */))) {
755 SUMA_S_Err("Failed to copy surface dset");
756 SUMA_RETURN(NULL);
757 }
758
759 SUMA_LH("Bandpass/detrending phase");
760 /* detrend */
761 nnort = THD_bandpass_vectors (SDSET_VECNUM(in_dset), SDSET_VECLEN(in_dset),
762 fvec, (float)TR, fbot, ftop, qdet, nref,
763 refvec);
764 if (!nnort) {
765 SUMA_S_Err("Bad bandpass call, going to hell now.");
766 SUMA_RETURN(NULL);
767 }
768 if (num_ort) *num_ort = nnort;
769
770 /* normalize */
771 for (i=0; i<SDSET_VECLEN(in_dset); ++i) {
772 THD_normalize( SDSET_VECNUM(in_dset) , fvec[i] ) ;
773 }
774
775 SUMA_LH("About to form output dset");
776 /* make a copy of the input dset */
777 o_dset = SUMA_MaskedCopyofDset(in_dset, NULL, NULL, 1, 0);
778
779 SUMA_LH("Now filling with output");
780 /* Now fill it with fvec*/
781 if (!SUMA_VecArray2Dset((void **)fvec,
782 o_dset,
783 NULL, -1,
784 NULL, -1,
785 -1,
786 SUMA_float)) {
787 SUMA_S_Err("Misery");
788 SUMA_RETURN(NULL);
789 }
790
791 /* cleanup the large array */
792 for (i=0; i<SDSET_VECLEN(in_dset); ++i) {
793 SUMA_free(fvec[i]); fvec[i]=NULL;
794 }
795 SUMA_free(fvec); fvec = NULL;
796
797 /* write out dset for good measure? */
798 if (LocalHead) {
799 SUMA_LH("Writing detrended dset to disk\n");
800 SUMA_WriteDset_s(FuncName, o_dset, SUMA_ASCII_NIML, 1, 1);
801 }
802
803 /* that is it, return */
804 SUMA_RETURN(o_dset);
805 }
806
807
SUMA_dot_product(SUMA_DSET * in_dset,double * ts,SUMA_DSET ** out_dsetp,NI_element * dotopt)808 SUMA_Boolean SUMA_dot_product(SUMA_DSET *in_dset,
809 double *ts,
810 SUMA_DSET **out_dsetp,
811 NI_element *dotopt)
812 {
813 static char FuncName[]={"SUMA_dot_product"};
814 SUMA_DSET *dot = NULL;
815 static SUMA_DSET *in_dset_last=NULL;
816 double *dcol=NULL;
817 float *fcol = NULL;
818 int ic=0, ir=0;
819 byte *bbv=NULL;
820 int *iiv=NULL;
821 float *ffv=NULL;
822 double *ddv=NULL;
823 int prec=1; /* NEVER CHANGE THIS DEFAULT */
824 char *s=NULL, *sname=NULL;
825 SUMA_VARTYPE vtp = SUMA_notypeset;
826 float par[3];
827 static SUMA_VARTYPE vtp_last = SUMA_notypeset;
828 NI_element *nelb=NULL;
829 SUMA_Boolean LocalHead = NOPE;
830
831 SUMA_ENTRY;
832
833 if ( !in_dset ||
834 !out_dsetp ||
835 !ts ) {
836 SUMA_S_Err("NULL input");
837 SUMA_RETURN(NOPE);
838 }
839
840 if (in_dset != in_dset_last) {
841 SUMA_LH("Checking consistency");
842 if (!SUMA_is_AllConsistentNumeric_dset(in_dset, &vtp)) {
843 SUMA_S_Err( "Input dataset is either not all numeric \n"
844 "or has inconsistent types");
845 SUMA_RETURN(NOPE);
846 }
847 in_dset_last = in_dset;
848 vtp_last = vtp;
849 } else {
850 vtp = vtp_last;
851 }
852
853 /* decide if we need new_dset or not */
854 if (*out_dsetp) {
855 SUMA_LH("Reusing dset");
856 dot = *out_dsetp;
857 if (SUMA_ColType2TypeCast(SUMA_TypeOfDsetColNumb(dot, 0))==SUMA_double) {
858 prec = 2;
859 dcol = (double*)dot->dnel->vec[0];
860 memset((void *)dcol, 0, SDSET_VECLEN(dot)*sizeof(double));
861 } else {
862 prec = 1;
863 fcol = (float *)dot->dnel->vec[0];
864 memset((void *)fcol, 0, SDSET_VECLEN(dot)*sizeof(float));
865 }
866 } else {
867 if (dotopt) {
868 NI_GET_INT(dotopt,"numeric_precision", prec);
869 if (!NI_GOT) prec = 1;
870 } else {
871 /* NEVER CHANGE THIS , prec is set at initialization */
872 }
873 if (!SDSET_LABEL(in_dset)) SUMA_LabelDset(in_dset, NULL);
874
875 sname = SUMA_append_string("dot.",SDSET_FILENAME(in_dset));
876 SUMA_LHv("Creating new dset (%s) with precision %d\n", sname, prec);
877 if (!(dot = SUMA_CreateDsetPointer(
878 sname, SUMA_NODE_BUCKET, NULL,
879 SDSET_IDMDOM(in_dset),
880 SDSET_VECLEN(in_dset) ) ) ){
881
882 }
883 SUMA_free(sname); sname=NULL;
884 SUMA_LabelDset(dot, NULL);
885
886 if (!SUMA_AddDsetNelCol (dot, NI_get_attribute(in_dset->inel,"COLMS_LABS"),
887 SUMA_NODE_INDEX,
888 (void *)SDSET_NODE_INDEX_COL(in_dset),
889 NULL,1 )) {
890 fprintf (stderr,"Error %s:\nFailed in SUMA_AddNelCol", FuncName);
891 exit(1);
892 }
893 if (prec == 2) {
894 dcol = (double *) SUMA_calloc(SDSET_VECLEN(in_dset), sizeof(double));
895 if (!SUMA_AddDsetNelCol(dot, "dotproduct", SUMA_NODE_DOUBLE,
896 (void *)dcol, NULL , 1)) {
897 SUMA_S_Err("Failed to add col");
898 SUMA_RETURN(NOPE);
899 }
900 SUMA_free(dcol); dcol = (double*)dot->dnel->vec[0];
901 } else {
902 fcol = (float *) SUMA_calloc(SDSET_VECLEN(in_dset), sizeof(float));
903 if (dotopt && NI_YES_ATTR(dotopt,"normalize_dset")) {
904 /* par[2] is not available yet. The 2 is just a place holder */
905 par[0] = SDSET_VECNUM(in_dset); par[1] = 1; par[2]=2;
906 if (!SUMA_AddDsetNelCol(dot, "XcorrCoef", SUMA_NODE_XCORR,
907 (void *)fcol, (void *)par , 1)) {
908 SUMA_S_Err("Failed to add col");
909 SUMA_RETURN(NOPE);
910 }
911 } else {
912 if (!SUMA_AddDsetNelCol(dot, "dotproduct", SUMA_NODE_FLOAT,
913 (void *)fcol, NULL , 1)) {
914 SUMA_S_Err("Failed to add col");
915 SUMA_RETURN(NOPE);
916 }
917 }
918 SUMA_free(fcol); fcol = (float*)dot->dnel->vec[0];
919 }
920 *out_dsetp = dot;
921 }
922
923 /* preprocess in_dset ? */
924 if (dotopt) {
925 if (NI_YES_ATTR(dotopt,"normalize_dset")) {
926 SUMA_LH("Dotting preprocessed dset");
927 /* WARNING: From this point on in_dset will point
928 to another dset. */
929 if (!(in_dset = SUMA_GetDotPreprocessedDset(in_dset, dotopt))) {
930 SUMA_S_Err("Failed to get preprocessed dset");
931 SUMA_RETURN(NOPE);
932 }
933 } else {
934 SUMA_LH("Dotting initial dset");
935 }
936 } else {
937 SUMA_LH("Dotting initial dset, no dotopt");
938 }
939
940 SUMA_LHv("Dotting %dx%d values at x%d precision...\n",
941 SDSET_VECNUM(in_dset), SDSET_VECLEN(in_dset), prec);
942 switch(vtp) {
943 case SUMA_byte:
944 for (ic=0; ic<SDSET_VECNUM(in_dset); ++ic) {
945 bbv = (byte*)in_dset->dnel->vec[ic];
946 if (prec == 2) {
947 for (ir=0; ir<SDSET_VECLEN(in_dset); ++ir)
948 dcol[ir] += ts[ic]*(double)bbv[ir];
949 } else {
950 for (ir=0; ir<SDSET_VECLEN(in_dset); ++ir)
951 fcol[ir] += (float)(ts[ic]*(double)bbv[ir]);
952 }
953 }
954 break;
955 case SUMA_int:
956 for (ic=0; ic<SDSET_VECNUM(in_dset); ++ic) {
957 iiv = (int*)in_dset->dnel->vec[ic];
958 if (prec == 2) {
959 for (ir=0; ir<SDSET_VECLEN(in_dset); ++ir)
960 dcol[ir] += ts[ic]*(double)iiv[ir];
961 } else {
962 for (ir=0; ir<SDSET_VECLEN(in_dset); ++ir)
963 fcol[ir] += (float)(ts[ic]*(double)iiv[ir]);
964 }
965 }
966 break;
967 case SUMA_float:
968 for (ic=0; ic<SDSET_VECNUM(in_dset); ++ic) {
969 ffv = (float*)in_dset->dnel->vec[ic];
970 if (prec == 2) {
971 for (ir=0; ir<SDSET_VECLEN(in_dset); ++ir)
972 dcol[ir] += ts[ic]*(double)ffv[ir];
973 } else {
974 for (ir=0; ir<SDSET_VECLEN(in_dset); ++ir)
975 fcol[ir] += (float)(ts[ic]*(double)ffv[ir]);
976 }
977 }
978 break;
979 case SUMA_double:
980 for (ic=0; ic<SDSET_VECNUM(in_dset); ++ic) {
981 ddv = (double*)in_dset->dnel->vec[ic];
982 if (prec == 2) {
983 for (ir=0; ir<SDSET_VECLEN(in_dset); ++ir)
984 dcol[ir] += ts[ic]*ddv[ir];
985 } else {
986 for (ir=0; ir<SDSET_VECLEN(in_dset); ++ir)
987 fcol[ir] += (float)(ts[ic]*ddv[ir]);
988 }
989 }
990 break;
991 case SUMA_complex:
992 SUMA_S_Err("No support for complex type here");
993 SUMA_RETURN(NOPE);
994 break;
995 default:
996 SUMA_S_Err("What kind of numeric type is this?");
997 SUMA_RETURN(NOPE);
998 break;
999 }
1000
1001 SUMA_LH("Updating range\n");
1002
1003 if(!(s = SUMA_CreateDsetColRangeCompString(dot, 0,
1004 SUMA_TypeOfDsetColNumb(dot,0)))) {
1005 SUMA_S_Err("Failed to calculate range");
1006 SUMA_RETURN(NOPE);
1007 }else {
1008 NI_element *nelb = SUMA_FindDsetAttributeElement(dot, "COLMS_RANGE");
1009 SUMA_LHv("New range string is %s\n", s);
1010 SUMA_AddColAtt_CompString(nelb, 0, s, SUMA_NI_CSS,0);
1011 SUMA_free(s); s=NULL;
1012 }
1013
1014 /* reset the stat params, just in case */
1015 if (dotopt && NI_YES_ATTR(dotopt,"normalize_dset") &&
1016 (nelb = SUMA_FindDsetAttributeElement(dot, "COLMS_STATSYM"))) {
1017 par[0] = SDSET_VECNUM(in_dset); par[1] = 1;
1018 NI_GET_INT(dotopt,"num_ort_parameters", par[2]);
1019 if (!NI_GOT || par[2] == -1.0) {
1020 SUMA_S_Errv("Failed to get ort parameters (%f), \n"
1021 "you can't trust the p values\n ",
1022 par[2] );
1023 par[2] = 2; /* something ... */
1024 }
1025 SUMA_LH("Adding ort stuff");
1026 SUMA_AddColAtt_CompString( nelb, 0,
1027 NI_stat_encode(NI_STAT_CORREL,
1028 par[0], par[1], par[2]),
1029 SUMA_NI_CSS, 0);
1030 }
1031 SUMA_LH("Heading out");
1032 if (0 && LocalHead) {
1033 ic = 0;
1034 for (ir=0; ir<SDSET_VECLEN(in_dset) && ic < 10; ++ir) {
1035 if (prec ==2) {
1036 if (dcol[ir]) {
1037 fprintf(SUMA_STDERR,"dcol[%d]=%f\n", ir, dcol[ir]);
1038 ++ic;
1039 }
1040 } else {
1041 if (fcol[ir]) {
1042 fprintf(SUMA_STDERR,"fcol[%d]=%f\n", ir, fcol[ir]);
1043 ++ic;
1044 }
1045 }
1046 }
1047 if (!ic) {
1048 SUMA_S_Note("No non zero rows encountered!");
1049 }
1050 ir = 1;ic = 0;
1051 DSET_WRITE_1D(dot, "file:dot.1D.dset", ic, ir);
1052 s = SUMA_WriteDset_s ("dot", dot, SUMA_ASCII_NIML, 1, 1);
1053 SUMA_free(s);
1054 }
1055
1056 SUMA_RETURN(YUP);
1057 }
1058
1059 /* *************************************************************************** */
1060 /* ******************* Functions to deal with Group InstaCorr **************** */
1061 /* *************************************************************************** */
1062
1063 /*! find dsets and their overlay planes that go with giset
1064 if target_name is not null, the function creates new dsets and their
1065 corresponding overlay planes, otherwise, it returns what had been created.
1066 */
SUMA_GICOR_Dsets(SUMA_SurfaceObject * SOv[],GICOR_setup * giset,char * target_name,DList * DsetList,SUMA_DSET * sdsetv[],SUMA_OVERLAYS * ov[],NI_element * nel)1067 SUMA_Boolean SUMA_GICOR_Dsets(SUMA_SurfaceObject *SOv[],
1068 GICOR_setup *giset,
1069 char *target_name,
1070 DList *DsetList,
1071 SUMA_DSET *sdsetv[],
1072 SUMA_OVERLAYS *ov[],
1073 NI_element *nel)
1074 {
1075 static char FuncName[]={"SUMA_GICOR_Dsets"};
1076 char *targetv[2]={NULL, NULL},
1077 *dset_namev[2]={NULL, NULL}, *atr=NULL;
1078 int i, ii, *Ti=NULL, ovind = 0, nvals=0, vv=0;
1079 static char *blab[6] = { "GIC_Delta" , "GIC_Zscore" ,
1080 "AAA_Delta" , "AAA_Zscore" ,
1081 "BBB_Delta" , "BBB_Zscore" } ;
1082 NI_str_array *labar=NULL ;
1083 SUMA_Boolean LocalHead = NOPE;
1084
1085 SUMA_ENTRY;
1086
1087 if (target_name &&
1088 giset->sdset_ID[0][0] != '\0' && giset->sdset_ID[1][0] != '\0') {
1089 SUMA_S_Warn("Hello anew from 3dGroupInCorr, \n"
1090 "Attempting to reuse previous setup...");
1091 goto CHECK_DSET_AND_OVERLAYS;
1092 }
1093
1094 if (target_name) { /* Brand new init, search/create by name */
1095 SUMA_LHv("Brand new init, target_name=%s\n", target_name);
1096 if (!nel) {
1097 SUMA_S_Err("Need GICOR setup nel for creating new dsets");
1098 SUMA_RETURN(NOPE);
1099 }
1100 /* Form the names of the dsets to be created */
1101 if (SOv[1]) { /* have two surfaces */
1102 targetv[0] = SUMA_append_string(target_name,".Left");
1103 targetv[1] = SUMA_append_string(target_name,".Right");
1104 } else {
1105 targetv[0] = SUMA_copy_string(target_name);
1106 }
1107
1108 /* Now create a dataset for each case */
1109 for (i=0; i<2; ++i) {
1110 if (targetv[i]) {
1111 SUMA_LHv("Working %s\n", targetv[i]);
1112 /* dset names */
1113 dset_namev[i] = SUMA_append_string(targetv[i], SOv[i]->idcode_str);
1114 sdsetv[i] = SUMA_CreateDsetPointer (dset_namev[i],
1115 SUMA_NODE_BUCKET,
1116 NULL,
1117 SOv[i]->idcode_str,
1118 SOv[i]->N_Node);
1119 sprintf(giset->sdset_ID[i],"%s", SDSET_ID(sdsetv[i]));
1120
1121 /* insert that element into DaList */
1122 if (!SUMA_InsertDsetPointer(&sdsetv[i], DsetList, 0)) {
1123 SUMA_SL_Err("Failed to insert dset into list");
1124 SUMA_free(dset_namev[i]); SUMA_free(targetv[i]);
1125 SUMA_RETURN(NOPE);
1126 }
1127 /* add the columns */
1128 Ti = (int *) SUMA_calloc(SDSET_VECLEN(sdsetv[i]), sizeof(int));
1129 for (ii=0; ii <SDSET_VECLEN(sdsetv[i]); ++ii) Ti[ii]=ii;
1130 SUMA_AddDsetNelCol (sdsetv[i], "node index",
1131 SUMA_NODE_INDEX, Ti, NULL, 1);
1132 SUMA_free(Ti); Ti=NULL;
1133
1134 atr = NI_get_attribute( nel , "target_nvals" ) ;
1135 if( atr == NULL ) SUMA_GIQUIT;
1136 nvals = (int)strtod(atr,NULL) ; nvals = MAX(1,nvals);
1137
1138 atr = NI_get_attribute( nel , "target_labels" ) ;
1139 if( atr != NULL )
1140 labar = NI_decode_string_list( atr , ";" ) ;
1141
1142 for( vv=0 ; vv < nvals ; vv++ ){
1143 if (labar != NULL && vv < labar->num) atr = labar->str[vv];
1144 else if (vv < 6) {
1145 atr = blab[vv];
1146 } else {
1147 atr = "What the hell is this?";
1148 }
1149 if (vv%2 == 0) { /* beta */
1150 SUMA_AddDsetNelCol (sdsetv[i], atr,
1151 SUMA_NODE_FLOAT, NULL, NULL, 1);
1152 } else { /* zscore */
1153 SUMA_AddDsetNelCol (sdsetv[i], atr,
1154 SUMA_NODE_ZSCORE, NULL, NULL, 1);
1155 }
1156 }
1157 if (labar) SUMA_free_NI_str_array(labar); labar=NULL;
1158
1159 SUMA_LHv("Creating overlays for %s\n", targetv[i]);
1160 /* create overlays */
1161 ov[i] = SUMA_CreateOverlayPointer (targetv[i], sdsetv[i],
1162 SOv[i]->idcode_str, NULL);
1163 if (!ov[i]) {
1164 SUMA_SL_Err("Failed in SUMA_CreateOverlayPointer.\n");
1165 SUMA_free(dset_namev[i]); SUMA_free(targetv[i]);
1166 SUMA_RETURN(NOPE);
1167 }
1168 ov[i]->ShowMode = SW_SurfCont_DsetViewCol;
1169 ov[i]->GlobalOpacity = 0.8;
1170 ov[i]->isBackGrnd = NOPE;
1171 ov[i]->OptScl->BrightFact = 0.5;
1172 ov[i]->OptScl->find = 0;
1173 ov[i]->OptScl->tind = 1;
1174 ov[i]->OptScl->bind = 0;
1175 ov[i]->OptScl->UseThr = 1; /* turn on threshold use */
1176 ov[i]->SymIrange = 1; /* Use symmetric range */
1177 ov[i]->OptScl->AutoIntRange = 0; /* Do not update range */
1178 ov[i]->OptScl->IntRange[0] = -0.5; /* set the range */
1179 ov[i]->OptScl->IntRange[1] = 0.5;
1180 ov[i]->OptScl->ThreshRange[0] = 2.0;
1181 ov[i]->OptScl->ThreshRange[1] = 0.0;
1182
1183 /* Now add the overlay to SOv[i]->Overlays */
1184 if (!SUMA_AddNewPlane ((SUMA_ALL_DO *)SOv[i], ov[i],
1185 SUMAg_DOv, SUMAg_N_DOv, 1)) {
1186 SUMA_SL_Crit("Failed in SUMA_AddNewPlane");
1187 SUMA_FreeOverlayPointer(ov[i]);
1188 SUMA_free(dset_namev[i]); SUMA_free(targetv[i]);
1189 SUMA_RETURN (NOPE);
1190 }
1191 SUMA_free(dset_namev[i]); dset_namev[i]=NULL;
1192 SUMA_free(targetv[i]); targetv[i]=NULL;
1193 }
1194 }
1195
1196 /* Done with brand new init */
1197 SUMA_RETURN(YUP);
1198 }
1199
1200 CHECK_DSET_AND_OVERLAYS:
1201 SUMA_LH("Checking Dset and Overlays.\n");
1202
1203 { /* just use what is in giset */
1204 if (giset->sdset_ID[0][0] == '\0') {
1205 SUMA_S_Err("No ID in sdset_ID. Unexpected happenstance");
1206 SUMA_RETURN(NOPE);
1207 }
1208 if (!(sdsetv[0] = SUMA_FindDset_s(giset->sdset_ID[0], DsetList))) {
1209 SUMA_S_Err("SDSET for 0 not found");
1210 SUMA_RETURN(NOPE);
1211 }
1212 if (giset->sdset_ID[1][0] != '\0') {
1213 if (!(sdsetv[1] = SUMA_FindDset_s(giset->sdset_ID[1], DsetList))) {
1214 SUMA_S_Err("SDSET for 1 not found");
1215 SUMA_RETURN(NOPE);
1216 }
1217 }
1218 /* fetch the overlays */
1219 for (i=0; i<2; ++i) {
1220 if (sdsetv[i]) {
1221 if (!(ov[i] = SUMA_Fetch_OverlayPointerByDset(
1222 (SUMA_ALL_DO *)SOv[i],
1223 sdsetv[i], &ovind))) {
1224 SUMA_S_Err("Failed to find overlay pointer");
1225 SUMA_RETURN(NOPE);
1226 }
1227 }
1228 }
1229 }
1230
1231 /* at this point, we have the relevant dsets in sdsetv,
1232 and their overlays in ov */
1233 SUMA_RETURN(YUP);
1234 }
1235
1236 /*! find surfaces appropriate for giset */
SUMA_GICOR_Surfaces(GICOR_setup * giset,SUMA_SurfaceObject * SOv[])1237 SUMA_Boolean SUMA_GICOR_Surfaces(GICOR_setup *giset, SUMA_SurfaceObject *SOv[])
1238 {
1239 static char FuncName[]={"SUMA_GICOR_Surfaces"};
1240
1241 SUMA_ENTRY;
1242
1243 if (!(SOv[0] = SUMA_FindSOp_inDOv_from_N_Node(
1244 giset->nnode_domain[0],
1245 giset->nnode_domain[1] ? SUMA_LEFT:SUMA_NO_SIDE,
1246 1, 1,
1247 SUMAg_DOv, SUMAg_N_DOv))) {
1248 SUMA_S_Errv("Could not find domain parent for a domain of %d nodes\n",
1249 giset->nnode_domain[0]);
1250 SUMA_RETURN(NOPE);
1251 }
1252
1253 if (giset->nnode_domain[1]) {
1254 if (!(SOv[1]=SUMA_FindSOp_inDOv_from_N_Node(
1255 giset->nnode_domain[1], SUMA_RIGHT,
1256 1, 1,
1257 SUMAg_DOv, SUMAg_N_DOv))) {
1258 SUMA_S_Errv("Could not find domain parent for a "
1259 "RH domain of %d nodes\n",
1260 giset->nnode_domain[1]);
1261 SUMA_RETURN(NOPE);
1262 }
1263 }
1264
1265 SUMA_RETURN(YUP);
1266 }
1267
1268
1269 /*!
1270 Function called from SUMA_niml, when GICorr sends a setup element
1271 This function parallels AFNI's GICOR_setup_func
1272 */
SUMA_GICOR_setup_func(NI_stream nsg,NI_element * nel)1273 SUMA_Boolean SUMA_GICOR_setup_func( NI_stream nsg , NI_element *nel )
1274 {
1275 static char FuncName[]={"SUMA_GICOR_setup_func"};
1276 GICOR_setup *giset = NULL;
1277 SUMA_DSET *sdsetv[2]={NULL, NULL};
1278 SUMA_OVERLAYS *ov[2]={NULL, NULL};
1279 SUMA_SurfaceObject *SOv[2]={NULL, NULL};
1280 char *pre=NULL;
1281 SUMA_Boolean LocalHead = NOPE;
1282
1283 SUMA_ENTRY;
1284
1285
1286 /* fetch the giset struct */
1287 giset = SUMAg_CF->giset;
1288 if( giset != NULL && giset->ready ) SUMA_RETURN(YUP) ;
1289
1290 if( giset == NULL ){
1291 giset = (GICOR_setup *)SUMA_calloc(1,sizeof(GICOR_setup)) ;
1292 SUMAg_CF->giset = giset;
1293 } else {
1294 memset(giset,0, sizeof(GICOR_setup)) ;
1295 }
1296
1297 if (!SUMA_init_GISET_setup(nsg, nel, giset, 0)) SUMA_GIQUIT;
1298
1299 /* Now find surfaces that can be the domain */
1300 if (!SUMA_GICOR_Surfaces(giset, SOv)) {
1301 SUMA_S_Err("Failed to find surfaces for giset");
1302 SUMA_RETURN(NOPE);
1303 }
1304
1305 /* Now create appropriate dsets */
1306 pre = NI_get_attribute( nel , "target_name" ) ;
1307 if( pre == NULL || *pre == '\0' ) pre = "GICorrelletto" ;
1308 if (!SUMA_GICOR_Dsets(SOv, giset, pre, SUMAg_CF->DsetList,
1309 sdsetv, ov, nel)) {
1310 SUMA_S_Err("Failed to find/create dsets for giset");
1311 SUMA_RETURN(NOPE);
1312 }
1313
1314
1315 giset->ready = 1 ;
1316
1317 if (LocalHead) {
1318 SUMA_Show_GISET(giset, NULL, 0);
1319 }
1320
1321 SUMA_RETURN(YUP) ;
1322 }
1323
1324 /*! Surface version of AFNI's GICOR_process_dataset*/
SUMA_GICOR_process_dataset(NI_element * nel)1325 SUMA_Boolean SUMA_GICOR_process_dataset( NI_element *nel )
1326 {
1327 static char FuncName[]={"SUMA_GICOR_process_dataset"};
1328 GICOR_setup *giset = SUMAg_CF->giset ;
1329 int vmul ; float thr ;
1330 int id=0, ic=0;
1331 char *seedstr=NULL, *Cside=NULL, ident[300]={""}, prefix[300]={""};
1332 SUMA_SurfaceObject *SOv[2]={NULL,NULL};
1333 SUMA_DSET *sdsetv[2]={NULL, NULL};
1334 SUMA_OVERLAYS *ov[2]={NULL, NULL};
1335 SUMA_Boolean LocalHead = NOPE;
1336
1337 SUMA_ENTRY;
1338
1339 if( nel == NULL || nel->vec_num < 2 ){ /* should never happen */
1340 ERROR_message("badly formatted dataset from 3dGroupInCorr!") ;
1341 SUMA_RETURN(NOPE) ;
1342 }
1343
1344 if (nel->vec_num % 2) {
1345 SUMA_SLP_Err("Number of sub-bricks not multiple of two!");
1346 SUMA_RETURN(NOPE) ;
1347 }
1348
1349 if( giset == NULL ||
1350 !giset->ready ){ /* should not happen */
1351
1352 if( giset != NULL ) giset->ready = 0 ;
1353 /* AFNI_misc_CB(im3d->vwid->func->gicor_pb,(XtPointer)im3d,NULL) ; */
1354 SUMA_SLP_Err(" ******* SUMA: *********\n"
1355 " 3dGrpInCorr sent data \n"
1356 " but setup isn't ready!\n" ) ;
1357 SUMA_RETURN(NOPE) ;
1358 }
1359
1360 /* get the surfaces, the dsets, and their overlays */
1361 if (!SUMA_GICOR_Surfaces(giset, SOv)) {
1362 SUMA_S_Err("Failed to find surfaces for giset");
1363 SUMA_RETURN(NOPE);
1364 }
1365 if (!SUMA_GICOR_Dsets(SOv, giset, NULL, SUMAg_CF->DsetList,
1366 sdsetv, ov, NULL)) {
1367 SUMA_S_Err("Failed to find/create dsets for giset");
1368 SUMA_RETURN(NOPE);
1369 }
1370
1371 /* copy NIML data into dataset */
1372
1373 SUMA_LH("Populating the dset in question, redisplay, etc.");
1374
1375 if (!SUMA_PopulateDsetsFromGICORnel(nel, giset, sdsetv)) {
1376 SUMA_S_Err("Failed to populate. Not fun.");
1377 SUMA_RETURN(NOPE);
1378 }
1379
1380 /* Add dsets to autosave list */
1381 if (!(seedstr = NI_get_attribute(nel,"boomerang_msg"))) {
1382 seedstr = "seed_XXXX";
1383 }
1384 for (id=0; id<2; ++id) {
1385 if (sdsetv[id] && SOv[id]) {
1386 snprintf(ident,298*sizeof(char), "filename:%s",
1387 SDSET_FILENAME(sdsetv[id]));
1388 if (SOv[id]->Side == SUMA_LEFT) Cside = ".lh";
1389 else if (SOv[id]->Side == SUMA_RIGHT) Cside = ".rh";
1390 else Cside = "";
1391 snprintf(prefix,298*sizeof(char),"seed_%s.GIC%s",
1392 seedstr, Cside);
1393 if (!(SUMA_Add_to_SaveList(&SUMAg_CF->SaveList, "sdset",
1394 ident, prefix))) {
1395 SUMA_S_Warnv("Failed to add to save list %s %s\n", ident, prefix);
1396 }
1397 }
1398 }
1399
1400 /* colorize and redisplay */
1401 for (id=0; id<2; ++id) {
1402 if (ov[id] && SOv[id]) {
1403 SOv[id]->SurfCont->curColPlane = ov[id];
1404 if (!SUMA_OpenCloseSurfaceCont(NULL, (SUMA_ALL_DO*)SOv[id], NULL)) {
1405 SUMA_SLP_Err("Cannot open Surface Controller!");
1406 SUMA_RETURN(NOPE);
1407 }
1408 if ( SOv[id]->SurfCont->SwitchDsetlst &&
1409 !SOv[id]->SurfCont->SwitchDsetlst->isShaded ) {
1410 SUMA_RefreshDsetList ((SUMA_ALL_DO*)SOv[id]);
1411 }
1412 SUMA_LHv("Colorizing %d\n", id);
1413 if (!SUMA_ColorizePlane (ov[id])) {
1414 SUMA_SLP_Err("Failed to colorize plane.\n");
1415 SUMA_RETURN(NOPE);
1416 }
1417 if (!SUMA_Remixedisplay ((SUMA_ALL_DO*)SOv[id])) {
1418 SUMA_SLP_Err("Failed to remix redisplay.\n");
1419 SUMA_RETURN(NOPE);
1420 }
1421 SUMA_LHv("Initializing %d\n", id);
1422 SUMA_UpdateColPlaneShellAsNeeded((SUMA_ALL_DO*)SOv[id]);
1423 }
1424 }
1425
1426 #if 0 /* This is no longer needed because of the modification
1427 of SE_Redisplay*All* . The last controller to be
1428 redrawn is now the one in which the pointer last
1429 resided. I leave this here in case the problem
1430 resurfaces again ... */
1431 /* if you have multiple controllers open, force a redisplay
1432 in the one where the click occurred, to force a reset of
1433 GL's State variables. Otherwise you may have problems
1434 selecting nodes repeatedly without going in and out of the window */
1435 if (SUMAg_CF->PointerLastInViewer >= 0 && SUMAg_N_SVv > 1) {
1436 SUMA_SurfaceViewer *sv = &(SUMAg_SVv[SUMAg_CF->PointerLastInViewer]);
1437 if (!sv->isShaded) { /* this should always be true ... */
1438 SUMA_OpenGLStateReset(SUMAg_DOv, SUMAg_N_DOv, sv);
1439 SUMA_handleRedisplay((XtPointer)sv->X->GLXAREA);
1440 }
1441 }
1442 #endif
1443
1444 if (LocalHead) {
1445 SUMA_ShowDset(sdsetv[0],0,NULL);
1446 SUMA_ShowDset(sdsetv[1],0,NULL);
1447 }
1448
1449 #if 0
1450 /* YOU CANNOT DO THIS: Although meshes are isotopic,
1451 a node index do not necessarily refer to a homologous
1452 anatomical region on both left and right surfaces. */
1453 if (SUMAg_CF->Dev) {
1454 if (SOv[0] && SOv[1] &&
1455 SOv[0]->N_Node == SOv[1]->N_Node) {
1456 static float *fv1=NULL, *fv0=NULL;
1457 static double dot[12];
1458 int ii=0;
1459
1460 if (!fv0) fv0= (float *)SUMA_calloc(SOv[0]->N_Node, sizeof(float));
1461 if (!fv1) fv1= (float *)SUMA_calloc(SOv[1]->N_Node, sizeof(float));
1462 SUMA_S_Note ( "Cross correlation of left with right "
1463 "hemisphere patterns:\n");
1464 for (ipair=0; ipair < nel->vec_num/2; ++ipair) {
1465 memcpy(fv0,(SDSET_VEC(sdsetv[0],ipair)),
1466 sizeof(float)*SDSET_VECLEN(sdsetv[0])) ;
1467 memcpy(fv1,(SDSET_VEC(sdsetv[1],ipair)),
1468 sizeof(float)*SDSET_VECLEN(sdsetv[1])) ;
1469 THD_normalize( SDSET_VECLEN(sdsetv[0]) , fv0 ) ;
1470 THD_normalize( SDSET_VECLEN(sdsetv[1]) , fv1 ) ;
1471 dot[ipair]=0.0;
1472 for (ii=0; ii<SDSET_VECLEN(sdsetv[0]); ++ii)
1473 dot[ipair] += (double)fv0[ii]*(double)fv1[ii];
1474 fprintf(SUMA_STDERR, "Pair %d: %f\t", ipair, dot[ipair]);
1475 }
1476 fprintf(SUMA_STDERR, "\n");
1477 }
1478 }
1479 #endif
1480
1481 SUMA_RETURN(YUP) ;
1482 }
1483
1484 /*!
1485 SUMA's version of AFNI_gicor_setref
1486 */
SUMA_AFNI_gicor_setref(SUMA_SurfaceObject * SO,int node)1487 int SUMA_AFNI_gicor_setref( SUMA_SurfaceObject *SO, int node )
1488 {
1489 static char FuncName[]={"SUMA_AFNI_gicor_setref"};
1490 NI_element *nel=NULL;
1491 char buf[256]={"bise"}, boom[64]={""}, *cside=NULL;
1492 GICOR_setup *giset = SUMAg_CF->giset ;
1493 THD_fvec3 iv,jv; THD_ivec3 kv;
1494 int ijk=-1,ii=0;
1495 SUMA_SurfaceObject *SOv[2]={NULL, NULL};
1496 SUMA_Boolean LocalHead = NOPE;
1497
1498 SUMA_ENTRY;
1499
1500 if (!SO) {
1501 SUMA_S_Err("NULL input");
1502 SUMA_RETURN(-1) ;
1503 }
1504
1505 if (node < 0) { /* OK, return */
1506 SUMA_LHv("node %d\n", node);
1507 SUMA_RETURN(0) ;
1508 }
1509
1510 if( giset == NULL ||
1511 !giset->ready ){ /* should not happen */
1512 SUMA_LHv("giset=%p, giset->ready=%d\n",
1513 giset, giset ? giset->ready:-111);
1514 if( giset != NULL ) giset->ready = 0 ;
1515 SUMA_RETURN(-1) ;
1516 }
1517
1518 /* change node index to proper ijk */
1519 if (!SUMA_GICOR_Surfaces(giset, SOv)) {
1520 SUMA_S_Err("Failed to find surfaces for giset");
1521 SUMA_RETURN(-1);
1522 }
1523
1524 if (SUMA_isRelated_SO(SO, SOv[0],1)) {
1525 ijk = node;
1526 } else if (SUMA_isRelated_SO(SO, SOv[1],1)) {
1527 ijk = node+giset->nnode_domain[0];
1528 } else {
1529 SUMA_SLP_Warn("Cannot change node to ijk");
1530 SUMA_RETURN(-1);
1531 }
1532
1533 cside = "";
1534 if (SO->Side == SUMA_LEFT) cside = "L";
1535 else if (SO->Side == SUMA_RIGHT) cside = "R";
1536 sprintf(boom,"%d%s", node, cside);
1537 SUMA_LHv("Node %d --> ijk %d\n", node, ijk);
1538
1539 /* if socket has gone bad, we're done */
1540
1541 if( NI_stream_goodcheck(giset->ns,1) < 1 ){
1542 SUMA_S_Note("Socket socks, toggling off connection");
1543 SUMAg_CF->Connected_v[SUMA_GICORR_LINE]=NOPE;
1544 if( giset != NULL ) giset->ready = 0 ;
1545 SUMA_RETURN(-1) ;
1546 }
1547
1548 /* find where we are */
1549
1550 /* INFO_message("DEBUG: AFNI_gicor_setref called: ijk=%d",ijk) ; */
1551
1552 if( giset->ivec != NULL ){
1553 ii = bsearch_int( ijk , giset->nvec , giset->ivec ) ;
1554 if( ii < 0 ){
1555 WARNING_message("GrpInCorr set point not in mask from 3dGroupInCorr") ;
1556 SUMA_RETURN(0) ; /* don't return an error */
1557 }
1558 }
1559
1560 /* send ijk node index to 3dGroupInCorr */
1561 nel = NI_new_data_element( "SETREF_ijk" , 0 ) ;
1562 NI_set_attribute(nel, "boomerang_msg", boom);
1563
1564 sprintf( buf , "%d" , ijk ) ;
1565 NI_set_attribute( nel , "index" , buf ) ;
1566
1567 sprintf( buf , "%g" , giset->seedrad ) ;
1568 NI_set_attribute( nel , "seedrad" , buf ) ;
1569
1570 sprintf( buf , "%d" , giset->ttest_opcode ) ;
1571 NI_set_attribute( nel , "ttest_opcode" , buf ) ;
1572
1573 ii = NI_write_element( giset->ns , nel , NI_TEXT_MODE ) ;
1574 NI_free_element( nel ) ;
1575 if( ii <= 0 ){
1576 ERROR_message("3dGroupInCorr connection has failed?!") ;
1577 SUMA_RETURN(-1) ;
1578 }
1579
1580 SUMA_RETURN(0) ;
1581 }
1582
1583
SUMA_Show_GISET(GICOR_setup * giset,FILE * out,int verb)1584 void SUMA_Show_GISET(GICOR_setup *giset, FILE *out, int verb) {
1585 static char FuncName[]={"SUMA_Show_GISET"};
1586 char *s=NULL;
1587
1588 SUMA_ENTRY;
1589
1590 s = SUMA_GISET_Info(giset, verb);
1591
1592 if (!out) {
1593 out = SUMA_STDOUT;
1594 }
1595
1596 fprintf(out, "%s\n", s);
1597
1598 SUMA_free(s); s = NULL;
1599
1600 SUMA_RETURNe;
1601 }
1602
SUMA_GISET_Info(GICOR_setup * giset,int verb)1603 char *SUMA_GISET_Info(GICOR_setup *giset, int verb) {
1604 static char FuncName[]={"SUMA_GISET_Info"};
1605 char *s=NULL;
1606 SUMA_STRING *SS=NULL;
1607
1608 SUMA_ENTRY;
1609
1610 SS = SUMA_StringAppend_va(NULL, NULL);
1611
1612 if (giset) {
1613 SS = SUMA_StringAppend(SS, " GICORR-setup:\n");
1614 SS = SUMA_StringAppend_va(SS, " ready: %d\n"
1615 " ndset: %d %d, nvec: %d\n"
1616 " ttestopcode: %d, vmul: %d\n"
1617 " seedrad: %f\n"
1618 " ns: %p\n"
1619 " session: %p, dset: %p (%s)\n"
1620 " nds:%d, nvox: %d\n"
1621 " nivec: %d, ivec: %p\n"
1622 " sdset_ID: %s, %s\n"
1623 " nnode_domain: %d, %d\n"
1624 " nnode_mask: %d %d\n",
1625 giset->ready,
1626 giset->ndset_A , giset->ndset_B , giset->nvec,
1627 giset->ttest_opcode , giset->vmul, giset->seedrad,
1628 giset->ns, giset->session, giset->dset,
1629 giset->dset ? DSET_PREFIX(giset->dset):"NULL",
1630 giset->nds, giset->nvox, giset->nivec, giset->ivec,
1631 giset->sdset_ID[0] ? giset->sdset_ID[0]:"NULL",
1632 giset->sdset_ID[1] ? giset->sdset_ID[1]:"NULL",
1633 giset->nnode_domain[0], giset->nnode_domain[1],
1634 giset->nnode_mask[0], giset->nnode_mask[1]);
1635 } else {
1636 SS = SUMA_StringAppend_va(SS, " GICORR-setup: NULL\n");
1637 }
1638
1639 SS = SUMA_StringAppend_va(SS, NULL);
1640 s = SS->s; SUMA_free(SS); SS= NULL;
1641 SUMA_RETURN(s);
1642 }
1643