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
9 static int BuildMethod = SUMA_OFFSETS2_NO_REC;
10
11 /* local protos */
12 static int SUMA_nstat_mean(SUMA_DSET *din, int icol, float *fin_orig,
13 byte *nmask, byte strict_mask, int nnode,
14 SUMA_OFFSET_STRUCT *OffS_out, float rhood,
15 SUMA_DSET *dout, float * fout, int ndbg);
16 static int SUMA_nstat_num (SUMA_DSET *din, int icol, float *fin_orig,
17 byte *nmask, byte strict_mask, int nnode,
18 SUMA_OFFSET_STRUCT *OffS_out, float rhood,
19 SUMA_DSET *dout, float * fout, int ndbg);
20 static int SUMA_nstat_mode(SUMA_DSET *din, int icol, float *fin_orig,
21 byte *nmask, byte strict_mask, int nnode,
22 SUMA_OFFSET_STRUCT *OffS_out, float rhood,
23 SUMA_DSET *dout, float * fout, int ndbg);
24 static int SUMA_world_stats_node_dbg_f(SUMA_OFFSET_STRUCT * OffS_out,
25 int node, int nval, int nnodes, char * lblcp,
26 byte * nmask, byte strict_mask,
27 float rhood, float * fin_orig, float * fout);
28 static int float_list_comp_mode( float_list *f, float *mode, int *nvals,
29 int *index );
30
SUMA_SurfClust_Set_Method(int m)31 void SUMA_SurfClust_Set_Method(int m)
32 {
33 BuildMethod = m;
34 return;
35 }
SUMA_SurfClust_Get_Method(void)36 int SUMA_SurfClust_Get_Method(void)
37 { return(BuildMethod); }
38
SUMA_FreeClustDatum(void * data)39 void SUMA_FreeClustDatum (void * data)
40 {
41 static char FuncName[]={"SUMA_FreeClustDatum"};
42 SUMA_CLUST_DATUM *Clust = NULL;
43 SUMA_ENTRY;
44
45 if (!data) SUMA_RETURNe;
46 Clust = (SUMA_CLUST_DATUM *)data;
47 if (Clust->NodeList) SUMA_free(Clust->NodeList);
48 if (Clust->ValueList) SUMA_free(Clust->ValueList);
49 SUMA_free(Clust);
50
51 SUMA_RETURNe;
52 }
53
54 /*!
55 \brief Calculate area of each node as one third of the sum of
56 the areas of incident triangles.
57 mask: if not null, then calculate areas for nodes in mask only
58 If you change this function make sure changes are also done
59 on RickR's compute_node_areas since the two functions use
60 the same principle.
61 */
SUMA_CalculateNodeAreas(SUMA_SurfaceObject * SO,byte * mask)62 float *SUMA_CalculateNodeAreas(SUMA_SurfaceObject *SO, byte *mask)
63 {
64 static char FuncName[]={"SUMA_CalculateNodeAreas"};
65 float *NodeAreas=NULL;
66 int *flist = NULL, i, c;
67
68 SUMA_ENTRY;
69
70 if (!SO) { SUMA_RETURN(NodeAreas); }
71 if (!SO->PolyArea || !SO->MF) {
72 if (!SUMA_SurfaceMetrics_eng(SO, "PolyArea|MemberFace", NULL, 0,
73 SUMAg_CF->DsetList)) {
74 fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_SurfaceMetrics.\n",
75 FuncName);
76 SUMA_RETURN(NodeAreas);
77 }
78 }
79
80 NodeAreas = (float *)SUMA_malloc(SO->N_Node*sizeof(float));
81 if (!NodeAreas) {
82 SUMA_SL_Crit ("Failed to allocate for NodeAreas");
83 SUMA_RETURN(NodeAreas);
84 }
85
86 for (i=0; i<SO->N_Node; ++i) {
87 NodeAreas[i] = 0.0;
88 if (!mask || mask[i]) {
89 flist = SO->MF->NodeMemberOfFaceSet[i];
90 for (c = 0; c < SO->MF->N_Memb[i]; c++) {
91 NodeAreas[i] += SO->PolyArea[flist[c]];
92 }
93 NodeAreas[i] /= 3.0;
94 }
95 }
96
97 SUMA_RETURN(NodeAreas);
98 }
99
100
101 /*!
102 \brief builds a cluster starting from some node
103
104 \param dothisnode (int) start building from this node
105 \param AddToThisClust (SUMA_CLUST_DATUM *)pointer to cluster to add to.
106 Send NULL when function is first called.
107 This is a non NULL when the function recurses.
108 \param ToBeAssigned (float *) if ToBeAssigned[i] then node i
109 (index into SO's nodelist) is considered in the clustering
110 and the value at that node is ToBeAssigned[i].
111 Vector is SO->N_Node elements long.
112 Gets modified during recursion.
113 \param N_ToBeAssigned (int *) pointer to number of value values
114 in ToBeAssigned.
115 Changes with recursive calls.
116 \param NodeArea(float *) Vector containing area of each node of surface
117 \param SO (SUMA_SurfaceObject *) the usual deal
118 \param Opt (SUMA_SURFCLUST_OPTIONS *) options for cluster building.
119
120 \sa recursive calls SUCK a lot of memory and are a nightmare to track with
121 function I/O. They are also slow and not necessary!
122
123 \sa SUMA_Build_Cluster_From_Node_NoRec
124 */
SUMA_Build_Cluster_From_Node(int dothisnode,SUMA_CLUST_DATUM * AddToThisClust,float * ToBeAssigned,int * N_TobeAssigned,float * NodeArea,SUMA_SurfaceObject * SO,SUMA_SURFCLUST_OPTIONS * Opt)125 SUMA_CLUST_DATUM * SUMA_Build_Cluster_From_Node(
126 int dothisnode, SUMA_CLUST_DATUM *AddToThisClust,
127 float *ToBeAssigned, int *N_TobeAssigned, float *NodeArea,
128 SUMA_SurfaceObject *SO, SUMA_SURFCLUST_OPTIONS *Opt)
129 {
130 static char FuncName[]={"SUMA_Build_Cluster_From_Node"};
131 SUMA_CLUST_DATUM *Clust = NULL;
132 static int ncall;
133 int il, jl, neighb, nnn, kkk;
134 SUMA_GET_OFFSET_STRUCT *OffS = NULL;
135 DList *offlist = NULL;
136 DListElmt *elm = NULL;
137 SUMA_OFFSET_LL_DATUM *dat=NULL;
138 int NewClust = 0;
139 SUMA_Boolean LocalHead = NOPE;
140
141 SUMA_ENTRY;
142 ++ncall;
143 if (dothisnode < 0) {
144 SUMA_SL_Err("Unexpected negative index.");
145 SUMA_RETURN(NULL);
146 }
147 if (!AddToThisClust) {
148 Clust = (SUMA_CLUST_DATUM *)SUMA_malloc(sizeof(SUMA_CLUST_DATUM));
149 Clust->N_Node = 0; Clust->totalarea = 0.0; /* Clust->rank = -1; */
150 Clust->totalvalue = 0.0; Clust->totalabsvalue = 0.0;
151 Clust->minvalue = ToBeAssigned[dothisnode]; Clust->minnode = dothisnode;
152 Clust->maxvalue = ToBeAssigned[dothisnode]; Clust->maxnode = dothisnode;
153 Clust->minabsvalue = SUMA_ABS(Clust->minvalue);
154 Clust->minabsnode = Clust->minnode;
155 Clust->maxabsvalue = SUMA_ABS(Clust->maxvalue);
156 Clust->maxabsnode = Clust->maxnode;
157 Clust->varvalue = 0.0; Clust->centralnode = 0;
158 Clust->weightedcentralnode = 0;
159 Clust->NodeList = (int *)SUMA_malloc((*N_TobeAssigned) * sizeof(int));
160 Clust->ValueList = (float *)SUMA_malloc((*N_TobeAssigned) * sizeof(float));
161 memset(Clust->com, 0, 3*sizeof(float));
162 memset(Clust->cen, 0, 3*sizeof(float));
163 if (!Clust->NodeList || !Clust->ValueList) {
164 SUMA_SL_Crit("Failed to allocate for NodeList or ValueList");
165 SUMA_free(Clust); Clust = NULL;
166 SUMA_RETURN(NULL);
167 }
168 NewClust = 1;
169 if (LocalHead)
170 fprintf (SUMA_STDERR,"%s: New Cluster %p, with node %d\n",
171 FuncName, Clust, dothisnode);
172 } else {
173 NewClust = 0;
174 Clust = AddToThisClust;
175 if (LocalHead)
176 fprintf (SUMA_STDERR,"%s: Reusing Cluster %p, with node %d\n",
177 FuncName, Clust, dothisnode);
178 }
179
180 /* Add node to cluster */
181 if (LocalHead)
182 fprintf(SUMA_STDERR,"%s: Adding node %d to cluster %p of %d nodes\n",
183 FuncName, dothisnode, Clust, Clust->N_Node);
184 Clust->NodeList[Clust->N_Node] = dothisnode;
185 Clust->totalarea += NodeArea[dothisnode];
186 Clust->totalvalue += ToBeAssigned[dothisnode];
187 Clust->totalabsvalue += (float)fabs((float)ToBeAssigned[dothisnode]);
188 if (ToBeAssigned[dothisnode] < Clust->minvalue) {
189 Clust->minvalue = ToBeAssigned[dothisnode];
190 Clust->minnode = dothisnode;
191 }
192 if (ToBeAssigned[dothisnode] > Clust->maxvalue) {
193 Clust->maxvalue = ToBeAssigned[dothisnode];
194 Clust->maxnode = dothisnode;
195 }
196 if (SUMA_ABS(ToBeAssigned[dothisnode]) < Clust->minabsvalue) {
197 Clust->minabsvalue = SUMA_ABS(ToBeAssigned[dothisnode]);
198 Clust->minabsnode = dothisnode;
199 }
200 if (SUMA_ABS(ToBeAssigned[dothisnode]) > Clust->maxabsvalue) {
201 Clust->maxabsvalue = SUMA_ABS(ToBeAssigned[dothisnode]);
202 Clust->maxabsnode = dothisnode;
203 }
204 Clust->ValueList[Clust->N_Node] = ToBeAssigned[dothisnode];
205 nnn = SO->NodeDim*dothisnode;
206 for (kkk = 0; kkk<SO->NodeDim; ++kkk) {
207 Clust->com[kkk] += ToBeAssigned[dothisnode] * SO->NodeList[nnn+kkk];
208 Clust->cen[kkk] += SO->NodeList[nnn+kkk];
209 }
210 ++Clust->N_Node;
211
212 /* mark it as assigned, an reduce the number of nodes left to assign*/
213 ToBeAssigned[dothisnode] = 0; --(*N_TobeAssigned);
214
215 if (BuildMethod == SUMA_OFFSETS2) {
216 /* Tres bad memory utilization due to recursive calls */
217 if (*N_TobeAssigned) {
218 /* look in its vicinity - bad memory usage due to recursive calls*/
219 OffS = SUMA_Initialize_getoffsets (SO->N_Node);
220 SUMA_getoffsets2 (dothisnode, SO, Opt->DistLim, OffS, NULL, 0);
221 #if 0
222 if (NewClust) {
223 FILE *fid=NULL;
224 char *s=NULL, tmp[50];
225 fid = fopen("offsets2.1D", "w");
226 if (!fid) {
227 SUMA_SL_Err("Could not open file for writing.\n"
228 "Check file permissions, disk space.\n");
229 } else {
230 s = SUMA_ShowOffset_Info(OffS, 0);
231 if (s) { fprintf(fid,"%s\n", s); SUMA_free(s); s = NULL;}
232 fclose(fid);
233 }
234 }
235 #endif
236 /* search to see if any are to be assigned */
237 if (Opt->DistLim >= 0.0) {
238 for (il=1; il<OffS->N_layers; ++il) {
239 /* starting at layer 1,
240 layer 0 is the node itself */
241 for (jl=0; jl<OffS->layers[il].N_NodesInLayer; ++jl) {
242 neighb = OffS->layers[il].NodesInLayer[jl];
243 if ( ToBeAssigned[neighb] &&
244 OffS->OffVect[neighb] <= Opt->DistLim) {
245 /* take that node into the cluster */
246 SUMA_Build_Cluster_From_Node(
247 neighb, Clust, ToBeAssigned,
248 N_TobeAssigned, NodeArea, SO, Opt);
249 }
250 }
251 }
252 } else { /* accept nodes connect by -((int)Opt->DistLim) edges or less*/
253 for (il=1; il <= -((int)Opt->DistLim); ++il) {
254 for (jl=0; jl<OffS->layers[il].N_NodesInLayer; ++jl) {
255 neighb = OffS->layers[il].NodesInLayer[jl];
256 if ( ToBeAssigned[neighb] ) {
257 /* take that node into the cluster */
258 SUMA_Build_Cluster_From_Node(
259 neighb, Clust, ToBeAssigned,
260 N_TobeAssigned, NodeArea, SO, Opt);
261 }
262 }
263 }
264 }
265 /* free this OffS structure (Note you can't recycle the same
266 structure because you are using many OffS at one because
267 of recursive calls */
268 if (OffS) SUMA_Free_getoffsets(OffS); OffS = NULL;
269 }
270 } else if (BuildMethod == SUMA_OFFSETS_LL) {
271 if (*N_TobeAssigned) {
272 /* look in its vicinity */
273 if (!(offlist = SUMA_getoffsets_ll (dothisnode, SO, Opt->DistLim,
274 NULL, 0))) {
275 SUMA_SL_Err("Failed to get offsets.\nNo cleanup done.");
276 SUMA_RETURN(NULL);
277 }
278 #if 0
279 if (NewClust) {
280 FILE *fid=NULL;
281 char *s=NULL, tmp[50];
282 fid = fopen("offsets_ll.1D", "w");
283 if (!fid) {
284 SUMA_SL_Err("Could not open file for writing.\n"
285 "Check file permissions, disk space.\n");
286 } else {
287 s = SUMA_ShowOffset_ll_Info(offlist, 0);
288 if (s) { fprintf(fid,"%s\n", s); SUMA_free(s); s = NULL;}
289 fclose(fid);
290 }
291 }
292 #endif
293
294 /* search to see if any are to be assigned, start at layer 1*/
295 elm = dlist_head(offlist);
296 dat = (SUMA_OFFSET_LL_DATUM *)elm->data;
297 if (dat->layer != 0) {
298 SUMA_SL_Err("Unexpected non zero layer for first element.");
299 SUMA_RETURN(NULL);
300 }
301 do {
302 elm = elm->next;
303 dat = (SUMA_OFFSET_LL_DATUM *)elm->data;
304 neighb = dat->ni;
305 if (ToBeAssigned[neighb]) {
306 if (dat->off <= Opt->DistLim) {
307 /* take that node into the cluster */
308 SUMA_Build_Cluster_From_Node( neighb, Clust, ToBeAssigned,
309 N_TobeAssigned, NodeArea, SO,
310 Opt);
311 } else if (dat->layer <= -(int)Opt->DistLim) {
312 /* take that node into the cluster */
313 SUMA_Build_Cluster_From_Node( neighb, Clust, ToBeAssigned,
314 N_TobeAssigned, NodeArea, SO,
315 Opt);
316 }
317 }
318 } while (elm != dlist_tail(offlist));
319 dlist_destroy(offlist); SUMA_free(offlist); offlist = NULL;
320 }
321 }
322
323
324 SUMA_RETURN(Clust);
325 }
326
327 /*!
328 A macro for SUMA_Build_Cluster_From_Node_NoRec
329 */
330 #define SUMA_ADD_NODE_TO_CLUST(dothisnode, Clust, NodeArea, ToBeAssigned, SO){ \
331 int nnn, kkk; \
332 if (LocalHead) fprintf(SUMA_STDERR,"%s: Adding node %d to cluster %p of %d nodes\nClustMax is %f, incoming value is %f\n", FuncName, dothisnode, Clust, Clust->N_Node, Clust->maxvalue, ToBeAssigned[dothisnode]); \
333 Clust->NodeList[Clust->N_Node] = dothisnode; \
334 Clust->totalarea += NodeArea[dothisnode]; \
335 Clust->totalvalue += ToBeAssigned[dothisnode]; \
336 Clust->totalabsvalue += (float)fabs((float)ToBeAssigned[dothisnode]); \
337 if (ToBeAssigned[dothisnode] < Clust->minvalue) { Clust->minvalue = ToBeAssigned[dothisnode]; Clust->minnode = dothisnode; } \
338 if (ToBeAssigned[dothisnode] > Clust->maxvalue) { Clust->maxvalue = ToBeAssigned[dothisnode]; Clust->maxnode = dothisnode; } \
339 if (SUMA_ABS(ToBeAssigned[dothisnode]) < Clust->minabsvalue) { Clust->minabsvalue = SUMA_ABS(ToBeAssigned[dothisnode]); Clust->minabsnode = dothisnode; }\
340 if (SUMA_ABS(ToBeAssigned[dothisnode]) > Clust->maxabsvalue) { Clust->maxabsvalue = SUMA_ABS(ToBeAssigned[dothisnode]); Clust->maxabsnode = dothisnode; }\
341 Clust->ValueList[Clust->N_Node] = ToBeAssigned[dothisnode]; \
342 nnn = SO->NodeDim*dothisnode; \
343 for (kkk = 0; kkk<SO->NodeDim; ++kkk) {\
344 Clust->com[kkk] += ToBeAssigned[dothisnode] * SO->NodeList[nnn+kkk]; \
345 Clust->cen[kkk] += SO->NodeList[nnn+kkk]; \
346 } \
347 ++Clust->N_Node; \
348 }
349
350 /*!
351 \brief builds a cluster starting from some node
352
353 \param dothisnode (int) start building from this node
354 \param ToBeAssigned (float *) if ToBeAssigned[i] then node i (index into
355 SO's nodelist) is considered in the clustering
356 and the value at that node is ToBeAssigned[i].
357 Vector is SO->N_Node elements long. Gets
358 modified in function.
359 \param N_ToBeAssigned (int *) pointer to number of value values in
360 ToBeAssigned. Gets modified in function.
361 \param NodeArea(float *) Vector containing area of each node of surface
362 \param SO (SUMA_SurfaceObject *) the usual deal
363 \param Opt (SUMA_SURFCLUST_OPTIONS *) options for cluster building.
364
365 \sa Based on recursive SUMA_Build_Cluster_From_Node
366 */
SUMA_Build_Cluster_From_Node_NoRec(int dothisnode,float * ToBeAssigned,int * N_TobeAssigned,float * NodeArea,SUMA_SurfaceObject * SO,SUMA_SURFCLUST_OPTIONS * Opt)367 SUMA_CLUST_DATUM * SUMA_Build_Cluster_From_Node_NoRec ( int dothisnode,
368 float *ToBeAssigned, int *N_TobeAssigned, float *NodeArea,
369 SUMA_SurfaceObject *SO, SUMA_SURFCLUST_OPTIONS *Opt )
370 {
371 static char FuncName[]={"SUMA_Build_Cluster_From_Node_NoRec"};
372 SUMA_CLUST_DATUM *Clust = NULL;
373 static int ncall;
374 static int N_Orig = -1;
375 int il, jl, neighb;
376 void *dtmp=NULL;
377 SUMA_GET_OFFSET_STRUCT *OffS = NULL;
378 DList *offlist = NULL, *candlist=NULL;
379 DListElmt *elm = NULL, *dothiselm=NULL;
380 SUMA_OFFSET_LL_DATUM *dat=NULL;
381 int NewClust = 0, NeedNewLine = 0;
382 byte *visited=NULL;
383 SUMA_Boolean LocalHead = NOPE;
384
385 SUMA_ENTRY;
386
387 ++ncall;
388
389 /* a trick to know when to initialize */
390 if (Opt->update < 0) { N_Orig = *N_TobeAssigned; Opt->update = -Opt->update; }
391
392 if (dothisnode < 0) {
393 SUMA_SL_Err("Unexpected negative index.");
394 SUMA_RETURN(NULL);
395 }
396
397 OffS = SUMA_Initialize_getoffsets (SO->N_Node);
398 Clust = (SUMA_CLUST_DATUM *)SUMA_malloc(sizeof(SUMA_CLUST_DATUM));
399 Clust->N_Node = 0; Clust->totalarea = 0.0; /* Clust->rank = -1; */
400 Clust->totalvalue = 0.0; Clust->totalabsvalue = 0.0;
401 Clust->minvalue = ToBeAssigned[dothisnode]; Clust->minnode = dothisnode;
402 Clust->maxvalue = ToBeAssigned[dothisnode]; Clust->maxnode = dothisnode;
403 Clust->minabsvalue = SUMA_ABS(Clust->minvalue);
404 Clust->minabsnode = Clust->minnode;
405 Clust->maxabsvalue = SUMA_ABS(Clust->maxvalue);
406 Clust->maxabsnode = Clust->maxnode;
407 Clust->varvalue = 0.0; Clust->centralnode = 0;
408 Clust->weightedcentralnode = 0;
409 Clust->NodeList = (int *)SUMA_malloc((*N_TobeAssigned) * sizeof(int));
410 Clust->ValueList = (float *)SUMA_malloc((*N_TobeAssigned) * sizeof(float));
411 memset(Clust->com, 0, 3*sizeof(float));
412 memset(Clust->cen, 0, 3*sizeof(float));
413 if (!Clust->NodeList || !Clust->ValueList || !OffS) {
414 SUMA_SL_Crit("Failed to allocate for NodeList or ValueList");
415 SUMA_free(Clust); Clust = NULL;
416 SUMA_RETURN(NULL);
417 }
418 candlist = (DList*)SUMA_malloc(sizeof(DList));
419 visited = (byte *)SUMA_calloc(SO->N_Node, sizeof(byte));
420 if (!visited || !candlist) {
421 SUMA_SL_Crit("Failed to allocate for visited or candlist");
422 SUMA_free(Clust); Clust = NULL;
423 SUMA_RETURN(NULL);
424 }
425 dlist_init(candlist, NULL);
426 /* Add node to cluster */
427 SUMA_ADD_NODE_TO_CLUST(dothisnode, Clust, NodeArea, ToBeAssigned, SO);
428 /* mark it as assigned, an reduce the number of nodes left to assign*/
429 ToBeAssigned[dothisnode] = 0; --(*N_TobeAssigned);
430 visited[dothisnode] = YUP;
431 dlist_ins_next(candlist, dlist_tail(candlist), (VOID_CAST)dothisnode);
432 while (*N_TobeAssigned && dlist_size(candlist)) {
433 /* look in its vicinity */
434 dothiselm = dlist_head(candlist);
435 dothisnode = (INT_CAST) dothiselm->data;
436 SUMA_getoffsets2 (dothisnode, SO, Opt->DistLim, OffS, NULL, 0);
437 /* remove node from candidate list */
438 dlist_remove(candlist, dothiselm, (void*)&dtmp);
439 /* search to see if any are to be assigned */
440 if (Opt->DistLim >= 0.0) {
441 for (il=1; il<OffS->N_layers; ++il) {
442 /* starting at layer 1, layer 0 is the node itself */
443 for (jl=0; jl<OffS->layers[il].N_NodesInLayer; ++jl) {
444 neighb = OffS->layers[il].NodesInLayer[jl];
445 if ( ToBeAssigned[neighb] &&
446 OffS->OffVect[neighb] <= Opt->DistLim) {
447 /* take that node into the cluster */
448 SUMA_ADD_NODE_TO_CLUST( neighb, Clust,
449 NodeArea, ToBeAssigned, SO);
450 /* mark it as assigned, an reduce the number
451 of nodes left to assign*/
452 ToBeAssigned[neighb] = 0; --(*N_TobeAssigned);
453 if (Opt->update) {
454 if (N_Orig - *N_TobeAssigned >= Opt->update) {
455 if (LocalHead)
456 fprintf( SUMA_STDERR,
457 "%s: tick (%d nodes processed)\n",
458 FuncName, N_Orig - *N_TobeAssigned);
459 else fprintf(SUMA_STDERR,".");
460 NeedNewLine = 1;
461
462 N_Orig = *N_TobeAssigned;
463 }
464 }
465 /* mark it as a candidate if it has not been visited as
466 a candidate before */
467 if (!visited[neighb]) {
468 dlist_ins_next(candlist, dlist_tail(candlist),
469 (VOID_CAST)neighb);
470 visited[neighb] = YUP;
471 }
472 }
473 }
474 }
475 } else {/* accept nodes connected by up to -((int)Opt->DistLim) edges */
476 for (il=1; il <= -((int)Opt->DistLim); ++il) {
477 for (jl=0; jl<OffS->layers[il].N_NodesInLayer; ++jl) {
478 neighb = OffS->layers[il].NodesInLayer[jl];
479 if ( ToBeAssigned[neighb] ) {
480 /* take that node into the cluster */
481 SUMA_ADD_NODE_TO_CLUST( neighb, Clust,
482 NodeArea, ToBeAssigned, SO);
483 /* mark it as assigned, an reduce the number
484 of nodes left to assign*/
485 ToBeAssigned[neighb] = 0; --(*N_TobeAssigned);
486 if (Opt->update) {
487 if (N_Orig - *N_TobeAssigned >= Opt->update) {
488 if (LocalHead)
489 fprintf( SUMA_STDERR,
490 "%s: tick (%d nodes processed)\n",
491 FuncName, N_Orig - *N_TobeAssigned);
492 else fprintf(SUMA_STDERR,".");
493 NeedNewLine = 1;
494 N_Orig = *N_TobeAssigned;
495 }
496 }
497 /* mark it as a candidate if it has not been visited as
498 a candidate before */
499 if (!visited[neighb]) {
500 dlist_ins_next(candlist, dlist_tail(candlist),
501 (VOID_CAST)neighb);
502 visited[neighb] = YUP;
503 }
504 }
505 }
506 }
507 }
508 /* recycle */
509 SUMA_Recycle_getoffsets (OffS);
510 }
511 /* free this OffS structure */
512 if (OffS) SUMA_Free_getoffsets(OffS); OffS = NULL;
513 if (NeedNewLine) fprintf(SUMA_STDERR,"\n");
514 /* destroy the list */
515 if (candlist) {
516 dlist_destroy(candlist); SUMA_free(candlist); candlist = NULL;
517 }
518 SUMA_RETURN(Clust);
519 }
520
SUMA_ClustCommandLineFromOpt(char * pname,SUMA_SurfaceObject * SO,SUMA_SURFCLUST_OPTIONS * Opt,char * filler)521 char *SUMA_ClustCommandLineFromOpt(char *pname, SUMA_SurfaceObject *SO,
522 SUMA_SURFCLUST_OPTIONS *Opt,
523 char *filler)
524 {
525 static char FuncName[]={"SUMA_ClustCommandLineFromOpt"};
526 SUMA_STRING *SS=NULL;
527 char *s = NULL;
528
529 SUMA_ENTRY;
530
531 if (!Opt) SUMA_RETURN(s);
532
533 SS = SUMA_StringAppend(NULL, NULL);
534
535 if (!pname) pname = "SurfClust";
536
537 if (SO) {
538 if (!(s = SUMA_SurfaceFileName(SO,1))) {
539 s = SUMA_copy_string("Qui_Etes_Vous.gii");
540 }
541 } else {
542 s = SUMA_copy_string("SOunknown");
543 }
544 SS = SUMA_StringAppend_va(SS,"%s -i %s -input %s %d "
545 "-rmm %f "
546 , pname, s, Opt->in_name,
547 Opt->labelcol, Opt->DistLim );
548 if (s) SUMA_free(s); s= NULL;
549 if (Opt->tind>=0) {
550 if (Opt->DoThreshold == SUMA_NO_THRESH) {
551 } else if (Opt->DoThreshold == SUMA_LESS_THAN) {
552 SS = SUMA_StringAppend_va(SS,"-thresh_col %d -thresh %f ",
553 Opt->tind, Opt->ThreshR[0]);
554 } else if (Opt->DoThreshold == SUMA_ABS_LESS_THAN) {
555 SS = SUMA_StringAppend_va(SS,"-thresh_col %d -athresh %f ",
556 Opt->tind, Opt->ThreshR[0]);
557 } else if (Opt->DoThreshold == SUMA_THRESH_OUTSIDE_RANGE) {
558 SS = SUMA_StringAppend_va(SS,"-thresh_col %d -ex_range %f %f",
559 Opt->tind, Opt->ThreshR[0], Opt->ThreshR[1]);
560 } else if (Opt->DoThreshold == SUMA_THRESH_INSIDE_RANGE) {
561 SS = SUMA_StringAppend_va(SS,"-thresh_col %d -in_range %f %f",
562 Opt->tind, Opt->ThreshR[0], Opt->ThreshR[1]);
563 } else {
564 SS = SUMA_StringAppend(SS,"NO_COMPARABLE_THRESHOLD ");
565 }
566 }
567 if (Opt->AreaLim < 0) {
568 SS = SUMA_StringAppend_va(SS,"-n %d ", (int)-Opt->AreaLim);
569 } else if (Opt->AreaLim > 0){
570 SS = SUMA_StringAppend_va(SS,"-amm2 %f ", Opt->AreaLim);
571 }
572 switch (Opt->SortMode) {
573 case SUMA_SORT_CLUST_NO_SORT:
574 SS = SUMA_StringAppend(SS,"-sort_none ");
575 break;
576 case SUMA_SORT_CLUST_BY_NUMBER_NODES:
577 SS = SUMA_StringAppend(SS,"-sort_n_nodes ");
578 break;
579 case SUMA_SORT_CLUST_BY_AREA:
580 SS = SUMA_StringAppend(SS,"-sort_area ");
581 break;
582 default:
583 SS = SUMA_StringAppend(SS,"-sort_none ");
584 break;
585 }
586 if (Opt->DoCentrality) {
587 SS = SUMA_StringAppend(SS,"-cent ");
588 } else {
589 SS = SUMA_StringAppend(SS,"-no_cent ");
590 }
591 if (filler) SS = SUMA_StringAppend(SS,filler);
592
593 SUMA_SS2S(SS,s);
594
595 SUMA_RETURN(s);
596 }
597
SUMA_AvgGradient(SUMA_SurfaceObject * SO,float ** FirstNeighbDist,float * nv,byte * mask,byte mask_zeros,SUMA_GRAD_SCALE_OPTS normopt)598 float *SUMA_AvgGradient(SUMA_SurfaceObject *SO, float **FirstNeighbDist,
599 float *nv, byte *mask,
600 byte mask_zeros, SUMA_GRAD_SCALE_OPTS normopt)
601 {
602 static char FuncName[]={"SUMA_AvgGradient"};
603 int i=0, k=0, ki=0, ng=0, nmsk=0;
604 float *gr=NULL, avg=0.0;
605 double gavg=0.0;
606 SUMA_Boolean LocalHead=NOPE;
607
608 SUMA_ENTRY;
609
610 if (!SO || !SO->FN || !nv || !FirstNeighbDist) {
611 SUMA_S_Err("NULL input");
612 SUMA_RETURN(gr);
613 }
614
615 if (mask_zeros) {
616 if (!mask) {
617 mask = (byte *)SUMA_malloc(SO->N_Node*sizeof(byte));
618 memset(mask,1,SO->N_Node*sizeof(byte));
619 }
620 for (i=0; i<SO->N_Node; ++i) if (!nv[i]) mask[i]=0;
621 }
622
623 if (!(gr = (float *)SUMA_calloc(SO->N_Node, sizeof(float)))) {
624 SUMA_S_Crit("Failed to allocate");
625 SUMA_RETURN(NULL);
626 }
627 gavg = 0.0; nmsk = 0;
628 for (i=0; i<SO->N_Node; ++i) {
629 if (SO->FN->NodeId[i] != i) {
630 SUMA_S_Errv("Unexpected node index mismatch (%d,%d)."
631 "Please report case to author.",
632 SO->FN->NodeId[i], i);
633 SUMA_free(gr); gr=NULL; SUMA_RETURN(NULL);
634 }
635 gr[i]=0.0;
636 if (IN_MASK(mask, i)) {
637 ng = 0; ++nmsk;
638 avg = nv[i]; gavg += nv[i];
639 for (k=0; k<SO->FN->N_Neighb[i]; ++k) {
640 ki = SO->FN->FirstNeighb[i][k];
641 if (IN_MASK(mask,ki)) {
642 gr[i] += (nv[i]-nv[ki])/FirstNeighbDist[i][k]; ++ng;
643 avg += nv[ki];
644 }
645 }
646 if (ng) {
647 gr[i] /= (float)ng; /* average gradient */
648 avg = avg / (ng+1.0); /* local intensity average */
649 switch(normopt) {
650 case SUMA_GMEAN_GRAD_SCALE:
651 case SUMA_NO_GRAD_SCALE:
652 break;
653 case SUMA_MEAN_GRAD_SCALE:
654 gr[i] = 100.0*gr[i]/avg;
655 break;
656 default:
657 ERROR_message("Bad normalization option");
658 break;
659 }
660 }
661 }
662 }
663 if (normopt == SUMA_GMEAN_GRAD_SCALE && nmsk > 0 && gavg != 0.0) {
664 gavg = gavg/(float)nmsk/100.0;
665 for (i=0; i<SO->N_Node; ++i) {
666 if (IN_MASK(mask,ki)) {
667 gr[i] = gr[i] / gavg;
668 }
669 }
670 }
671
672 SUMA_RETURN(gr);
673 }
674
SUMA_DsetAvgGradient(SUMA_SurfaceObject * SO,float ** FirstNeighbDist,SUMA_DSET * din,byte * maskp,byte mask_by_zeros,SUMA_GRAD_SCALE_OPTS normopt)675 SUMA_DSET * SUMA_DsetAvgGradient(
676 SUMA_SurfaceObject *SO, float **FirstNeighbDist, SUMA_DSET *din,
677 byte *maskp, byte mask_by_zeros, SUMA_GRAD_SCALE_OPTS normopt)
678 {
679 static char FuncName[]={"SUMA_DsetAvgGradient"};
680 SUMA_DSET *dout = NULL;
681 int *icols = NULL, N_icols = -1, *ind = NULL, n_incopy=-1, masked_only=0;
682 int k = -1, N_mask=-1;
683 void *ncoli=NULL;
684 char *lblcp=NULL, *s=NULL;
685 float *fin_orig=NULL, *fout = NULL, fp = -1.0, **NeighbDist=NULL;
686 byte *mask=NULL;
687 SUMA_Boolean LocalHead = NOPE;
688
689 SUMA_ENTRY;
690
691
692 if (!SO || !din) {
693 SUMA_S_Err("Bad input");
694 SUMA_RETURN(NULL);
695 }
696 if (!SO->FN) {
697 if (!SUMA_SurfaceMetrics_eng(SO, "EdgeList|PolyArea|MemberFace",NULL,
698 0, SUMAg_CF->DsetList)) {
699 fprintf (SUMA_STDERR,
700 "Error %s: Failed in SUMA_SurfaceMetrics.\n", FuncName);
701 }
702 }
703
704 if (maskp) mask=maskp;
705
706 if (!FirstNeighbDist) {
707 NeighbDist = SUMA_CalcNeighbDist (SO);
708 } else {
709 NeighbDist = FirstNeighbDist;
710 }
711 if (!NeighbDist) {
712 SUMA_S_Err("NULL dists");
713 SUMA_RETURN(NULL);
714 }
715
716 /* what columns can we process ?*/
717 icols = SUMA_FindNumericDataDsetCols(din, &N_icols);
718 if (N_icols <= 0) {
719 SUMA_SL_Err("No approriate data columns in dset");
720 SUMA_RETURN(NOPE);
721 }
722 SUMA_LHv("Have %d numeric columns of input.\n", N_icols);
723 if (!(ind = SDSET_NODE_INDEX_COL(din))) {
724 SUMA_S_Note("Trying to populate the node index element");
725 if (!SUMA_PopulateDsetNodeIndexNel(din, 0)) {
726 SUMA_S_Err("Failed to populate NodeIndex Nel");
727 SUMA_RETURN(NULL);
728 }
729 }
730 /* Create a dset, at least as big as din*/
731 if ((ind = SDSET_NODE_INDEX_COL(din))) {
732 if (!masked_only) {
733 /* preserve all rows */
734 ncoli =
735 SUMA_Copy_Part_Column(ind,
736 NI_rowtype_find_code(SUMA_ColType2TypeCast(SUMA_NODE_INDEX)),
737 SDSET_VECLEN(din), NULL, masked_only, &n_incopy);
738 } else {
739 ncoli =
740 SUMA_Copy_Part_Column(ind,
741 NI_rowtype_find_code(SUMA_ColType2TypeCast(SUMA_NODE_INDEX)),
742 SDSET_VECLEN(din), mask, masked_only, &n_incopy);
743 }
744 if (!ncoli) {
745 SUMA_SL_Err("No index data got copied.");
746 SUMA_RETURN(NULL);
747 }
748 dout = SUMA_CreateDsetPointer("AvgGradient", SUMA_NODE_BUCKET, NULL,
749 SDSET_IDMDOM(din), n_incopy);
750 if (!SUMA_AddDsetNelCol (dout, NI_get_attribute(din->inel,"COLMS_LABS"), SUMA_NODE_INDEX, ncoli, NULL ,1)) {
751 SUMA_SL_Crit("Failed in SUMA_AddDsetNelCol");
752 SUMA_FreeDset((void*)dout); dout = NULL;
753 SUMA_RETURN(NULL);
754 }
755 if (ncoli) SUMA_free(ncoli); ncoli = NULL;
756 } else {
757 SUMA_S_Err( "Do not have node indices in input dset!\n"
758 " and could not create one.");
759 SUMA_RETURN(NULL);
760 }
761
762 /* Now, for each code, do the dance */
763 for (k=0; k < N_icols; ++k) {
764 /* get a float copy of the data column */
765 if (!(fin_orig = SUMA_DsetCol2FloatFullSortedColumn(din, icols[k], &mask,
766 0.0, SO->N_Node,
767 &N_mask, k==0?YUP:NOPE))){
768 SUMA_S_Err("Failed to extract");
769 SUMA_FreeDset(dout); dout=NULL;
770 if (!maskp) SUMA_free(mask); mask=NULL;
771 SUMA_RETURN(dout);
772 }
773
774 if (LocalHead) {
775 s = SUMA_ShowMeSome(fin_orig, SUMA_float, SO->N_Node, 10, NULL);
776 SUMA_LHv("fin_orig:\n%s\n", s); SUMA_free(s);
777 }
778
779
780 /* Now I have the data column, nice and solid , do the stats */
781 fout = SUMA_AvgGradient(SO, NeighbDist, fin_orig, mask,
782 mask_by_zeros, normopt);
783 if (!fout) {
784 SUMA_SL_Crit("Failed to compute gradient fout!");
785 SUMA_RETURN(NULL);
786 }
787
788 lblcp = SUMA_DsetColLabelCopy(din, icols[k], 1);
789 lblcp = SUMA_append_replace_string("D[", lblcp, "", 2);
790 lblcp = SUMA_append_replace_string(lblcp, "]", "", 1);
791 if (!SUMA_AddDsetNelCol (dout, lblcp, SUMA_NODE_FLOAT,
792 NULL, NULL ,1)) {
793 SUMA_S_Crit("Failed to add dset column");
794 SUMA_RETURN(NULL);
795 }
796 SUMA_free(lblcp); lblcp=NULL;
797 SUMA_LHv("Sticking column %d in dset (fout[0]=%f)\n", k, fout[0]);
798 if (!SUMA_Vec2DsetCol (dout, k, (void *)fout, SUMA_float,
799 masked_only, mask)) {
800 SUMA_S_Err("Failed to store output");
801 SUMA_free(fin_orig); fin_orig = NULL; SUMA_free(fout); fout = NULL;
802 if (maskp) SUMA_free(mask); mask=NULL;
803 SUMA_FreeDset(dout); dout=NULL;
804 SUMA_RETURN(dout);
805 }
806
807 if (fin_orig) SUMA_free(fin_orig); fin_orig = NULL;
808 if (fout) SUMA_free(fout); fout = NULL;
809 } /* for k */
810
811 if (!maskp) SUMA_free(mask); mask = NULL;
812 if (NeighbDist != FirstNeighbDist) {
813 SUMA_free2D((char **)NeighbDist, SO->FN->N_Node);
814 }
815 SUMA_RETURN(dout);
816 }
817
818
819 /* Find the extreme points of a dataset */
SUMA_DsetExtrema(SUMA_SurfaceObject * SO,float ** FirstNeighbDist,SUMA_DSET * din,SUMA_DSET * dgrad,float r,float fthresh,float gthresh,byte * maskp,byte mask_by_zeros,SUMA_EXTREMA_DIRECTIONS dir,char * tout)820 SUMA_DSET *SUMA_DsetExtrema(
821 SUMA_SurfaceObject *SO, float **FirstNeighbDist,
822 SUMA_DSET *din, SUMA_DSET *dgrad, float r, float fthresh, float gthresh,
823 byte *maskp, byte mask_by_zeros, SUMA_EXTREMA_DIRECTIONS dir,
824 char *tout)
825 {
826 static char FuncName[]={"SUMA_DsetExtrema"};
827 SUMA_DSET *dout = NULL;
828 int *icols = NULL, N_icols = -1, *ind = NULL, n_incopy=-1, masked_only=0;
829 int k = -1, N_mask=-1, *isrt=NULL, dbg_node=-1;
830 int *mout=NULL, n_peak=0, ni=0, cand=-1, il=0, jl=0, neighb=0;
831 void *ncoli=NULL;
832 char *lblcp=NULL, *s=NULL, Cside=' ', *Sdir=NULL;
833 float *fin=NULL, *gin=NULL, **NeighbDist=NULL, *finsrt=NULL;
834 byte *mask=NULL, *ex_mask=NULL;
835 SUMA_GET_OFFSET_STRUCT *OffS = NULL;
836 FILE *ftable=NULL;
837 SUMA_Boolean LocalHead = NOPE;
838
839 SUMA_ENTRY;
840
841 if (!SO || !din) {
842 SUMA_S_Err("Bad input");
843 SUMA_RETURN(NULL);
844 }
845 if (!SO->FN) {
846 if (!SUMA_SurfaceMetrics_eng(SO, "EdgeList|PolyArea|MemberFace",NULL,
847 0, SUMAg_CF->DsetList)) {
848 fprintf (SUMA_STDERR,
849 "Error %s: Failed in SUMA_SurfaceMetrics.\n", FuncName);
850 }
851 }
852
853 if (SO->Side == SUMA_LEFT) Cside = 'L';
854 else if (SO->Side == SUMA_RIGHT) Cside = 'R';
855
856 switch (dir) {
857 case SUMA_MAXIMUS:
858 Sdir = "max";
859 break;
860 case SUMA_MINIMUS:
861 Sdir = "min";
862 break;
863 case SUMA_EXTREMUS:
864 Sdir = "ext";
865 break;
866 default:
867 SUMA_S_Errv("Bad news for dir %d\n", dir);
868 SUMA_RETURN(NULL);
869 }
870
871 if (tout) {
872 if (!(ftable = fopen(tout,"w"))) {
873 SUMA_S_Errv("Could not open %s for writing\n", tout);
874 SUMA_RETURN(NULL);
875 }
876 fprintf(ftable,
877 "#Output of SurfExtrema. Columns are comma separated.\n"
878 "#Col.1: Node index \n"
879 "#Col.2: Hemisphere side if applicable\n"
880 "#Col.3: Rank of %s\n"
881 "#Col.4: Value at node\n"
882 "#Col.5: Gradient at node\n"
883 "#Col.6: X Y Z coordinates of node. \n"
884 "# Use proper -sv option to relate to volume coords.\n"
885 "#Col.7: Sub-brick index of input\n"
886 "#Col.8: Sub-brick label of input\n"
887 "#To select a set of columns, you can use cut. "
888 "# For example, say you want cols 1,2, and 4:\n"
889 "# \\cut -f '1,2,4' -d ',' %s\n"
890 "# If you just want nodeindex with hemisphere\n"
891 "# label stuck to it (from 3dGroupInCorr's batch mode):\n"
892 "# \\cut -f '1,2' -d ',' %s | \\sed 's/ *, *//g'\n",
893 Sdir, tout, tout);
894 }
895 if (maskp) mask=maskp;
896
897 if (!FirstNeighbDist) {
898 NeighbDist = SUMA_CalcNeighbDist (SO);
899 } else {
900 NeighbDist = FirstNeighbDist;
901 }
902 if (!NeighbDist) {
903 SUMA_S_Err("NULL dists");
904 SUMA_RETURN(NULL);
905 }
906
907 /* what columns can we process ?*/
908 icols = SUMA_FindNumericDataDsetCols(din, &N_icols);
909 if (N_icols <= 0) {
910 SUMA_SL_Err("No approriate data columns in dset");
911 SUMA_RETURN(NOPE);
912 }
913 SUMA_LHv("Have %d numeric columns of input.\n", N_icols);
914 if (!(ind = SDSET_NODE_INDEX_COL(din))) {
915 SUMA_S_Note("Trying to populate the node index element");
916 if (!SUMA_PopulateDsetNodeIndexNel(din, 0)) {
917 SUMA_S_Err("Failed to populate NodeIndex Nel");
918 SUMA_RETURN(NULL);
919 }
920 }
921 /* Create a dset, at least as big as din*/
922 if ((ind = SDSET_NODE_INDEX_COL(din))) {
923 if (!masked_only) {
924 /* preserve all rows */
925 ncoli =
926 SUMA_Copy_Part_Column(ind,
927 NI_rowtype_find_code(SUMA_ColType2TypeCast(SUMA_NODE_INDEX)),
928 SDSET_VECLEN(din), NULL, masked_only, &n_incopy);
929 } else {
930 ncoli =
931 SUMA_Copy_Part_Column(ind,
932 NI_rowtype_find_code(SUMA_ColType2TypeCast(SUMA_NODE_INDEX)),
933 SDSET_VECLEN(din), mask, masked_only, &n_incopy);
934 }
935 if (!ncoli) {
936 SUMA_SL_Err("No index data got copied.");
937 SUMA_RETURN(NULL);
938 }
939 dout = SUMA_CreateDsetPointer("Extrema", SUMA_NODE_BUCKET, NULL,
940 SDSET_IDMDOM(din), n_incopy);
941 if (!SUMA_AddDsetNelCol (dout, NI_get_attribute(din->inel,"COLMS_LABS"), SUMA_NODE_INDEX, ncoli, NULL ,1)) {
942 SUMA_SL_Crit("Failed in SUMA_AddDsetNelCol");
943 SUMA_FreeDset((void*)dout); dout = NULL;
944 SUMA_RETURN(NULL);
945 }
946 if (ncoli) SUMA_free(ncoli); ncoli = NULL;
947 } else {
948 SUMA_S_Err( "Do not have node indices in input dset!\n"
949 " and could not create one.");
950 SUMA_RETURN(NULL);
951 }
952
953 OffS = SUMA_Initialize_getoffsets (SO->N_Node);
954
955 /* Now, for each column, do the dance */
956 for (k=0; k < N_icols; ++k) {
957 lblcp = SUMA_DsetColLabelCopy(din, icols[k], 1);
958
959 /* get a float copy of the data column */
960 if (!(fin = SUMA_DsetCol2FloatFullSortedColumn(din, icols[k], &mask,
961 0.0, SO->N_Node,
962 &N_mask, k==0?YUP:NOPE))){
963 SUMA_S_Err("Failed to extract");
964 SUMA_FreeDset(dout); dout=NULL;
965 if (!maskp) SUMA_free(mask); mask=NULL;
966 SUMA_RETURN(dout);
967 }
968
969 if (LocalHead) {
970 s = SUMA_ShowMeSome(fin, SUMA_float, SO->N_Node, 10, NULL);
971 SUMA_LHv("fin:\n%s\n", s); SUMA_free(s);
972 }
973 if (dgrad) {
974 if (!(gin = SUMA_DsetCol2FloatFullSortedColumn(dgrad, icols[k], &mask,
975 0.0, SO->N_Node,
976 &N_mask, k==0?YUP:NOPE))){
977 SUMA_S_Err("Failed to extract");
978 SUMA_FreeDset(dout); dout=NULL;
979 if (!maskp) SUMA_free(mask); mask=NULL;
980 SUMA_RETURN(dout);
981 }
982 } else {
983 gin = SUMA_AvgGradient(SO, NeighbDist, fin, mask,
984 mask_by_zeros, SUMA_MEAN_GRAD_SCALE);
985 if (!gin) {
986 SUMA_SL_Crit("Failed to compute gradient gin!");
987 SUMA_RETURN(NULL);
988 }
989 }
990
991 /* Now do some work */
992 /* sort the values in fin (lowest first)*/
993 finsrt = (float *)SUMA_malloc(SO->N_Node*sizeof(float));
994 memcpy(finsrt, fin, SO->N_Node*sizeof(float));
995
996 if (dir == SUMA_EXTREMUS) { /* Take absolute value */
997 for (ni=0; ni<SO->N_Node; ++ni) finsrt[ni] = SUMA_ABS(finsrt[ni]);
998 }
999
1000 isrt = SUMA_z_qsort(finsrt, SO->N_Node);
1001 if (LocalHead) {
1002 s = SUMA_ShowMeSome(isrt, SUMA_int, SO->N_Node, 10, NULL);
1003 SUMA_LHv("isrt:\n%s\n", s); SUMA_free(s);
1004 }
1005 /* find top candidate */
1006 ex_mask = (byte *)SUMA_malloc(SO->N_Node*sizeof(byte));
1007 if (!mask) {
1008 memset(ex_mask,1,sizeof(byte)*SO->N_Node);
1009 } else {
1010 memcpy(ex_mask, mask, sizeof(byte)*SO->N_Node);
1011 }
1012 mout = (int*)SUMA_calloc(SO->N_Node, sizeof(int));
1013 n_peak = 0;
1014 ni = SO->N_Node;
1015 cand = -1;
1016 while(ni > 0) {
1017 --ni;
1018 #if 0
1019 if (isrt[ni] == dbg_node) {
1020 fprintf(SUMA_STDERR,"ex_mask[%d]=%d, gin=%f (%f), fin=%f (%f)\n",
1021 isrt[ni], ex_mask[isrt[ni]],
1022 gin[isrt[ni]], gthresh, finsrt[ni], fthresh);
1023 }
1024 #endif
1025 if ( ex_mask[isrt[ni]] && (
1026 ( dir == SUMA_MAXIMUS && /* Maxima */
1027 gin[isrt[ni]] >= gthresh &&
1028 finsrt[ni] >= fthresh ) ||
1029 ( dir == SUMA_MINIMUS && /* minima */
1030 gin[isrt[ni]] <= gthresh &&
1031 finsrt[ni] <= fthresh ) ||
1032 ( dir == SUMA_EXTREMUS && /* absolute */
1033 SUMA_ABS(gin[isrt[ni]]) >= gthresh &&
1034 SUMA_ABS(fin[isrt[ni]]) >= fthresh )
1035 )
1036 ) { /* have good one */
1037 ++n_peak;
1038 cand = isrt[ni];
1039 mout[cand] = n_peak; /* store peak number */
1040 ex_mask[cand] = 0; /* mark it as no longer eligible */
1041 SUMA_LHv("Node %d %c is peak number %d: Val %f, Grad %f,"
1042 " Coords: %f %f %f "
1043 " input %d:%s\n",
1044 cand, Cside, n_peak, fin[isrt[ni]],
1045 gin[isrt[ni]],
1046 SO->NodeList[SO->NodeDim*cand],
1047 SO->NodeList[SO->NodeDim*cand+1],
1048 SO->NodeList[SO->NodeDim*cand+2],
1049 icols[k], lblcp);
1050 if (ftable) {
1051 fprintf(ftable,"%d , %c , %d , %f , %f , %f %f %f , %d , %s\n",
1052 cand, Cside, n_peak, fin[isrt[ni]], gin[isrt[ni]],
1053 SO->NodeList[SO->NodeDim*cand],
1054 SO->NodeList[SO->NodeDim*cand+1],
1055 SO->NodeList[SO->NodeDim*cand+2],
1056 icols[k], lblcp);
1057 }
1058 /* mark all nodes within r as out of contest */
1059 SUMA_getoffsets2 (cand, SO, r, OffS, NULL, 0);
1060 if (r > 0.0) {
1061 for (il=1; il<OffS->N_layers; ++il) {
1062 /* starting at layer 1, layer 0 is the node itself */
1063 for (jl=0; jl<OffS->layers[il].N_NodesInLayer; ++jl) {
1064 neighb = OffS->layers[il].NodesInLayer[jl];
1065 if (OffS->OffVect[neighb] <= r) {
1066 #if 0
1067 if (neighb == dbg_node) {
1068 fprintf(SUMA_STDERR,
1069 "node %d is at dist %f (%f) from candidate node %d\n",
1070 neighb, OffS->OffVect[neighb], r, cand);
1071 }
1072 #endif
1073 ex_mask[neighb] = 0; /* take it out */
1074 }
1075 }
1076 }
1077 } else {
1078 SUMA_S_Err("Not ready to deal with negative r");
1079 SUMA_RETURN(NULL);
1080 }
1081 }
1082 SUMA_Recycle_getoffsets (OffS);
1083 }
1084
1085 lblcp = SUMA_append_replace_string(Sdir, lblcp, "[", 2);
1086 lblcp = SUMA_append_replace_string(lblcp, "]", "", 1);
1087 if (!SUMA_AddDsetNelCol (dout, lblcp, SUMA_NODE_INT,
1088 NULL, NULL ,1)) {
1089 SUMA_S_Crit("Failed to add dset column");
1090 SUMA_RETURN(NULL);
1091 }
1092 SUMA_free(lblcp); lblcp=NULL;
1093 SUMA_LHv("Sticking column %d in dset (mout[0]=%d)\n", k, mout[0]);
1094 if (!SUMA_Vec2DsetCol (dout, k, (void *)mout, SUMA_int,
1095 masked_only, mask)) {
1096 SUMA_S_Err("Failed to store output");
1097 SUMA_free(fin); fin = NULL;
1098 SUMA_free(finsrt); finsrt = NULL;
1099 SUMA_free(mout); mout = NULL;
1100 if (maskp) SUMA_free(mask); mask=NULL;
1101 SUMA_FreeDset(dout); dout=NULL;
1102 SUMA_RETURN(dout);
1103 }
1104
1105 if (mout) SUMA_free(mout); mout=NULL;
1106 if (ex_mask) SUMA_free(ex_mask); ex_mask=NULL;
1107 if (fin) SUMA_free(fin); fin = NULL;
1108 if (finsrt) SUMA_free(finsrt); finsrt = NULL;
1109 if (isrt) SUMA_free(isrt); isrt = NULL;
1110 } /* for k */
1111
1112 if (ftable) fclose(ftable); ftable = NULL;
1113 if (!maskp) SUMA_free(mask); mask = NULL;
1114 if (NeighbDist != FirstNeighbDist) {
1115 SUMA_free2D((char **)NeighbDist, SO->FN->N_Node);
1116 }
1117 if (OffS) SUMA_Free_getoffsets(OffS); OffS = NULL;
1118
1119 SUMA_RETURN(dout);
1120
1121
1122
1123 }
1124
1125
1126 /*!
1127 \brief Finds locally extreme nodes on the surface
1128
1129 \param dothisnode (int) start building from this node
1130 \param ToBeAssigned (float *)
1131 if ToBeAssigned[i] then node i (index into SO's nodelist)
1132 is considered in the clustering and the value at that
1133 node is ToBeAssigned[i]. Vector is SO->N_Node elements
1134 long. Gets modified in function.
1135 \param N_ToBeAssigned (int *) pointer to number of value values
1136 in ToBeAssigned. Gets modified in function.
1137 \param NodeArea(float *) Vector containing area of each node of surface
1138 \param SO (SUMA_SurfaceObject *) the usual deal
1139 \param Opt (SUMA_SURFCLUST_OPTIONS *) options for cluster building.
1140
1141 */
SUMA_Build_Cluster_From_Node_Extrema(int dothisnode,float * ToBeAssigned,int * N_TobeAssigned,float * NodeArea,SUMA_SurfaceObject * SO,SUMA_SURFCLUST_OPTIONS * Opt)1142 SUMA_CLUST_DATUM * SUMA_Build_Cluster_From_Node_Extrema (
1143 int dothisnode,
1144 float *ToBeAssigned, int *N_TobeAssigned, float *NodeArea,
1145 SUMA_SurfaceObject *SO, SUMA_SURFCLUST_OPTIONS *Opt )
1146 {
1147 static char FuncName[]={"SUMA_Build_Cluster_From_Node_Extrema"};
1148 SUMA_CLUST_DATUM *Clust = NULL;
1149 static int ncall;
1150 static int N_Orig = -1;
1151 int il, jl, neighb;
1152 void *dtmp=NULL;
1153 SUMA_GET_OFFSET_STRUCT *OffS = NULL;
1154 DList *offlist = NULL, *candlist=NULL;
1155 DListElmt *elm = NULL, *dothiselm=NULL;
1156 SUMA_OFFSET_LL_DATUM *dat=NULL;
1157 int NewClust = 0;
1158 byte *visited=NULL;
1159 SUMA_Boolean LocalHead = NOPE;
1160
1161 SUMA_ENTRY;
1162
1163 ++ncall;
1164 /* a trick to know when to initialize */
1165 if (Opt->update < 0) { N_Orig = *N_TobeAssigned; Opt->update = -Opt->update; }
1166
1167 if (dothisnode < 0) {
1168 SUMA_SL_Err("Unexpected negative index.");
1169 SUMA_RETURN(NULL);
1170 }
1171
1172 OffS = SUMA_Initialize_getoffsets (SO->N_Node);
1173 Clust = (SUMA_CLUST_DATUM *)SUMA_malloc(sizeof(SUMA_CLUST_DATUM));
1174 Clust->N_Node = 0; Clust->totalarea = 0.0; /* Clust->rank = -1; */
1175 Clust->totalvalue = 0.0; Clust->totalabsvalue = 0.0;
1176 Clust->minvalue = ToBeAssigned[dothisnode]; Clust->minnode = dothisnode;
1177 Clust->maxvalue = ToBeAssigned[dothisnode]; Clust->maxnode = dothisnode;
1178 Clust->varvalue = 0.0; Clust->centralnode = 0;
1179 Clust->weightedcentralnode = 0;
1180 Clust->NodeList = (int *)SUMA_malloc((*N_TobeAssigned) * sizeof(int));
1181 Clust->ValueList = (float *)SUMA_malloc((*N_TobeAssigned) * sizeof(float));
1182 if (!Clust->NodeList || !Clust->ValueList || !OffS) {
1183 SUMA_SL_Crit("Failed to allocate for NodeList or ValueList");
1184 SUMA_free(Clust); Clust = NULL;
1185 SUMA_RETURN(NULL);
1186 }
1187 memset(Clust->com, 0, 3*sizeof(float));
1188 memset(Clust->cen, 0, 3*sizeof(float));
1189 candlist = (DList*)SUMA_malloc(sizeof(DList));
1190 visited = (byte *)SUMA_calloc(SO->N_Node, sizeof(byte));
1191 if (!visited || !candlist) {
1192 SUMA_SL_Crit("Failed to allocate for visited or candlist");
1193 SUMA_free(Clust); Clust = NULL;
1194 SUMA_RETURN(NULL);
1195 }
1196 dlist_init(candlist, NULL);
1197 /* Add node to cluster */
1198 SUMA_ADD_NODE_TO_CLUST(dothisnode, Clust, NodeArea, ToBeAssigned, SO);
1199 /* mark it as assigned, an reduce the number of nodes left to assign*/
1200 ToBeAssigned[dothisnode] = 0; --(*N_TobeAssigned);
1201 visited[dothisnode] = YUP;
1202 dlist_ins_next(candlist, dlist_tail(candlist), (VOID_CAST)dothisnode);
1203 while (*N_TobeAssigned && dlist_size(candlist)) {
1204 /* look in its vicinity */
1205 dothiselm = dlist_head(candlist);
1206 dothisnode = (INT_CAST) dothiselm->data;
1207 SUMA_getoffsets2 (dothisnode, SO, Opt->DistLim, OffS, NULL, 0);
1208 /* remove node from candidate list */
1209 dlist_remove(candlist, dothiselm, (void*)&dtmp);
1210 /* search to see if any are to be assigned */
1211 if (Opt->DistLim >= 0.0) {
1212 for (il=1; il<OffS->N_layers; ++il) {
1213 /* starting at layer 1, layer 0 is the node itself */
1214 for (jl=0; jl<OffS->layers[il].N_NodesInLayer; ++jl) {
1215 neighb = OffS->layers[il].NodesInLayer[jl];
1216 if ( ToBeAssigned[neighb] &&
1217 OffS->OffVect[neighb] <= Opt->DistLim) {
1218 /* take that node into the cluster */
1219 SUMA_ADD_NODE_TO_CLUST( neighb, Clust,
1220 NodeArea, ToBeAssigned, SO);
1221 /* mark it as assigned, an reduce the number
1222 of nodes left to assign*/
1223 ToBeAssigned[neighb] = 0; --(*N_TobeAssigned);
1224 if (Opt->update) {
1225 if (N_Orig - *N_TobeAssigned >= Opt->update) {
1226 if (LocalHead)
1227 fprintf( SUMA_STDERR,
1228 "%s: tick (%d nodes processed)\n",
1229 FuncName, N_Orig - *N_TobeAssigned);
1230 else fprintf(SUMA_STDERR,".");
1231
1232 N_Orig = *N_TobeAssigned;
1233 }
1234 }
1235 /* mark it as a candidate if it has not been visited as
1236 a candidate before */
1237 if (!visited[neighb]) {
1238 dlist_ins_next(candlist, dlist_tail(candlist),
1239 (VOID_CAST)neighb);
1240 visited[neighb] = YUP;
1241 }
1242 }
1243 }
1244 }
1245 } else {/* accept nodes connected by up to -((int)Opt->DistLim) edges */
1246 for (il=1; il <= -((int)Opt->DistLim); ++il) {
1247 for (jl=0; jl<OffS->layers[il].N_NodesInLayer; ++jl) {
1248 neighb = OffS->layers[il].NodesInLayer[jl];
1249 if ( ToBeAssigned[neighb] ) {
1250 /* take that node into the cluster */
1251 SUMA_ADD_NODE_TO_CLUST( neighb, Clust,
1252 NodeArea, ToBeAssigned, SO);
1253 /* mark it as assigned, an reduce the number
1254 of nodes left to assign*/
1255 ToBeAssigned[neighb] = 0; --(*N_TobeAssigned);
1256 if (Opt->update) {
1257 if (N_Orig - *N_TobeAssigned >= Opt->update) {
1258 if (LocalHead)
1259 fprintf( SUMA_STDERR,
1260 "%s: tick (%d nodes processed)\n",
1261 FuncName, N_Orig - *N_TobeAssigned);
1262 else fprintf(SUMA_STDERR,".");
1263
1264 N_Orig = *N_TobeAssigned;
1265 }
1266 }
1267 /* mark it as a candidate if it has not been visited as
1268 a candidate before */
1269 if (!visited[neighb]) {
1270 dlist_ins_next(candlist, dlist_tail(candlist),
1271 (VOID_CAST)neighb);
1272 visited[neighb] = YUP;
1273 }
1274 }
1275 }
1276 }
1277 }
1278 /* recycle */
1279 SUMA_Recycle_getoffsets (OffS);
1280 }
1281 /* free this OffS structure */
1282 if (OffS) SUMA_Free_getoffsets(OffS); OffS = NULL;
1283 if (Opt->update) fprintf(SUMA_STDERR,"\n");
1284 /* destroy the list */
1285 if (candlist) {
1286 dlist_destroy(candlist);
1287 SUMA_free(candlist); candlist = NULL;
1288 }
1289 SUMA_RETURN(Clust);
1290 }
1291
1292 /*!
1293 \brief Finds clusters of data on the surface.
1294
1295 \param SO (SUMA_SurfaceObject *) pointer to surface in question
1296 \param ni (int *) vector of N_Ni node indices that are to be considered
1297 for clustering
1298 \param nv (float *) vector of N_Ni node values (corresponding to ni)
1299 Only non 0 values are considered for clustering
1300 \param N_ni (int) number of values in ni and nv
1301 \param dothisnode (int) index of node (into SO's nodelist, not ni)
1302 to start/proceed from
1303 \param Opt (SUMA_SURFCLUST_OPTIONS *) structure containing clustering options
1304 \param AddToThisClust (SUMA_CLUST_DATUM *) add to this cluster
1305 \param NodeArea (float *) SO->N_Node vector of node areas.
1306 \return ClustList (DList *) list of clusters, in no particular order.
1307 Processing of the list is to be done later.
1308
1309 \sa SUMA_Build_Cluster_From_Node
1310 */
SUMA_FindClusters(SUMA_SurfaceObject * SO,int * ni,float * nv,int N_ni,int dothisnode,SUMA_SURFCLUST_OPTIONS * Opt,float * NodeArea)1311 DList *SUMA_FindClusters ( SUMA_SurfaceObject *SO, int *ni,
1312 float *nv, int N_ni,
1313 int dothisnode, SUMA_SURFCLUST_OPTIONS *Opt,
1314 float *NodeArea)
1315 {
1316 static char FuncName[]={"SUMA_FindClusters"};
1317 DList *list=NULL;
1318 DListElmt *elm=NULL;
1319 float *ToBeAssigned=NULL;
1320 float mean;
1321 int N_n, nc, i, kk, PureNothing=0;
1322 SUMA_CLUST_DATUM *Clust = NULL;
1323 SUMA_Boolean LocalHead = NOPE;
1324
1325 SUMA_ENTRY;
1326
1327 if (!SO || !nv || !ni) {
1328 SUMA_S_Err("Bad parameters");
1329 SUMA_RETURN(list);
1330 }
1331
1332 if (!NodeArea) {
1333 if (!SO->NodeAreas) SO->NodeAreas = SUMA_CalculateNodeAreas(SO, NULL);
1334 NodeArea = SO->NodeAreas;
1335 }
1336
1337 if (!NodeArea) {
1338 SUMA_S_Err("No node areas and failed to get them from SO\n");
1339 SUMA_RETURN(list);
1340 }
1341
1342 if (dothisnode == -1) { /* initialize */
1343 SUMA_LH("Initializing");
1344 nc = 0;
1345 /* initialize the list */
1346 list = (DList *)SUMA_malloc(sizeof(DList));
1347 dlist_init(list, SUMA_FreeClustDatum);
1348 ToBeAssigned = (float *)SUMA_calloc(SO->N_Node, sizeof(float));
1349 N_n = N_ni;
1350 for (i=0; i<N_n; ++i) {
1351 if (!nv[i]) {
1352 ++PureNothing;
1353 }
1354 ToBeAssigned[ni[i]] = nv[i];
1355 }
1356 if (Opt->update) {
1357 fprintf( SUMA_STDERR,
1358 "%s: Have %d nodes to work with. %d nodes have 0 value.\n",
1359 FuncName, N_n, PureNothing);
1360 }
1361 }
1362
1363 while (N_n - PureNothing > 0) {
1364 dothisnode = -1;
1365 for (i=0;i<SO->N_Node; ++i) {
1366 if (ToBeAssigned[i]) {
1367 dothisnode = i; continue;
1368 }
1369 }
1370
1371 if (dothisnode < 0) {
1372 SUMA_S_Err("Not expected here. dothisnode < 0");
1373 SUMA_RETURN(list);
1374 }
1375
1376 if (BuildMethod == SUMA_OFFSETS2_NO_REC) {
1377 SUMA_LH("Building clusts, no rec");
1378 Clust = SUMA_Build_Cluster_From_Node_NoRec( dothisnode,
1379 ToBeAssigned,
1380 &N_n, NodeArea,
1381 SO, Opt);
1382 } else if ( BuildMethod == SUMA_OFFSETS2 ||
1383 BuildMethod == SUMA_OFFSETS_LL ) {
1384 Clust = SUMA_Build_Cluster_From_Node( dothisnode, NULL,
1385 ToBeAssigned,
1386 &N_n, NodeArea,
1387 SO, Opt);
1388 } else {
1389 SUMA_S_Errv("No Such Method (%d)!\n", BuildMethod);
1390 SUMA_DUMP_TRACE("%s", FuncName);
1391 SUMA_RETURN(list);
1392 }
1393 if (!Clust) {
1394 SUMA_SL_Err("Failed in SUMA_Build_Cluster_From_Node*");
1395 SUMA_RETURN(list);
1396 }
1397 if (LocalHead) fprintf( SUMA_STDERR,
1398 "%s: Cluster %p is finished, %d nodes\n",
1399 FuncName, Clust, Clust->N_Node);
1400
1401 if ( (Opt->AreaLim > 0 && Clust->totalarea < Opt->AreaLim) ||
1402 (Opt->NodeLim > 0 && Clust->N_Node < Opt->NodeLim) ||
1403 (Opt->NodeLim < 0 && Opt->AreaLim < 0 &&
1404 Clust->N_Node < -Opt->AreaLim) ) { /* The last option
1405 is to allow users to use a negative
1406 area in the interface to have the effect
1407 of limiting by number of nodes */
1408 SUMA_LH("Cluster less than area (or node number) limit");
1409 SUMA_FreeClustDatum((void *)Clust); Clust = NULL;
1410 } else {
1411 mean = Clust->totalvalue/((float)Clust->N_Node);
1412 for (kk=0; kk < Clust->N_Node; ++kk) {
1413 Clust->varvalue +=
1414 (Clust->ValueList[kk] - mean) * (Clust->ValueList[kk] - mean);
1415 }
1416 if (Clust->N_Node > 1) Clust->varvalue /= (Clust->N_Node - 1);
1417 else Clust->varvalue = 0.0;
1418 /* reallocate to save space */
1419 Clust->NodeList =
1420 (int *)SUMA_realloc(Clust->NodeList, sizeof(int)*Clust->N_Node);
1421 Clust->ValueList =
1422 (float *)SUMA_realloc( Clust->ValueList,
1423 sizeof(float)*Clust->N_Node);
1424 if (!Clust->NodeList || !Clust->ValueList) {
1425 SUMA_SL_Crit("Failed to reallocate for NodeList or ValueList");
1426 SUMA_RETURN(NULL);
1427 }
1428 /* find the central node */
1429 if (Opt->DoCentrality) {
1430 if (Opt->update) {
1431 SUMA_SL_Note( "Looking for central nodes...\n"
1432 "(use -no_cent to skip this slow step)");
1433 }
1434 if (!SUMA_ClusterCenterofMass (SO, Clust, 1)) {
1435 SUMA_SL_Err("Failed to find central node");
1436 SUMA_RETURN(list);
1437 }
1438 }
1439
1440 dlist_ins_next(list, dlist_tail(list), (void *)Clust);
1441 ++nc;
1442 }
1443 }
1444
1445 if (N_n == 0) {
1446 if (LocalHead) fprintf( SUMA_STDERR,
1447 "%s: No more nodes to consider, cleaning up.\n",
1448 FuncName);
1449 if (ToBeAssigned) SUMA_free(ToBeAssigned); ToBeAssigned = NULL;
1450 }
1451
1452 SUMA_RETURN(list);
1453 }
1454
1455 /*! Show the ViewState structure */
SUMA_Show_SurfClust_list(DList * list,FILE * Out,int detail,char * params,char * opts)1456 SUMA_Boolean SUMA_Show_SurfClust_list(
1457 DList *list, FILE *Out, int detail, char *params, char *opts)
1458 {
1459 static char FuncName[]={"SUMA_Show_SurfClust_list"};
1460 char *s = NULL;
1461
1462 SUMA_ENTRY;
1463
1464 if (Out == NULL) Out = stdout;
1465
1466 s = SUMA_Show_SurfClust_list_Info(list, detail, params, opts);
1467 if (!s) {
1468 SUMA_SL_Err("Failed in SUMA_Show_SurfClust_list_Info");
1469 SUMA_RETURN(NOPE);
1470 } else {
1471 fprintf(Out, "%s", s);
1472 SUMA_free(s); s = NULL;
1473 }
1474
1475 SUMA_RETURN(YUP);
1476 }
1477
1478 /*! Show the SurfClust list */
SUMA_Show_SurfClust_list_Info(DList * list,int detail,char * params,char * opts)1479 char *SUMA_Show_SurfClust_list_Info(DList *list, int detail, char *params,
1480 char *opts)
1481 {
1482 static char FuncName[]={"SUMA_Show_SurfClust_list_Info"};
1483 int i, ic, max;
1484 SUMA_STRING *SS = NULL;
1485 DListElmt *elmt=NULL;
1486 SUMA_CLUST_DATUM *cd=NULL;
1487 char *s=NULL, *pad_str, str[20];
1488 int lc[]= { 6, 6, 9, 9, 9, 6, 6, 9, 6, 9, 6, 9, 9, 9, 8, 9, 8,
1489 9, 9, 9, 9, 9, 9 };
1490 char Col[][12] = {
1491 {"# Rank"}, {"num Nd"}, {"Area"}, {"Mean"},
1492 {"|Mean|"},{"Cent"}, {"W Cent"},{"Min V"},
1493 {"Min Nd"}, {"Max V"}, {"Max Nd"} , {"Var"}, {"SEM"},
1494 {"Min |V|"}, {"|Min| Nd"}, {"Max |V|"}, {"|Max| Nd"},
1495 {"COM x"}, {"COM y"}, {"COM z"}, {"Cent x"}, {"Cent y"}, {"Cent z"} };
1496
1497 SUMA_ENTRY;
1498
1499 SS = SUMA_StringAppend (NULL, NULL);
1500
1501 if (!list) {
1502 SS = SUMA_StringAppend (SS,"NULL cluster list.\n");
1503 SUMA_SS2S(SS,s);
1504 SUMA_RETURN(s);
1505 }
1506
1507 if (!list->size) {
1508 SS = SUMA_StringAppend (SS,"Empty cluster list.\n");
1509 SUMA_SS2S(SS,s);
1510 SUMA_RETURN(s);
1511 }else{
1512 if (!opts || !strstr(opts,"No1DColHead")) {
1513 SS = SUMA_StringAppend (SS,"#Col. 0 = Rank \n"
1514 "#Col. 1 = Number of nodes\n"
1515 "#Col. 2 = Total area (units^2)\n"
1516 "#Col. 3 = Mean value\n"
1517 "#Col. 4 = Mean absolute value\n"
1518 "#Col. 5 = Central node\n"
1519 "#Col. 6 = Weighted central node\n"
1520 "#Col. 7 = Minimum value\n"
1521 "#Col. 8 = Minimum node\n"
1522 "#Col. 9 = Maximum value\n"
1523 "#Col. 10 = Maximum node\n"
1524 "#Col. 11 = Variance\n"
1525 "#Col. 12 = Standard error of the mean\n"
1526 "#Col. 13 = Minimum |value|\n"
1527 "#Col. 14 = |Minimum| node\n"
1528 "#Col. 15 = Maximum |value|\n"
1529 "#Col. 16 = |Maximum| node\n"
1530 "#Col. 17 = Center of Mass x\n"
1531 "#Col. 18 = Center of Mass y\n"
1532 "#Col. 19 = Center of Mass z\n"
1533 "#Col. 20 = Centroid x\n"
1534 "#Col. 21 = Centroid y\n"
1535 "#Col. 22 = Centroid z\n"
1536 );
1537 }
1538 SS = SUMA_StringAppend_va (SS,"#Command history:\n"
1539 "#%s\n", params);
1540 for (ic=0; ic<23; ++ic) {
1541 if (ic == 0) sprintf(str,"%s", Col[ic]);
1542 else sprintf(str,"%s", Col[ic]);
1543 pad_str = SUMA_pad_string(str, ' ', lc[ic], 0);
1544 SS = SUMA_StringAppend_va (SS,"%s ", pad_str);
1545 SUMA_free(pad_str);
1546 }
1547 SS = SUMA_StringAppend_va (SS,"\n");
1548 if (detail == 1) {
1549 SS = SUMA_StringAppend_va (
1550 SS,"#Other columns: list of 5 first nodes in ROI.\n");
1551 }
1552 if (detail == 2) {
1553 SS = SUMA_StringAppend_va (
1554 SS,"#Other columns: list all nodes in ROI.\n");
1555 }
1556 if (detail > 0) {
1557 SS = SUMA_StringAppend_va (
1558 SS,"#A total of %d cluster%s were found.\n",
1559 list->size, SUMA_COUNTER_PLURAL(list->size));
1560 }
1561 }
1562
1563 elmt = NULL;
1564 ic = 1;
1565 do {
1566 if (!elmt) elmt = dlist_head(list); else elmt = elmt->next;
1567 if (!elmt) SS = SUMA_StringAppend_va (
1568 SS,"#%d%s cluster element is NULL!\n",
1569 ic, SUMA_COUNTER_SUFFIX(ic));
1570 else {
1571 cd = (SUMA_CLUST_DATUM *)elmt->data;
1572 if (detail > 0) SS = SUMA_StringAppend_va (
1573 SS,"#%d%s cluster\n",
1574 ic, SUMA_COUNTER_SUFFIX(ic));
1575 SS = SUMA_StringAppend_va (SS,
1576 "%6d %6d %9.2f"
1577 " %9.3f %9.3f"
1578 " %6d %6d"
1579 " %9.3f %6d"
1580 " %9.3f %6d"
1581 " %9.3f %9.3f"
1582 " %9.3f %8d"
1583 " %9.3f %8d"
1584 " %9.3f %9.3f %9.3f"
1585 " %9.3f %9.3f %9.3f"
1586 , ic, cd->N_Node, cd->totalarea
1587 , cd->totalvalue/((float)cd->N_Node)
1588 , cd->totalabsvalue/((float)cd->N_Node)
1589 , cd->centralnode, cd->weightedcentralnode
1590 , cd->minvalue, cd->minnode
1591 , cd->maxvalue, cd->maxnode
1592 , cd->varvalue, sqrt(cd->varvalue/cd->N_Node)
1593 , cd->minabsvalue, cd->minabsnode
1594 , cd->maxabsvalue, cd->maxabsnode
1595 , cd->com[0]/cd->totalvalue
1596 , cd->com[1]/cd->totalvalue
1597 , cd->com[2]/cd->totalvalue
1598 , cd->cen[0]/((float)cd->N_Node)
1599 , cd->cen[1]/((float)cd->N_Node)
1600 , cd->cen[2]/((float)cd->N_Node)
1601 );
1602 if (detail > 0) {
1603 if (detail == 1) {
1604 if (cd->N_Node < 5) max = cd->N_Node; else max = 5;
1605 } else max = cd->N_Node;
1606 for (i=0;i<max; ++i)
1607 SS = SUMA_StringAppend_va (SS,"%d\t", cd->NodeList[i]);
1608 }
1609 SS = SUMA_StringAppend(SS,"\n");
1610 }
1611 ++ic;
1612 } while (elmt != dlist_tail(list));
1613
1614 SUMA_SS2S(SS,s);
1615
1616 SUMA_RETURN (s);
1617 }
1618
1619 /*! Turn SurfClust list nimly*/
SUMA_SurfClust_list_2_nel(DList * list,int detail,char * params,char * opts)1620 NI_element *SUMA_SurfClust_list_2_nel(DList *list, int detail, char *params,
1621 char *opts)
1622 {
1623 static char FuncName[]={"SUMA_SurfClust_list_2_nel"};
1624 int i, ic;
1625 DListElmt *elmt=NULL;
1626 SUMA_CLUST_DATUM *cd=NULL;
1627 NI_element *nel=NULL;
1628 char *s=NULL;
1629 float *fv=NULL;
1630 int *iv=NULL;
1631 int tlc[]= { NI_INT, NI_INT, NI_FLOAT, NI_FLOAT,
1632 NI_FLOAT, NI_INT, NI_INT, NI_FLOAT,
1633 NI_INT, NI_FLOAT, NI_INT, NI_FLOAT, NI_FLOAT,
1634 NI_FLOAT, NI_INT, NI_FLOAT, NI_INT };
1635 char Col[][12] = {
1636 {"Rank"}, {"num Nd"}, {"Area"}, {"Mean"},
1637 {"|Mean|"},{"Cent"}, {"W Cent"},{"Min V"},
1638 {"Min Nd"}, {"Max V"}, {"Max Nd"} , {"Var"}, {"SEM"},
1639 {"Min |V|"}, {"|Min| Nd"}, {"Max |V|"}, {"|Max| Nd"} };
1640
1641 SUMA_Boolean LocalHead = NOPE;
1642
1643 SUMA_ENTRY;
1644
1645
1646 if (!list) {
1647 SUMA_RETURN(nel);
1648 }
1649
1650 nel = NI_new_data_element("SurfClust",list->size);
1651 s = SUMA_copy_string(Col[0]);
1652 for (i=1;i<17;++i) {
1653 s = SUMA_append_replace_string(s,Col[i],";",1);
1654 }
1655 NI_set_attribute(nel,"ColumnLabels",s);
1656 NI_set_attribute(nel,"CommandLine",params);
1657
1658
1659 for (i=0;i<17;++i) {
1660 NI_add_column(nel, tlc[i], NULL);
1661 elmt = NULL;
1662 ic = 0;
1663 do {
1664 if (!elmt) elmt = dlist_head(list); else elmt = elmt->next;
1665 if (!elmt) {
1666 SUMA_S_Warnv(" cluster %d element is NULL!\n", ic);
1667 } else {
1668 switch(tlc[i]) {
1669 case NI_FLOAT:
1670 fv = (float*)nel->vec[i];
1671 break;
1672 case NI_INT:
1673 iv = (int*)nel->vec[i];
1674 break;
1675 default:
1676 SUMA_S_Errv("Not ready for type %d, col %d\n", tlc[i], i);
1677 break;
1678 }
1679 cd = (SUMA_CLUST_DATUM *)elmt->data;
1680 switch(i) {
1681 case 0:
1682 iv[ic] = ic+1; break;
1683 case 1:
1684 iv[ic] = cd->N_Node; break;
1685 case 2:
1686 fv[ic] = cd->totalarea; break;
1687 case 3:
1688 fv[ic] = cd->totalvalue/((float)cd->N_Node); break;
1689 case 4:
1690 fv[ic] = cd->totalabsvalue/((float)cd->N_Node); break;
1691 case 5:
1692 iv[ic] = cd->centralnode; break;
1693 case 6:
1694 iv[ic] = cd->weightedcentralnode; break;
1695 case 7:
1696 fv[ic] = cd->minvalue; break;
1697 case 8:
1698 iv[ic] = cd->minnode; break;
1699 case 9:
1700 fv[ic] = cd->maxvalue; break;
1701 case 10:
1702 iv[ic] = cd->maxnode; break;
1703 case 11:
1704 fv[ic] = cd->varvalue; break;
1705 case 12:
1706 fv[ic] = sqrt(cd->varvalue/cd->N_Node); break;
1707 case 13:
1708 fv[ic] = cd->minabsvalue; break;
1709 case 14:
1710 iv[ic] = cd->minabsnode; break;
1711 case 15:
1712 fv[ic] = cd->maxabsvalue; break;
1713 case 16:
1714 iv[ic] = cd->maxabsnode; break;
1715 default:
1716 SUMA_S_Errv("Not ready for column %d\n", i);
1717 break;
1718 }
1719 }
1720 ++ic;
1721 } while (elmt != dlist_tail(list));
1722 }
1723
1724 if (LocalHead) SUMA_ShowNel(nel);
1725 SUMA_RETURN (nel);
1726 }
1727
SUMA_ClustList2Mask(DList * list,int NodeMax)1728 byte *SUMA_ClustList2Mask(DList *list, int NodeMax)
1729 {
1730 static char FuncName[]={"SUMA_ClustList2Mask"};
1731 int j, cnt, crank;
1732 DListElmt *elmt=NULL;
1733 SUMA_CLUST_DATUM *cd=NULL;
1734 byte *ismask=NULL;
1735 SUMA_Boolean LocalHead = NOPE;
1736
1737 SUMA_ENTRY;
1738
1739 if (!list || list->size == 0) {
1740 SUMA_SL_Err("NULL or empty list");
1741 SUMA_RETURN(ismask);
1742 }
1743 /* which nodes in mask ? */
1744 ismask = (byte *)SUMA_calloc(NodeMax, sizeof(byte));
1745 /* you need to allocate for NodeMax (SO->N_Node),
1746 to be safe, otherwise you will need
1747 to search for the highest node index.. */
1748 elmt = NULL; cnt = 0; crank = 0;
1749 do {
1750 ++crank;
1751 if (!elmt) elmt = dlist_head(list);
1752 else elmt = elmt->next;
1753 cd = (SUMA_CLUST_DATUM *)elmt->data;
1754 for (j=0; j<cd->N_Node; ++j) {
1755 if(LocalHead) fprintf (SUMA_STDERR,"nic=%d\t", cd->NodeList[j]);
1756 ismask[cd->NodeList[j]] = (byte)crank;
1757 ++cnt;
1758 }
1759 } while (elmt != dlist_tail(list));
1760 if (LocalHead)
1761 fprintf(SUMA_STDERR,"%s:\n%d nodes in cluster list.\n", FuncName, cnt);
1762
1763 SUMA_RETURN(ismask);
1764
1765 }
1766
1767 /*! Masks a data set by a clustering list*/
SUMA_MaskDsetByClustList(SUMA_DSET * idset,SUMA_SurfaceObject * SO,DList * list,SUMA_Boolean FullList,char * leName)1768 SUMA_DSET *SUMA_MaskDsetByClustList(
1769 SUMA_DSET *idset, SUMA_SurfaceObject *SO,
1770 DList *list, SUMA_Boolean FullList, char *leName)
1771 {
1772 static char FuncName[]={"SUMA_MaskDsetByClustList"};
1773 int i, j;
1774 SUMA_DSET *dset = NULL;
1775 int *ni=NULL, N_ni, cnt;
1776 SUMA_CLUST_DATUM *cd=NULL;
1777 byte *ismask=NULL, *rowmask=NULL, *colmask = NULL;
1778 SUMA_Boolean LocalHead = NOPE;
1779
1780 SUMA_ENTRY;
1781
1782 if (!list || list->size == 0) {
1783 SUMA_SL_Err("NULL or empty list");
1784 SUMA_RETURN(dset);
1785 }
1786
1787 /* which nodes in mask ? */
1788 if (!(ismask = SUMA_ClustList2Mask(list, SO->N_Node))) {
1789 SUMA_S_Err("No nodes in list or bad list");
1790 SUMA_RETURN(dset);
1791 }
1792
1793 /* now form a rowmask vector to parallel rows in idset */
1794 rowmask = (byte *)SUMA_calloc(SDSET_VECLEN(idset), sizeof(byte));
1795 colmask = (byte *)SUMA_calloc(SDSET_VECNUM(idset), sizeof(byte));
1796
1797 /* get the node index column in idset */
1798 ni = SUMA_GetNodeDef(idset);
1799 N_ni = SDSET_VECLEN(idset);
1800 if (!ni) {
1801 SUMA_SL_Err("Failed to find node index column");
1802 SUMA_RETURN(NULL);
1803 }
1804 /* now, fill rowmask */
1805 for (i=0; i<SDSET_VECLEN(idset); ++i) {
1806 if (ismask[ni[i]]) {
1807 rowmask[i]=1;
1808 if(LocalHead) fprintf (SUMA_STDERR,"%d,%d\t", ni[i], i);
1809 }
1810 }
1811 /* fill colmask*/
1812 for (i=0; i<SDSET_VECNUM(idset); ++i) {
1813 if (SUMA_isDsetColumn_inferred(idset, i)) {
1814 colmask[i]=0;
1815 if (LocalHead)
1816 fprintf( SUMA_STDERR,
1817 "%s: Column %d will not be written "
1818 "because it is inferred.\n",
1819 FuncName, i);
1820 } else colmask[i]=1;
1821 }
1822
1823 dset = SUMA_MaskedCopyofDset(idset, rowmask, colmask, !FullList, 1);
1824 if (!dset) {
1825 SUMA_SL_Err("Failed to create masked copy of input");
1826 SUMA_RETURN(NULL);
1827 }
1828
1829 if (rowmask) SUMA_free(rowmask); rowmask = NULL;
1830 if (colmask) SUMA_free(colmask); colmask = NULL;
1831 if (ismask) SUMA_free(ismask); ismask = NULL;
1832
1833 SUMA_RETURN (dset);
1834 }
1835
1836 /*! Turn the clusters to a cluster dataset mask*/
SUMA_SurfClust_list_2_DsetMask(SUMA_SurfaceObject * SO,DList * list,SUMA_Boolean FullList,char * leName)1837 SUMA_DSET *SUMA_SurfClust_list_2_DsetMask(
1838 SUMA_SurfaceObject *SO, DList *list,
1839 SUMA_Boolean FullList, char *leName)
1840 {
1841 static char FuncName[]={"SUMA_SurfClust_list_2_DsetMask"};
1842 int i, ic, max, j, rank;
1843 DListElmt *elmt=NULL;
1844 SUMA_DSET *dset = NULL;
1845 int *NodeIndex=NULL, N_Node, *Val = NULL;
1846 SUMA_CLUST_DATUM *cd=NULL;
1847 SUMA_Boolean LocalHead = NOPE;
1848
1849 SUMA_ENTRY;
1850
1851 if (!list || list->size == 0) {
1852 SUMA_SL_Err("NULL or empty list");
1853 SUMA_RETURN(dset);
1854 }
1855 if (FullList) N_Node = SO->N_Node;
1856 else {
1857 elmt = NULL; N_Node = 0;
1858 do {
1859 if (!elmt) elmt = dlist_head(list);
1860 else elmt = elmt->next;
1861 cd = (SUMA_CLUST_DATUM *)elmt->data;
1862 N_Node += cd->N_Node;
1863 } while (elmt != dlist_tail(list));
1864 }
1865 NodeIndex = (int *)SUMA_malloc(N_Node*sizeof(int));
1866 Val = (int *)SUMA_malloc(N_Node * sizeof(int));
1867 if (!NodeIndex || !Val){
1868 SUMA_SL_Crit("Failed to allocate NodeIndex and or Val");
1869 SUMA_RETURN(dset);
1870 }
1871 if (FullList) {
1872 for (i=0; i<N_Node; ++i) {
1873 NodeIndex[i] = i; Val[i] = 0;
1874 }
1875 elmt = NULL; rank = 1;
1876 do {
1877 if (!elmt) elmt = dlist_head(list);
1878 else elmt = elmt->next;
1879 cd = (SUMA_CLUST_DATUM *)elmt->data;
1880 for (j=0; j<cd->N_Node; ++j) {
1881 Val[cd->NodeList[j]] = rank;
1882 }
1883 ++rank;
1884 } while (elmt != dlist_tail(list));
1885 } else {
1886 elmt = NULL; rank = 1; i = 0;
1887 do {
1888 if (!elmt) elmt = dlist_head(list);
1889 else elmt = elmt->next;
1890 cd = (SUMA_CLUST_DATUM *)elmt->data;
1891 for (j=0; j<cd->N_Node; ++j) {
1892 Val[i] = rank;
1893 NodeIndex[i] = cd->NodeList[j];
1894 ++i;
1895 }
1896 ++rank;
1897 } while (elmt != dlist_tail(list));
1898 }
1899
1900 SUMA_LH("Creating dset pointer");
1901 dset = SUMA_CreateDsetPointer(
1902 leName, /* usually the filename */
1903 SUMA_NODE_ROI, /* mix and match */
1904 NULL, /* no idcode, let function create one from filename*/
1905 NULL, /* no domain str specified */
1906 N_Node /* Number of nodes allocated for */
1907 ); /* DO NOT free dset, it is store in DsetList */
1908
1909 /* form the dataset */
1910 SUMA_LH("Adding NodeDef column ...");
1911 if (!SUMA_AddDsetNelCol (
1912 dset,
1913 "le Node Def",
1914 SUMA_NODE_INDEX,
1915 (void *)NodeIndex,
1916 NULL,
1917 1
1918 )) {
1919 fprintf (stderr,"Error %s:\nFailed in SUMA_AddNelCol", FuncName);
1920 SUMA_RETURN(NULL);
1921 }
1922
1923 if (!SUMA_AddDsetNelCol ( dset, "Cluster Rank",
1924 SUMA_NODE_INT, (void *)Val, NULL ,1)) {
1925 fprintf (stderr,"Error %s:\nFailed in SUMA_AddNelCol", FuncName);
1926 SUMA_RETURN (NULL);
1927 }
1928
1929 if (NodeIndex) SUMA_free(NodeIndex); NodeIndex = NULL;
1930 if (Val) SUMA_free(Val); Val = NULL;
1931
1932 SUMA_RETURN (dset);
1933 }
1934
1935 /*!
1936 \brief Finds a node that best approximates the property of the center of
1937 mass on the surface.
1938 The real center of mass of a curved surface is rarely on the surface,
1939 so this is an attempt at localizing a node that is central to an ROI
1940 and that can reflect the weight, or activity, distribution in that ROI.
1941
1942 \param SO (SUMA_SurfaceObject *)
1943 \param cd (SUMA_CLUST_DATUM *) basically the output of
1944 SUMA_SUMA_Build_Cluster_From_Node or
1945 an element of the list returned by
1946 SUMA_FindClusters
1947 Function will fill centralnode and weightedcentralnode
1948 \param UseSurfDist (int) 0: use distances along the surface
1949 (approximated by distances on the graph)
1950 1: use Euclidian distances.
1951 \return ans (int) : NOPE failed, YUP succeeded.
1952
1953 */
SUMA_ClusterCenterofMass(SUMA_SurfaceObject * SO,SUMA_CLUST_DATUM * cd,int UseSurfDist)1954 int SUMA_ClusterCenterofMass (SUMA_SurfaceObject *SO, SUMA_CLUST_DATUM *cd,
1955 int UseSurfDist)
1956 {
1957 static char FuncName[]={"SUMA_ClusterCenterofMass"};
1958 SUMA_GET_OFFSET_STRUCT *OffS = NULL;
1959 SUMA_OFFSET_STRUCT OS;
1960 int *CoverThisNode = NULL, i, c, ni, nc, centralnode, weightedcentralnode,
1961 nanch, k, WeightByValue = 1;
1962 float Uia[3], dia, Uca[3], dca, s[3], s_w[3], *DistVec=NULL,
1963 *weightedcentrality=NULL, minweightedcentrality = 0.0, *centrality=NULL,
1964 mincentrality = 0.0, mtotal_w, mi_w, fac, mtotal, mi;
1965 static int ncall;
1966 SUMA_Boolean LocalHead = NOPE;
1967
1968 SUMA_ENTRY;
1969
1970 centralnode = -1;
1971 weightedcentralnode = -1;
1972 if (!SO || !cd) { SUMA_RETURN(NOPE); }
1973 if (!cd->N_Node || !cd->NodeList) { SUMA_RETURN(NOPE); }
1974 OffS = SUMA_Initialize_getoffsets (SO->N_Node);
1975 if (!OffS) {
1976 SUMA_SL_Err("Failed to allocate for OffS");
1977 SUMA_RETURN(NOPE);
1978 }
1979
1980 CoverThisNode = (int *)SUMA_calloc(SO->N_Node, sizeof(int));
1981 if (!CoverThisNode) {
1982 SUMA_SL_Crit("Failed to allocate for CoverThisNode");
1983 SUMA_RETURN(NOPE);
1984 }
1985 if (cd->N_Node == SO->N_Node) {
1986 /* not much to do, return 0 */
1987 SUMA_SL_Note("Cluster spans entire surface.\nNo central node.\n");
1988 cd->centralnode = 0;
1989 cd->weightedcentralnode = 0;
1990 SUMA_RETURN(YUP);
1991 }
1992
1993 for (i=0; i<cd->N_Node; ++i) { CoverThisNode[cd->NodeList[i]] = 1; }
1994 nanch = cd->NodeList[0]; /* anchor node */
1995 SUMA_getoffsets2 (cd->NodeList[0], SO, 0, OffS, CoverThisNode, cd->N_Node);
1996 /* help me with a nicer structure (for sanity)*/
1997 if (!SUMA_GetOffset2Offset (OffS, &OS)) {
1998 SUMA_SL_Err("Failed in SUMA_GetOffset2Offset");
1999 SUMA_RETURN(NOPE);
2000 }
2001 if (OffS) SUMA_Free_getoffsets(OffS); OffS = NULL;
2002 /* put the distance in an easy to search array */
2003 DistVec = (float *) SUMA_calloc(SO->N_Node, sizeof(float));
2004 centrality = (float *)SUMA_malloc(OS.N_Neighb * sizeof(float));
2005 weightedcentrality = (float *)SUMA_malloc(OS.N_Neighb * sizeof(float));
2006 if (!DistVec || !centrality || !weightedcentrality) {
2007 SUMA_SL_Crit("Failed to allocate.");
2008 SUMA_RETURN(NOPE);
2009 }
2010 for (c=0; c < OS.N_Neighb; ++c) {
2011 DistVec[OS.Neighb_ind[c]] = OS.Neighb_dist[c];
2012 }
2013
2014 /* Now calculate the center of massity of each node
2015 This is no center of mass in the proper definition of the term*/
2016
2017 #if 1
2018 if (cd->N_Node == 1) {
2019 centralnode = cd->NodeList[0];
2020 weightedcentralnode = cd->NodeList[0];
2021 } else {
2022 for (c=0; c < OS.N_Neighb; ++c) {
2023 nc = OS.Neighb_ind[c]; /* Node index into SO of center node */
2024 s[0] = s[1] = s[2] = 0.0; s_w[0] = s_w[1] = s_w[2] = 0.0;
2025 centrality[c] = 0.0; weightedcentrality[c] = 0.0;
2026 mtotal_w = 0.0; /* recalculated each time, not a big deal ... */
2027 for (i=0; i<cd->N_Node; ++i) {
2028 mi_w = cd->ValueList[i];
2029 mtotal_w += mi_w;
2030 ni = cd->NodeList[i];/*Node index into SO of other node in cluster*/
2031 SUMA_UNIT_VEC((&(SO->NodeList[3*ni])), (&(SO->NodeList[3*nanch])),
2032 Uia, dia);
2033 if (UseSurfDist) {
2034 for (k=0; k<3;++k) {
2035 fac = Uia[k] * DistVec[ni];
2036 s_w[k] += mi_w * fac;
2037 s[k] += fac;
2038 } } else {
2039 for (k=0; k<3;++k) {
2040 fac = Uia[k] * dia;
2041 s_w[k] += mi_w * fac;
2042 s[k] += fac;
2043 } }
2044 }
2045 SUMA_UNIT_VEC((&(SO->NodeList[3*nc])), (&(SO->NodeList[3*nanch])),
2046 Uca, dca);
2047 if (UseSurfDist) {
2048 for (k=0; k<3;++k) {
2049 fac = Uca[k] * DistVec[nc];
2050 s_w[k] -= mtotal_w * fac;
2051 s[k] -= cd->N_Node * fac;
2052 } } else {
2053 for (k=0; k<3;++k) {
2054 fac = Uca[k] * dca;
2055 s_w[k] -= mtotal_w * fac;
2056 s[k] -= cd->N_Node * fac;
2057 } }
2058
2059 SUMA_NORM_VEC(s, 3, centrality[c]);
2060 SUMA_NORM_VEC(s_w, 3, weightedcentrality[c]);
2061 SUMA_LH("call %d, Centrality of node %d / %d = %f , weighted %f\n",
2062 ncall, nc, nanch, centrality[c], weightedcentrality[c]);
2063 if (c == 0) {
2064 mincentrality = centrality[c];
2065 centralnode = nc;
2066 minweightedcentrality = weightedcentrality[c];
2067 weightedcentralnode = nc;
2068 } else {
2069 if (centrality[c] < mincentrality) {
2070 mincentrality = centrality[c];
2071 centralnode = nc;
2072 }
2073 if (weightedcentrality[c] < minweightedcentrality) {
2074 minweightedcentrality = weightedcentrality[c];
2075 weightedcentralnode = nc;
2076 }
2077 }
2078
2079 }
2080 }
2081 cd->centralnode = centralnode;
2082 cd->weightedcentralnode = weightedcentralnode;
2083 SUMA_LH("Central node chosen %d, weighted %d\n",
2084 cd->centralnode, cd->weightedcentralnode);
2085 #else
2086 if (cd->N_Node == 1) {
2087 centralnode = cd->NodeList[0];
2088 } else {
2089 for (c=0; c < OS.N_Neighb; ++c) {
2090 nc = OS.Neighb_ind[c]; /* Node index into SO of center node */
2091 s[0] = s[1] = s[2] = 0.0;
2092 centrality[c] = 0.0;
2093 mtotal = 0.0;/*recalculated each time, not a big deal ... */
2094 for (i=0; i<cd->N_Node; ++i) {
2095 if (WeightByValue) mi = cd->ValueList[i];
2096 else mi = 1.0;
2097 mtotal += mi;
2098 ni = cd->NodeList[i];/*Node index into SO of other node in cluster*/
2099 SUMA_UNIT_VEC((&(SO->NodeList[3*ni])), (&(SO->NodeList[3*nanch])),
2100 Uia, dia);
2101 if (UseSurfDist) {
2102 for (k=0; k<3;++k) s[k] += mi * Uia[k] * DistVec[ni];
2103 } else { for (k=0; k<3;++k) s[k] += mi * Uia[k] * dia; }
2104 }
2105 SUMA_UNIT_VEC((&(SO->NodeList[3*nc])), (&(SO->NodeList[3*nanch])),
2106 Uca, dca);
2107 if (UseSurfDist) {
2108 for (k=0; k<3;++k) s[k] -= mtotal * Uca[k] * DistVec[nc];
2109 } else { for (k=0; k<3;++k) s[k] -= mtotal * Uca[k] * dca; }
2110
2111 SUMA_NORM_VEC(s, 3, centrality[c]);
2112 SUMA_LH("call %d, Centrality of node %d / %d = %f\n",
2113 ncall, nc, nanch, centrality[c]);
2114 if (c == 0) { mincentrality = centrality[c]; centralnode = nc; }
2115 else if (centrality[c] < mincentrality) {
2116 mincentrality = centrality[c]; centralnode = nc;
2117 }
2118 }
2119 }
2120 cd->centralnode = centralnode;
2121 SUMA_LH("Central node chosen %d\n", cd->centralnode);
2122 cd->weightedcentralnode = -1;
2123 #endif
2124
2125 if (CoverThisNode) SUMA_free(CoverThisNode); CoverThisNode = NULL;
2126 if (OS.Neighb_dist) SUMA_free(OS.Neighb_dist); OS.Neighb_dist = NULL;
2127 if (OS.Neighb_ind) SUMA_free(OS.Neighb_ind); OS.Neighb_ind = NULL;
2128 if (DistVec) SUMA_free(DistVec); DistVec = NULL;
2129 if (centrality) SUMA_free(centrality); centrality = NULL;
2130 if (weightedcentrality) SUMA_free(weightedcentrality);
2131 weightedcentrality = NULL;
2132
2133 ++ncall;
2134 SUMA_RETURN(YUP);
2135 }
2136
SUMA_Sort_ClustersList(DList * list,SUMA_SURF_CLUST_SORT_MODES SortMode)2137 SUMA_Boolean SUMA_Sort_ClustersList (DList *list,
2138 SUMA_SURF_CLUST_SORT_MODES SortMode)
2139 {
2140 static char FuncName[]={"SUMA_Sort_ClustersList"};
2141 DListElmt *elmt=NULL, *max_elmt=NULL, *elmt_comp=NULL;
2142 SUMA_CLUST_DATUM *cd=NULL, *cd_comp=NULL, *cd_max=NULL;
2143 int r = 0;
2144 SUMA_Boolean LocalHead = NOPE;
2145
2146 SUMA_ENTRY;
2147
2148 if (!list) { SUMA_RETURN(NOPE); }
2149 switch (SortMode) {
2150 case SUMA_SORT_CLUST_NOT_SET:
2151 SUMA_S_Err("Why are you calling me?");
2152 SUMA_RETURN(NOPE);
2153 break;
2154 case SUMA_SORT_CLUST_NO_SORT:
2155 SUMA_RETURN(YUP);
2156 break;
2157 case SUMA_SORT_CLUST_BY_NUMBER_NODES:
2158 case SUMA_SORT_CLUST_BY_AREA:
2159 elmt = NULL;
2160 do {
2161 if (!elmt) elmt = dlist_head(list);
2162 else elmt = elmt->next;
2163 cd = (SUMA_CLUST_DATUM *)elmt->data;
2164 /* compare to all ahead of this element */
2165 if (elmt != dlist_tail(list)) {
2166 max_elmt = elmt;
2167 cd_max = (SUMA_CLUST_DATUM *)max_elmt->data; elmt_comp = NULL;
2168 do {
2169 if (!elmt_comp) elmt_comp = elmt->next;
2170 else elmt_comp = elmt_comp->next;
2171 cd_comp = (SUMA_CLUST_DATUM *)elmt_comp->data;
2172 if (SortMode == SUMA_SORT_CLUST_BY_NUMBER_NODES) {
2173 if (cd_comp->N_Node > cd_max->N_Node) {
2174 max_elmt = elmt_comp;
2175 cd_max = (SUMA_CLUST_DATUM *)max_elmt->data;
2176 }
2177 }else if (SortMode == SUMA_SORT_CLUST_BY_AREA) {
2178 if (cd_comp->totalarea > cd_max->totalarea) {
2179 max_elmt = elmt_comp;
2180 cd_max = (SUMA_CLUST_DATUM *)max_elmt->data;
2181 }
2182 }
2183 } while (elmt_comp != dlist_tail(list));
2184 if (max_elmt != elmt) { /* max is the real deal */
2185 dlist_remove(list, max_elmt, (void *)(&cd_max));
2186 dlist_ins_prev(list, elmt, (void *)cd_max);
2187 elmt = elmt->prev; /* continue from that one */
2188 }
2189 }
2190 } while (elmt != dlist_tail(list));
2191 SUMA_RETURN(YUP);
2192 break;
2193 default:
2194 break;
2195 }
2196
2197 SUMA_RETURN(YUP);
2198 }
2199
SUMA_NodeClustNumber(SUMA_OVERLAYS * Sover,int node,SUMA_SurfaceObject * SO,SUMA_CLUST_DATUM ** cdp)2200 int SUMA_NodeClustNumber(SUMA_OVERLAYS *Sover, int node,
2201 SUMA_SurfaceObject *SO,
2202 SUMA_CLUST_DATUM **cdp)
2203 {
2204 static char FuncName[]={"SUMA_NodeClustNumber"};
2205 int icl=0;
2206
2207 SUMA_ENTRY;
2208
2209 if (cdp) *cdp = NULL;
2210 if (!Sover || node < 0) SUMA_RETURN(0);
2211 if (!SO) SO = SUMA_SO_of_ColPlane(Sover);
2212 if (!SO || !Sover->ClustOfNode
2213 || !Sover->ClustList) SUMA_RETURN(0);
2214 if (node >= SO->N_Node) SUMA_RETURN(0);
2215 if (cdp && Sover->ClustOfNode[node]) {
2216 *cdp = (SUMA_CLUST_DATUM *)dlist_ith_elmt_data(
2217 Sover->ClustList, Sover->ClustOfNode[node]-1);
2218 if (*cdp == NULL) {
2219 SUMA_S_Errv("Null datum for node %d's cluster %d.\n"
2220 "This should not happen\n",
2221 node, Sover->ClustOfNode[node]);
2222 SUMA_RETURN(0);
2223 }
2224 }
2225 SUMA_RETURN(Sover->ClustOfNode[node]);
2226 }
2227
2228
SUMA_create_SurfClust_Opt(char * init_for)2229 SUMA_SURFCLUST_OPTIONS *SUMA_create_SurfClust_Opt(char *init_for)
2230 {
2231 static char FuncName[]={"SUMA_create_SurfClust_Opt"};
2232 SUMA_SURFCLUST_OPTIONS *Opt=NULL;
2233 int i;
2234
2235 SUMA_ENTRY;
2236
2237 Opt = (SUMA_SURFCLUST_OPTIONS *)
2238 SUMA_calloc(1, sizeof(SUMA_SURFCLUST_OPTIONS));
2239
2240 /* DO NOT CHANGE THESE DEFAULTS, see init_for below */
2241 Opt->out_prefix = NULL;
2242 Opt->oform = SUMA_ASCII_NIML;
2243 Opt->DistLim = -1.5;
2244 Opt->AreaLim = -1.0;
2245 Opt->NodeLim = -1;
2246 Opt->in_name = NULL;
2247 Opt->nodecol = -1;
2248 Opt->labelcol = -1;
2249 Opt->OutROI = NOPE;
2250 Opt->OutClustDset = NOPE;
2251 Opt->FullROIList = NOPE;
2252 Opt->WriteFile = NOPE;
2253 Opt->DoThreshold = SUMA_NO_THRESH;
2254 Opt->ThreshR[0] = Opt->ThreshR[1] = 0.0;
2255 Opt->tind = 0;
2256 Opt->prepend_node_index = NOPE;
2257 Opt->update = 0;
2258 Opt->SortMode = SUMA_SORT_CLUST_NOT_SET;
2259 Opt->DoCentrality = 1;
2260
2261 if (!init_for || !strcmp(init_for, "SurfClust")) {
2262 /* Do nothing stuff above is good to go */
2263 } else if (!strcmp(init_for, "InteractiveClust")) {
2264 Opt->update = -(30000);
2265 Opt->AreaLim = -1.0;
2266 Opt->NodeLim = -1;
2267 Opt->DoCentrality = 0;
2268 Opt->DistLim = -1;
2269 Opt->SortMode = SUMA_SORT_CLUST_BY_AREA;
2270 } else {
2271 SUMA_S_Errv("Don't know how to init for %s\n"
2272 "Returning with defaults for SurfClust program", init_for);
2273 }
2274
2275 SUMA_RETURN(Opt);
2276 }
2277
SUMA_free_SurfClust_Opt(SUMA_SURFCLUST_OPTIONS * Opt)2278 SUMA_SURFCLUST_OPTIONS *SUMA_free_SurfClust_Opt(SUMA_SURFCLUST_OPTIONS *Opt)
2279 {
2280 static char FuncName[]={"SUMA_free_SurfClust_Opt"};
2281 SUMA_ENTRY;
2282 if (Opt) {
2283 if (Opt->out_prefix) SUMA_free(Opt->out_prefix);
2284 SUMA_free(Opt);
2285 }
2286 SUMA_RETURN(NULL);
2287 }
2288
2289 /* copy macro to function */
SUMA_world_stats_node_dbg_f(SUMA_OFFSET_STRUCT * OffS_out,int node,int nval,int nnodes,char * lblcp,byte * nmask,byte strict_mask,float rhood,float * fin_orig,float * fout)2290 static int SUMA_world_stats_node_dbg_f(SUMA_OFFSET_STRUCT * OffS_out,
2291 int node, int nval, int nnodes, char * lblcp,
2292 byte * nmask, byte strict_mask,
2293 float rhood, float * fin_orig, float * fout)
2294 {
2295 static char FuncName[]={"SUMA_world_stats_node_dbg_f"};
2296 int j;
2297
2298 if ( !OffS_out ){
2299 SUMA_S_Err("Ihr Idioten!\nThis debug function is for the offset method!");
2300 return 1;
2301 }
2302
2303 {
2304 FILE *mf ;
2305 char *mfname = SUMA_append_replace_num("SUMA_SurfLocalstat", "_node%ddbg_Col_", node, SUMA_int, 0);
2306 int nj;
2307 mfname = SUMA_append_replace_string(mfname, lblcp, "", 1);
2308 mfname = SUMA_append_replace_string(mfname, ".1D.dset", "", 1);
2309 SUMA_NICEATE_FILENAME(mfname,'\0');
2310 mf=fopen(mfname,"w");
2311 if (!mf) { SUMA_S_Errv("Failed to open %s for writing.\n", mfname); }
2312 else {
2313 fprintf(mf, "#Node %d in mask, total of %d neighbs of which %d went into output of (%f).\n",
2314 node, OffS_out[node].N_Neighb, nval-1, fout[node]);
2315 if (nmask) {
2316 fprintf(mf, "#nmask in use, nmask[node]=%d, strict_mask = %d\n", nmask[node], strict_mask);
2317 fprintf(mf, "#Col. 0: Node index of neighbor (1st row is debug node itself)\n");
2318 fprintf(mf, "#Col. 1: Graph distance of neighbor from debug node\n");
2319 fprintf(mf, "#Col. 2: Neighbor value\n");
2320 fprintf(mf, "#Col. 3: nmask value of neighbor (see strict_mask flag also)\n");
2321 fprintf(mf, "%6d\t%+2.3f\t%+2.3f\t%2d\n", node, 0.0, fin_orig[node], nmask[node]);
2322 } else {
2323 fprintf(mf, "#No masking\n");
2324 fprintf(mf, "#Col. 0: Node index of neighbor (1st row is debug node itself)\n");
2325 fprintf(mf, "#Col. 1: Graph distance of neighbor from debug node\n");
2326 fprintf(mf, "#Col. 2: Neighbor value\n");
2327 fprintf(mf, "%6d\t%+2.3f\t%+2.3f\n", node, 0.0, fin_orig[node]);
2328 }
2329 if (!nmask) {
2330 for (j=0; j<OffS_out[node].N_Neighb; ++j) {
2331 nj = OffS_out[node].Neighb_ind[j];
2332 if (OffS_out[node].Neighb_dist[j] <= rhood) { fprintf(mf, "%6d\t%+2.3f\t%+2.3f\n", nj, OffS_out[node].Neighb_dist[j], fin_orig[nj]); }
2333 }/* for j*/
2334 } else {
2335 for (j=0; j<OffS_out[node].N_Neighb; ++j) {
2336 nj = OffS_out[node].Neighb_ind[j];
2337 if (nmask[nj] || !strict_mask) {
2338 if (OffS_out[node].Neighb_dist[j] <= rhood) { fprintf(mf, "%6d\t%+2.3f\t%+2.3f\t%2d\n", nj, OffS_out[node].Neighb_dist[j], fin_orig[nj], nmask[nj]); }
2339 }
2340 }/* for j*/
2341 }
2342 SUMA_S_Notev("Node %d in mask, total of %d neighbs of which %d went into output of (%f).\nSee also %s\n",
2343 node, OffS_out[node].N_Neighb, nval-1, fout[node], mfname);
2344 SUMA_free(mfname); mfname=NULL;
2345 fclose(mf); mf=NULL;
2346 }
2347 }
2348
2349 return 0;
2350 }
2351
2352 #define SUMA_WORLD_STATS_NODE_DBG { \
2353 if (n == ndbg ) { \
2354 if (OffS_out) { \
2355 FILE *mf ; \
2356 char *mfname = SUMA_append_replace_num("SUMA_SurfLocalstat", "_node%ddbg_Col_", n, SUMA_int, 0); \
2357 mfname = SUMA_append_replace_string(mfname, lblcp, "", 1); \
2358 mfname = SUMA_append_replace_string(mfname, ".1D.dset", "", 1); \
2359 SUMA_NICEATE_FILENAME(mfname,'\0'); \
2360 mf=fopen(mfname,"w"); \
2361 if (!mf) { SUMA_S_Errv("Failed to open %s for writing.\n", mfname); } \
2362 else { \
2363 fprintf(mf, "#Node %d in mask, total of %d neighbs of which %d went into output of (%f).\n", \
2364 n, OffS_out[n].N_Neighb, nval-1, fout[n]); \
2365 if (nmask) {\
2366 fprintf(mf, "#nmask in use, nmask[n]=%d, strict_mask = %d\n", nmask[n], strict_mask); \
2367 fprintf(mf, "#Col. 0: Node index of neighbor (1st row is debug node itself)\n"); \
2368 fprintf(mf, "#Col. 1: Graph distance of neighbor from debug node\n"); \
2369 fprintf(mf, "#Col. 2: Neighbor value\n"); \
2370 fprintf(mf, "#Col. 3: nmask value of neighbor (see strict_mask flag also)\n"); \
2371 fprintf(mf, "%6d\t%+2.3f\t%+2.3f\t%2d\n", n, 0.0, fin_orig[n], nmask[n]); \
2372 } else { \
2373 fprintf(mf, "#No masking\n"); \
2374 fprintf(mf, "#Col. 0: Node index of neighbor (1st row is debug node itself)\n"); \
2375 fprintf(mf, "#Col. 1: Graph distance of neighbor from debug node\n"); \
2376 fprintf(mf, "#Col. 2: Neighbor value\n"); \
2377 fprintf(mf, "%6d\t%+2.3f\t%+2.3f\n", n, 0.0, fin_orig[n]); \
2378 } \
2379 if (!nmask) { \
2380 for (j=0; j<OffS_out[n].N_Neighb; ++j) { \
2381 nj = OffS_out[n].Neighb_ind[j]; \
2382 if (OffS_out[n].Neighb_dist[j] <= rhood) { fprintf(mf, "%6d\t%+2.3f\t%+2.3f\n", nj, OffS_out[n].Neighb_dist[j], fin_orig[nj]); }\
2383 }/* for j*/ \
2384 } else { \
2385 for (j=0; j<OffS_out[n].N_Neighb; ++j) { \
2386 nj = OffS_out[n].Neighb_ind[j]; \
2387 if (nmask[nj] || !strict_mask) { \
2388 if (OffS_out[n].Neighb_dist[j] <= rhood) { fprintf(mf, "%6d\t%+2.3f\t%+2.3f\t%2d\n", nj, OffS_out[n].Neighb_dist[j], fin_orig[nj], nmask[nj]); }\
2389 }\
2390 }/* for j*/ \
2391 } \
2392 SUMA_S_Notev("Node %d in mask, total of %d neighbs of which %d went into output of (%f).\nSee also %s\n", \
2393 n, OffS_out[n].N_Neighb, nval-1, fout[n], mfname);\
2394 SUMA_free(mfname); mfname=NULL;\
2395 fclose(mf); mf=NULL; \
2396 } \
2397 } else { \
2398 SUMA_S_Err("Ihr Idioten!\nThis debug macro is for the offset method!"); \
2399 { \
2400 FILE *mf ; \
2401 char *mfname = SUMA_append_replace_num("SUMA_SurfLocalstat", "_node%ddbg_Col_", n, SUMA_int, 0); \
2402 mfname = SUMA_append_replace_string(mfname, lblcp, "", 1); \
2403 mfname = SUMA_append_replace_string(mfname, ".1D.dset", "", 1); \
2404 SUMA_NICEATE_FILENAME(mfname,'\0'); \
2405 mf=fopen(mfname,"w"); \
2406 if (!mf) { SUMA_S_Errv("Failed to open %s for writing.\n", mfname); } \
2407 else { \
2408 fprintf(mf, "#Node %d in mask, total of %d neighbs in approx mask.\n", \
2409 n, nval);\
2410 for (j=0; j<SO->N_Node; ++j) { if (fwhm_mask[j]) fprintf(mf, "%d\n", j); } \
2411 } \
2412 SUMA_free(mfname); mfname=NULL;\
2413 fclose(mf); mf=NULL; \
2414 } \
2415 }\
2416 } \
2417 }
2418
2419 static double FWHM_MinArea = -1.0;
SUMA_GetFWHM_MinArea(void)2420 double SUMA_GetFWHM_MinArea(void)
2421 {
2422 return(FWHM_MinArea);
2423 }
SUMA_SetFWHM_MinArea(double MinArea)2424 void SUMA_SetFWHM_MinArea(double MinArea)
2425 {
2426 FWHM_MinArea = MinArea;
2427 return;
2428 }
2429
2430 #define FAST_APPROX 1
SUMA_CalculateLocalStats(SUMA_SurfaceObject * SO,SUMA_DSET * din,byte * nmask,byte strict_mask,float rhood,SUMA_OFFSET_STRUCT * UseThisOffset,int ncode,int * code,SUMA_DSET * UseThisDout,int ndbg,SUMA_SurfaceObject * SOf)2431 SUMA_DSET *SUMA_CalculateLocalStats(
2432 SUMA_SurfaceObject *SO, SUMA_DSET *din,
2433 byte *nmask, byte strict_mask,
2434 float rhood, SUMA_OFFSET_STRUCT *UseThisOffset,
2435 int ncode, int *code,
2436 SUMA_DSET *UseThisDout, int ndbg,
2437 SUMA_SurfaceObject *SOf)
2438 {
2439 static char FuncName[]={"SUMA_CalculateLocalStats"};
2440 SUMA_DSET *dout = NULL;
2441 int *icols = NULL, N_icols = -1, *ind = NULL, n_incopy=-1;
2442 int ic = -1, k = -1, n = -1, nj=-1, jj=-1, N_nmask=-1, nval=-1, j=-1;
2443 void *ncoli=NULL;
2444 int masked_only = 1;
2445 char *lblcp=NULL;
2446 float *fin_orig=NULL, *fout = NULL, fp = -1.0;
2447 byte *fwhm_mask=NULL;
2448 SUMA_OFFSET_STRUCT *OffS_out=NULL;
2449 byte *bfull=NULL;
2450 float *NodeAreaVec=NULL, ZoneArea, MinZoneArea;
2451 int ipl;
2452 float Eq[4];
2453 float *SegDist = NULL;
2454 int *mask_record = NULL;
2455 DList *striplist_vec[3];
2456 SUMA_Boolean LocalHead = NOPE;
2457
2458 SUMA_ENTRY;
2459
2460 if (ncode <=0 || !din || rhood <= 0.0f || !code || !SO) {
2461 SUMA_S_Errv("Bad input: SO=%p, din=%p, nmask=%p, "
2462 "rhood=%f, ncode=%d,code=%p, UseThisDout=%p\n",
2463 SO, din, nmask, rhood, ncode, code, UseThisDout);
2464 SUMA_RETURN(NULL);
2465 }
2466
2467 /* what columns can we process ?*/
2468 icols = SUMA_FindNumericDataDsetCols(din, &N_icols);
2469
2470 if (N_icols <= 0) {
2471 SUMA_SL_Err("No approriate data columns in dset");
2472 SUMA_RETURN(NULL);
2473 }
2474
2475
2476 if (UseThisDout) dout = UseThisDout;
2477 else {
2478 if (!(ind = SDSET_NODE_INDEX_COL(din))) {
2479 SUMA_S_Note("Trying to populate the node index element");
2480 if (!SUMA_PopulateDsetNodeIndexNel(din, 0)) {
2481 SUMA_S_Err("Failed to populate NodeIndex Nel");
2482 SUMA_RETURN(NULL);
2483 }
2484 }
2485 /* Create a dset, at least as big as din*/
2486 if ((ind = SDSET_NODE_INDEX_COL(din))) {
2487 if (!masked_only) {
2488 /* preserve all rows */
2489 ncoli =
2490 SUMA_Copy_Part_Column(ind,
2491 NI_rowtype_find_code(SUMA_ColType2TypeCast(SUMA_NODE_INDEX)),
2492 SDSET_VECLEN(din), NULL, masked_only, &n_incopy);
2493 } else {
2494 ncoli =
2495 SUMA_Copy_Part_Column(ind,
2496 NI_rowtype_find_code(SUMA_ColType2TypeCast(SUMA_NODE_INDEX)),
2497 SDSET_VECLEN(din), nmask, masked_only, &n_incopy);
2498 }
2499 if (!ncoli) {
2500 SUMA_SL_Err("No index data got copied.");
2501 SUMA_RETURN(NULL);
2502 }
2503 dout = SUMA_CreateDsetPointer("LocalStat", SUMA_NODE_BUCKET, NULL,
2504 SDSET_IDMDOM(din), n_incopy);
2505 if (!SUMA_AddDsetNelCol (dout, NI_get_attribute(din->inel,"COLMS_LABS"), SUMA_NODE_INDEX, ncoli, NULL ,1)) {
2506 SUMA_SL_Crit("Failed in SUMA_AddDsetNelCol");
2507 SUMA_FreeDset((void*)dout); dout = NULL;
2508 SUMA_RETURN(NULL);
2509 }
2510 if (lblcp) SUMA_free(lblcp); lblcp = NULL;
2511 if (ncoli) SUMA_free(ncoli); ncoli = NULL;
2512 } else {
2513 SUMA_S_Err( "Do not have node indices in input dset!\n"
2514 " and could not create one.");
2515 SUMA_RETURN(NULL);
2516 }
2517 }
2518
2519 /* some checks? Some day? */
2520 if (dout == UseThisDout) {
2521 if (SDSET_VECLEN(dout) != SDSET_VECLEN(din)) {
2522 SUMA_S_Errv("Mismatch in recycled dset (%d rows)"
2523 " and input dset (%d rows)\n",
2524 SDSET_VECLEN(dout), SDSET_VECLEN(din));
2525 SUMA_FreeDset((void*)dout); dout = NULL;
2526 SUMA_RETURN(NULL);
2527 }
2528 if (SDSET_VECNUM(dout) != SDSET_VECNUM(din)) {
2529 SUMA_S_Errv("Mismatch in recycled dset (%d cols) "
2530 "and input dset (%d cols)\n",
2531 SDSET_VECNUM(dout), SDSET_VECNUM(din));
2532 SUMA_FreeDset((void*)dout); dout = NULL;
2533 SUMA_RETURN(NULL);
2534 }
2535 }
2536
2537 /* calculate the areas per node, if need be (a little overhead)*/
2538 MinZoneArea = SUMA_GetFWHM_MinArea() ;
2539 if (MinZoneArea > 0.0) {
2540 NodeAreaVec = SUMA_CalculateNodeAreas(SO, NULL);
2541 }
2542
2543 /* Now, for each code, do the dance */
2544 for (ic=0; ic <ncode; ++ic) {
2545 for (k=0; k < N_icols; ++k) {
2546 /* get a float copy of the data column */
2547 fin_orig = SUMA_DsetCol2Float (din, icols[k], 1);
2548 if (!fin_orig) {
2549 SUMA_SL_Crit("Failed to get copy of column. Woe to thee!");
2550 SUMA_RETURN(NULL);
2551 }
2552 fout = (float *)SUMA_calloc(SO->N_Node, sizeof(float));
2553 if (!fout) {
2554 SUMA_SL_Crit("Failed to allocate fout!");
2555 SUMA_RETURN(NULL);
2556 }
2557 /* make sure column is not sparse, one value per node */
2558 if (k==0) {
2559 SUMA_LH( "Special case k = 0, going to"
2560 " SUMA_MakeSparseColumnFullSorted");
2561 bfull = NULL;
2562 if (!SUMA_MakeSparseColumnFullSorted(
2563 &fin_orig, SDSET_VECFILLED(din), 0.0,
2564 &bfull, din, SO->N_Node)) {
2565 SUMA_S_Err("Failed to get full column vector");
2566 SUMA_RETURN(NULL);
2567 }
2568 if (bfull) {
2569 SUMA_LH(
2570 "Something was filled in SUMA_MakeSparseColumnFullSorted\n" );
2571 if (nmask) { /* combine bfull with nmask */
2572 SUMA_LH( "Merging masks\n" );
2573 for (jj=0; jj < SO->N_Node; ++jj) {
2574 if (nmask[jj] && !bfull[jj]) nmask[jj] = 0;
2575 }
2576 } else { nmask = bfull; }
2577 }
2578 if (nmask) {
2579 N_nmask = 0;
2580 for (n=0; n<SO->N_Node; ++n) { if (nmask[n]) ++ N_nmask; }
2581 SUMA_LHv("FWHMing with node mask (%d nodes in mask)\n", N_nmask);
2582 if (!N_nmask) {
2583 SUMA_S_Warn("Empty mask, nothing to do");
2584 }
2585 }
2586 /* now calculate the neighbor offset structure */
2587 if (!UseThisOffset && code[ic] != NSTAT_FWHMx) { /* no need for
2588 Offset with FWHMx, zones must be big, too big for this method */
2589 SUMA_LH("Calculating OffS_out ...");
2590 OffS_out = SUMA_FormNeighbOffset (SO, rhood, NULL, nmask, -1.0);
2591 } else {
2592 OffS_out = UseThisOffset;
2593 }
2594 } else {
2595 SUMA_LH( "going to SUMA_MakeSparseColumnFullSorted");
2596 if (!SUMA_MakeSparseColumnFullSorted(&fin_orig, SDSET_VECFILLED(din), 0.0, NULL, din, SO->N_Node)) {
2597 SUMA_S_Err("Failed to get full column vector");
2598 SUMA_RETURN(NULL);
2599 }
2600 /* no need for reworking nmask and bfull for each column...*/
2601 }
2602 /* Now I have the data column, nice and solid , do the stats */
2603 switch (code[ic]) {
2604 case NSTAT_MEAN:
2605 {
2606 if( SUMA_nstat_mean(din, icols[k], fin_orig,
2607 nmask, strict_mask, SO->N_Node,
2608 OffS_out, rhood, dout, fout, ndbg) )
2609 SUMA_RETURN(NULL);
2610 break;
2611 }
2612 case NSTAT_MODE:
2613 {
2614 if( SUMA_nstat_mode(din, icols[k], fin_orig,
2615 nmask, strict_mask, SO->N_Node,
2616 OffS_out, rhood, dout, fout, ndbg) )
2617 SUMA_RETURN(NULL);
2618 break;
2619 }
2620 case NSTAT_NUM:
2621 {
2622 if( SUMA_nstat_num(din, icols[k], fin_orig,
2623 nmask, strict_mask, SO->N_Node,
2624 OffS_out, rhood, dout, fout, ndbg) )
2625 SUMA_RETURN(NULL);
2626 break;
2627 }
2628 case -666: /* used to be NSTAT_FWHMx: */
2629 {
2630 lblcp = SUMA_DsetColLabelCopy(din, icols[k], 1);
2631 lblcp = SUMA_append_replace_string("fwhm_", lblcp, "", 2);
2632 if (!SUMA_AddDsetNelCol (dout, lblcp, SUMA_NODE_FLOAT,
2633 (void *)fout, NULL ,1)) {
2634 SUMA_S_Crit("Failed to add dset column");
2635 SUMA_RETURN(NULL);
2636 }
2637 /* form a mask for fwhm function */
2638 if (!(fwhm_mask =
2639 (byte *)SUMA_calloc(SO->N_Node, sizeof(byte)))) {
2640 SUMA_S_Crit("Failed to allocate fwhm_mask");
2641 SUMA_RETURN(NULL);
2642 }
2643 /* create the slice strips, don't want to do this repeatedly */
2644 if (SUMA_Get_UseSliceFWHM()) {
2645 for (ipl=0; ipl<3;++ipl) {
2646 /* get the intersection strips, start alond the various directions */
2647 Eq[0] = Eq[1]=Eq[2]=Eq[3] = 0.0;
2648 Eq[ipl] = 1.0; Eq[3] = -SO->Center[ipl]; /* 0==Saggittal, 1==Coronal, 2==Axial */
2649 SUMA_LHv("Kill me!\nEq:[%f %f %f %f], step: %f\n", Eq[0], Eq[1], Eq[2], Eq[3], SO->EL->AvgLe);
2650 if (!(striplist_vec[ipl] = SUMA_SliceAlongPlane(SO, Eq, SO->EL->AvgLe))) {
2651 SUMA_S_Err("Failed to slice along plane");
2652 SUMA_RETURN(NULL);
2653 }
2654 }
2655 }
2656 if (!nmask) {
2657 SUMA_LH("No mask");
2658 for (n=0; n < SO->N_Node; ++n) {
2659 /* build thy fwhm mask (must have a clean mask here ) */
2660 fwhm_mask[n] = 1; nval = 1;
2661 if (NodeAreaVec) ZoneArea = NodeAreaVec[n]; else ZoneArea = -1.0;
2662 for (j=0; j<OffS_out[n].N_Neighb; ++j) {
2663 nj = OffS_out[n].Neighb_ind[j];
2664 if (OffS_out[n].Neighb_dist[j] <= rhood) {
2665 if (NodeAreaVec) ZoneArea += NodeAreaVec[nj];
2666 fwhm_mask[nj] = 1; ++nval;
2667 }
2668 }/* for j*/
2669 if (ZoneArea > 0.0 && ZoneArea < MinZoneArea) {
2670 if (n==ndbg) {
2671 SUMA_S_Notev("At node %d,\ndata for FWHM estimate spans %g mm2, need %g mm2 for minimum. Bad for node.\n",
2672 n, ZoneArea, MinZoneArea);
2673 }
2674 fout[n] = -1.0;
2675 } else {
2676 if (n==ndbg) {
2677 SUMA_S_Notev("At node %d,\ndata for FWHM estimate spans %g mm2, need %g mm2 for minimum. Looks good.\n",
2678 n, ZoneArea, MinZoneArea);
2679 }
2680 if (n==ndbg) SUMA_SetDbgFWHM(1);
2681 if (SUMA_Get_UseSliceFWHM()) {
2682 fout[n] = SUMA_estimate_slice_FWHM_1dif( SO, fin_orig, fwhm_mask, 1, NULL, striplist_vec);
2683 } else {
2684 fout[n] = SUMA_estimate_FWHM_1dif( SO, fin_orig, fwhm_mask, 1);
2685 }
2686 if (n==ndbg) SUMA_SetDbgFWHM(0);
2687 }
2688 SUMA_WORLD_STATS_NODE_DBG;
2689 /* reset mask */
2690 fwhm_mask[n] = 0;
2691 for (j=0; j<OffS_out[n].N_Neighb; ++j) {
2692 nj = OffS_out[n].Neighb_ind[j]; fwhm_mask[nj] = 0;
2693 }
2694 } /* for n */
2695 } else {
2696 SUMA_LH("Have mask");
2697 if (!strict_mask) {
2698 SUMA_S_Warn("For fwhm, masking must be STRICT!\nProceeding with foolishness.");
2699 }
2700 for (n=0; n < SO->N_Node; ++n) {
2701 if (nmask[n]) {
2702 /* build thy fwhm mask (must have a clean mask here ) */
2703 fwhm_mask[n] = 1; nval = 1;
2704 if (NodeAreaVec) ZoneArea = NodeAreaVec[n]; else ZoneArea = -1.0;
2705 for (j=0; j<OffS_out[n].N_Neighb; ++j) {
2706 nj = OffS_out[n].Neighb_ind[j];
2707 if (nmask[nj] || !strict_mask) {
2708 if (OffS_out[n].Neighb_dist[j] <= rhood) {
2709 if (NodeAreaVec) ZoneArea += NodeAreaVec[nj];
2710 fwhm_mask[nj] = 1; ++nval;
2711 }
2712 }
2713 }/* for j*/
2714 if (ZoneArea > 0.0 && ZoneArea < MinZoneArea) {
2715 if (n==ndbg) {
2716 SUMA_S_Notev("At node %d,\ndata for FWHM estimate spans %g mm2, need %g mm2 for minimum accepted\n",
2717 n, ZoneArea, MinZoneArea);
2718 }
2719 fout[n] = -1.0;
2720 } else {
2721 if (n==ndbg) {
2722 SUMA_S_Notev("At node %d,\ndata for FWHM estimate spans %g mm2, need %g mm2 for minimum. Looks good.\n",
2723 n, ZoneArea, MinZoneArea);
2724 }
2725 if (n==ndbg) SUMA_SetDbgFWHM(1);
2726 if (SUMA_Get_UseSliceFWHM()) {
2727 fout[n] = SUMA_estimate_slice_FWHM_1dif( SO, fin_orig, fwhm_mask, 1, NULL, striplist_vec);
2728 } else {
2729 fout[n] = SUMA_estimate_FWHM_1dif( SO, fin_orig, fwhm_mask, 1);
2730 }
2731 if (n==ndbg) SUMA_SetDbgFWHM(0);
2732 }
2733 SUMA_WORLD_STATS_NODE_DBG;
2734 /* reset mask */
2735 fwhm_mask[n] = 0;
2736 for (j=0; j<OffS_out[n].N_Neighb; ++j) {
2737 nj = OffS_out[n].Neighb_ind[j]; fwhm_mask[nj] = 0;
2738 }
2739 } else {
2740 fout[n] = 0.0; nval = 0;/* Non, rien de rien */
2741 }
2742 } /* for n */
2743 }
2744 if (SUMA_Get_UseSliceFWHM()) {
2745 for (ipl=0; ipl<3; ++ipl) {
2746 SUMA_FREE_DLIST(striplist_vec[ipl]);
2747 striplist_vec[ipl] = NULL; }
2748 }
2749 SUMA_free(fwhm_mask); fwhm_mask = NULL;
2750 SUMA_free(lblcp); lblcp = NULL;
2751 break;
2752 }
2753 case NSTAT_FWHMx:
2754 {
2755 lblcp = SUMA_DsetColLabelCopy(din, icols[k], 1);
2756 lblcp = SUMA_append_replace_string("fwhm_", lblcp, "", 2);
2757 if (!SUMA_AddDsetNelCol (dout, lblcp, SUMA_NODE_FLOAT,
2758 (void *)fout, NULL ,1)) {
2759 SUMA_S_Crit("Failed to add dset column");
2760 SUMA_RETURN(NULL);
2761 }
2762 /* form a mask for fwhm function */
2763 if (!(fwhm_mask = (byte *)SUMA_calloc(SO->N_Node, sizeof(byte)))) {
2764 SUMA_S_Crit("Failed to allocate fwhm_mask");
2765 SUMA_RETURN(NULL);
2766 }
2767
2768 #if FAST_APPROX
2769 if (!SOf && SO->isSphere) SOf = SO;
2770 SegDist = SUMA_SegmentDistortion(SO, SOf);
2771 /* this function should return a vector of 1s if SO == SOf */
2772 mask_record = (int *)SUMA_calloc(SO->N_Node,sizeof(int));
2773 if (!SegDist || !mask_record) { SUMA_S_Crit("Failed to allocate"); SUMA_RETURN(NULL);}
2774 nval = -1; /* Must initialize nval this way for SUMA_APPROX_NEIGHBORS */
2775 #endif
2776 if (!nmask) {
2777 SUMA_LH("No mask");
2778 for (n=0; n < SO->N_Node; ++n) {
2779 #if FAST_APPROX
2780 SUMA_APPROX_NEIGHBORS(SO, SOf, n, rhood, SegDist, mask_record, fwhm_mask, nval);
2781 #else
2782 /* build thy fwhm mask (must have a clean mask here ) */
2783 nval = SUMA_ApproxNeighbors(SO, SOf, n, rhood, fwhm_mask);
2784 #endif
2785
2786 if (n==ndbg) SUMA_SetDbgFWHM(1);
2787 if (nval > 6) fout[n] = SUMA_estimate_FWHM_1dif( SO, fin_orig, fwhm_mask, 1);
2788 else fout[n] = -1.0;
2789 if (n==ndbg) SUMA_SetDbgFWHM(0);
2790 SUMA_WORLD_STATS_NODE_DBG;
2791
2792 } /* for n */
2793 } else {
2794 SUMA_LH("Have mask");
2795 if (!strict_mask) {
2796 SUMA_S_Warn("For fwhm, masking must be STRICT!\nProceeding with foolishness.");
2797 }
2798 for (n=0; n < SO->N_Node; ++n) {
2799 if (nmask[n]) {
2800 #if FAST_APPROX
2801 SUMA_APPROX_NEIGHBORS(SO, SOf, n, rhood, SegDist, mask_record, fwhm_mask, nval);
2802 #else
2803 /* build thy fwhm mask (must have a clean mask here ) */
2804 nval = SUMA_ApproxNeighbors(SO, SOf, n, rhood, fwhm_mask);
2805 #endif
2806 if (nval > 6) {
2807 /* Now qualify the fwhm_mask to include only nodes in nmask, if strict_mask */
2808 if (strict_mask) {
2809 for (j=0; j<SO->N_Node; ++j) {
2810 fwhm_mask[j] *= nmask[j];
2811 }
2812 }
2813 if (n==ndbg) SUMA_SetDbgFWHM(1);
2814 fout[n] = SUMA_estimate_FWHM_1dif( SO, fin_orig, fwhm_mask, 1);
2815 if (n==ndbg) SUMA_SetDbgFWHM(0);
2816 } else {
2817 fout[n] = -1.0;
2818 }
2819 SUMA_WORLD_STATS_NODE_DBG;
2820 } else {
2821 fout[n] = 0.0; nval = 0;/* Non, rien de rien */
2822 }
2823 } /* for n */
2824 }
2825 #if FAST_APPROX
2826 if (SegDist) SUMA_free(SegDist); SegDist = NULL;
2827 if (mask_record) SUMA_free(mask_record); mask_record = NULL;
2828 #else
2829 SUMA_ApproxNeighbors(NULL, NULL, -1, 0, NULL); /* cleanup ApproxNeighbors */
2830 #endif
2831 SUMA_free(fwhm_mask); fwhm_mask = NULL;
2832 SUMA_free(lblcp); lblcp = NULL;
2833 break;
2834 }
2835 default:
2836 {
2837 SUMA_S_Errv("Should not be here, not ready to deal with this stat (%d)\n", code[ic]);
2838 SUMA_RETURN(NULL);
2839 break;
2840 }
2841 }
2842
2843 /* add this column to the output dset */
2844 if (!SUMA_Float2DsetCol (dout, icols[k], fout, 1, nmask)) {
2845 SUMA_S_Err("Failed to update dset's values");
2846 SUMA_RETURN(NULL);
2847 }
2848
2849 if (fin_orig) SUMA_free(fin_orig); fin_orig = NULL;
2850 if (fout) SUMA_free(fout); fout = NULL;
2851 } /* for k */
2852 /* Pre Dec 06 stupidity: if (bfull == nmask) { if (nmask) SUMA_free(nmask); nmask = NULL; bfull = NULL; } */
2853 if (bfull) SUMA_free(bfull); bfull = NULL;
2854
2855 }/* for ic */
2856 if (NodeAreaVec) SUMA_free(NodeAreaVec); NodeAreaVec = NULL;
2857 if (!UseThisOffset) OffS_out = SUMA_free_NeighbOffset (SO, OffS_out);
2858
2859 SUMA_RETURN(dout);
2860 }
2861
2862
2863 /* if strict mask, skip neighbors not in mask (generally 1) */
SUMA_nstat_mean(SUMA_DSET * din,int icol,float * fin_orig,byte * nmask,byte strict_mask,int nnode,SUMA_OFFSET_STRUCT * OffS_out,float rhood,SUMA_DSET * dout,float * fout,int ndbg)2864 static int SUMA_nstat_mean(SUMA_DSET *din, int icol, float *fin_orig,
2865 byte *nmask, byte strict_mask, int nnode,
2866 SUMA_OFFSET_STRUCT *OffS_out, float rhood,
2867 SUMA_DSET *dout, float * fout, int ndbg)
2868 {
2869 static char FuncName[]={"SUMA_nstat_mean"};
2870 char *lblcp=NULL;
2871 int node, bind, bnode; /* node index, neighb index, neighb node */
2872 int fp=-1, nval=0; /* float value, num values */
2873 SUMA_Boolean LocalHead=NOPE;
2874
2875 SUMA_ENTRY;
2876
2877 /* verify that we have the expected inputs */
2878 if( !din || !dout || !OffS_out || !fout ) {
2879 SUMA_S_Errv("Bad nstat_mean data inputs: din=%p, dout=%p, "
2880 "OffS_out=%p, fout=%p\n",
2881 din, dout, OffS_out, fout);
2882 SUMA_RETURN(1);
2883 }
2884
2885 if( icol < 0 || rhood < 0.0f || nnode < 0 ) {
2886 SUMA_S_Errv("Bad nstat_mean val inputs: icol=%d, rhood=%f, nnode=%d\n",
2887 icol, rhood, nnode);
2888 SUMA_RETURN(1);
2889 }
2890
2891 lblcp = SUMA_DsetColLabelCopy(din, icol, 1);
2892 lblcp = SUMA_append_replace_string("mean_", lblcp, "", 2);
2893 if (!SUMA_AddDsetNelCol(dout, lblcp, SUMA_NODE_FLOAT,
2894 (void *)fout, NULL ,1)) {
2895 SUMA_S_Crit("Failed to add dset column");
2896 SUMA_RETURN(1);
2897 }
2898
2899 /* SUMA_LH ends in '}', so skip ';' or enclose in {} */
2900 if (nmask) { SUMA_LH("Have mask"); }
2901 else { SUMA_LH("No mask"); }
2902
2903 for (node=0; node < nnode; ++node) {
2904 /* skip node (pass on value) if we are not in mask */
2905 if( nmask && !nmask[node] ) {
2906 fout[node] = fin_orig[node];
2907 continue;
2908 }
2909
2910 fp = fin_orig[node];
2911 nval = 1;
2912 for (bind=0; bind<OffS_out[node].N_Neighb; ++bind) {
2913 bnode = OffS_out[node].Neighb_ind[bind];
2914
2915 /* if strict mask, skip neighbors not in mask */
2916 if ( strict_mask && nmask && !nmask[bnode] )
2917 continue;
2918 if (OffS_out[node].Neighb_dist[bind] <= rhood) {
2919 fp += fin_orig[bnode]; ++nval; }
2920 }/* for bind*/
2921
2922 fout[node] = fp/(float)nval;
2923
2924 if (node == ndbg )
2925 SUMA_world_stats_node_dbg_f(OffS_out, node, nval, nnode, lblcp,
2926 nmask, strict_mask, rhood, fin_orig, fout);
2927
2928 } /* for node */
2929 SUMA_free(lblcp); lblcp = NULL;
2930
2931 SUMA_RETURN(0);
2932 }
2933
2934
2935 /* count the number of nodes applied in any computation */
2936 /* if strict mask, skip neighbors not in mask (generally 1) */
SUMA_nstat_num(SUMA_DSET * din,int icol,float * fin_orig,byte * nmask,byte strict_mask,int nnode,SUMA_OFFSET_STRUCT * OffS_out,float rhood,SUMA_DSET * dout,float * fout,int ndbg)2937 static int SUMA_nstat_num( SUMA_DSET *din, int icol, float *fin_orig,
2938 byte *nmask, byte strict_mask, int nnode,
2939 SUMA_OFFSET_STRUCT *OffS_out, float rhood,
2940 SUMA_DSET *dout, float * fout, int ndbg)
2941 {
2942 static char FuncName[]={"SUMA_nstat_num"};
2943 char *lblcp=NULL;
2944 int node, bind, bnode; /* node index, neighb index, neighb node */
2945 int nval=0; /* num values */
2946 SUMA_Boolean LocalHead=NOPE;
2947
2948 SUMA_ENTRY;
2949
2950 /* verify that we have the expected inputs */
2951 if( !din || !dout || !OffS_out || !fout ) {
2952 SUMA_S_Errv("Bad nstat_mean data inputs: din=%p, dout=%p, "
2953 "OffS_out=%p, fout=%p\n",
2954 din, dout, OffS_out, fout);
2955 SUMA_RETURN(1);
2956 }
2957
2958 if( icol < 0 || rhood < 0.0f || nnode < 0 ) {
2959 SUMA_S_Errv("Bad nstat_mean val inputs: icol=%d, rhood=%f, nnode=%d\n",
2960 icol, rhood, nnode);
2961 SUMA_RETURN(1);
2962 }
2963
2964 lblcp = SUMA_DsetColLabelCopy(din, icol, 1);
2965 lblcp = SUMA_append_replace_string("mean_", lblcp, "", 2);
2966 if (!SUMA_AddDsetNelCol(dout, lblcp, SUMA_NODE_FLOAT,
2967 (void *)fout, NULL ,1)) {
2968 SUMA_S_Crit("Failed to add dset column");
2969 SUMA_RETURN(1);
2970 }
2971
2972 /* SUMA_LH ends in '}', so skip ';' or enclose in {} */
2973 if (nmask) { SUMA_LH("Have mask"); }
2974 else { SUMA_LH("No mask"); }
2975
2976 for (node=0; node < nnode; ++node) {
2977 /* skip node (pass on value) if we are not in mask */
2978 if( nmask && !nmask[node] ) {
2979 fout[node] = fin_orig[node];
2980 continue;
2981 }
2982
2983 nval = 1;
2984 for (bind=0; bind<OffS_out[node].N_Neighb; ++bind) {
2985 bnode = OffS_out[node].Neighb_ind[bind];
2986
2987 /* if strict mask, skip neighbors not in mask */
2988 if ( strict_mask && nmask && !nmask[bnode] )
2989 continue;
2990 if (OffS_out[node].Neighb_dist[bind] <= rhood)
2991 ++nval;
2992
2993 }/* for bind*/
2994
2995 fout[node] = (float)nval;
2996
2997 if (node == ndbg )
2998 SUMA_world_stats_node_dbg_f(OffS_out, node, nval, nnode, lblcp,
2999 nmask, strict_mask, rhood, fin_orig, fout);
3000
3001 } /* for node */
3002 SUMA_free(lblcp); lblcp = NULL;
3003
3004 SUMA_RETURN(0);
3005 }
3006
3007
3008 /* if strict mask, skip neighbors not in mask (generally 1) */
SUMA_nstat_mode(SUMA_DSET * din,int icol,float * fin_orig,byte * nmask,byte strict_mask,int nnode,SUMA_OFFSET_STRUCT * OffS_out,float rhood,SUMA_DSET * dout,float * fout,int ndbg)3009 static int SUMA_nstat_mode(SUMA_DSET *din, int icol, float *fin_orig,
3010 byte *nmask, byte strict_mask, int nnode,
3011 SUMA_OFFSET_STRUCT *OffS_out, float rhood,
3012 SUMA_DSET *dout, float * fout, int ndbg)
3013 {
3014 static char FuncName[]={"SUMA_nstat_mode"};
3015 float_list flist;
3016 char *lblcp=NULL;
3017 int node, bind, bnode; /* node index, neighb index, neighb node */
3018 int nvals=0, index=0;
3019 SUMA_Boolean LocalHead=NOPE;
3020
3021 SUMA_ENTRY;
3022
3023 /* verify that we have the expected inputs */
3024 if( !din || !dout || !OffS_out || !fout ) {
3025 SUMA_S_Errv("Bad nstat_mean data inputs: din=%p, dout=%p, "
3026 "OffS_out=%p, fout=%p\n",
3027 din, dout, OffS_out, fout);
3028 SUMA_RETURN(1);
3029 }
3030
3031 if( icol < 0 || rhood < 0.0f || nnode < 0 ) {
3032 SUMA_S_Errv("Bad nstat_mean val inputs: icol=%d, rhood=%f, nnode=%d\n",
3033 icol, rhood, nnode);
3034 SUMA_RETURN(1);
3035 }
3036
3037 lblcp = SUMA_DsetColLabelCopy(din, icol, 1);
3038 lblcp = SUMA_append_replace_string("mean_", lblcp, "", 2);
3039 if (!SUMA_AddDsetNelCol(dout, lblcp, SUMA_NODE_FLOAT,
3040 (void *)fout, NULL ,1)) {
3041 SUMA_S_Crit("Failed to add dset column");
3042 SUMA_RETURN(1);
3043 }
3044
3045 /* prepare float list for mode, size can be modified later */
3046 init_float_list(&flist, 256);
3047
3048 /* SUMA_LH ends in '}', so skip ';' or enclose in {} */
3049 if (nmask) { SUMA_LH("Have mask"); }
3050 else { SUMA_LH("No mask"); }
3051
3052 for (node=0; node < nnode; ++node) {
3053 /* skip node (pass on value) if we are not in mask */
3054 if( nmask && !nmask[node] ) {
3055 fout[node] = fin_orig[node];
3056 continue;
3057 }
3058
3059 /* make list of values, sort, find most common */
3060 clear_float_list(&flist);
3061
3062 add_to_float_list(&flist, fin_orig[node], 256);
3063 for (bind=0; bind<OffS_out[node].N_Neighb; ++bind) {
3064 bnode = OffS_out[node].Neighb_ind[bind];
3065
3066 /* if strict mask, skip neighbors not in mask */
3067 if ( strict_mask && nmask && !nmask[bnode] )
3068 continue;
3069 if (OffS_out[node].Neighb_dist[bind] <= rhood) {
3070 add_to_float_list(&flist, fin_orig[bnode], 256);
3071 }
3072 }/* for bind*/
3073
3074 float_list_comp_mode( &flist, fout+node, &nvals, &index );
3075
3076 if (node == ndbg ) {
3077 SUMA_world_stats_node_dbg_f(OffS_out, node, fout[node], nnode, lblcp,
3078 nmask, strict_mask, rhood, fin_orig, fout);
3079 fprintf(SUMA_STDERR, "++ dnode %d: mode %f, nvals %d, index %d\n",
3080 node, fout[node], nvals, index);
3081 }
3082
3083 } /* for node */
3084 SUMA_free(lblcp); lblcp = NULL;
3085
3086 SUMA_RETURN(0);
3087 }
3088
3089 /*----------------------------------------------------------------------
3090 * float_list_comp_mode - sort list and compute mode
3091 *
3092 * return 0 : on success
3093 * -1 : on error
3094 *----------------------------------------------------------------------
3095 */
float_list_comp_mode(float_list * f,float * mode,int * nvals,int * index)3096 static int float_list_comp_mode( float_list *f, float *mode, int *nvals,
3097 int *index )
3098 {
3099 static char FuncName[]={"float_list_comp_mode"};
3100 float fcur;
3101 int ncur, c;
3102
3103 SUMA_ENTRY;
3104
3105 if( f->num <= 0 ) SUMA_RETURN(0);
3106
3107 /* start by sorting the list */
3108 if( f->num > 1 ) qsort_float(f->num, f->list);
3109
3110 /* init default results */
3111 *nvals = ncur = 1;
3112 *mode = fcur = f->list[0];
3113 *index = 0;
3114
3115 for ( c = 1; c < f->num; c++ )
3116 {
3117 if ( f->list[c] == fcur )
3118 ncur ++;
3119 else /* found a new entry to count */
3120 {
3121 if ( ncur > *nvals ) /* keep track of any new winner */
3122 {
3123 *mode = fcur;
3124 *nvals = ncur;
3125 *index = c - ncur; /* note first occurrence */
3126 }
3127
3128 fcur = f->list[c];
3129 ncur = 1;
3130 }
3131 }
3132
3133 if ( ncur > *nvals ) /* keep track of any new winner */
3134 {
3135 *mode = fcur;
3136 *nvals = ncur;
3137 *index = c - ncur; /* note first occurrence */
3138 }
3139
3140 SUMA_RETURN(0);
3141 }
3142