00001 #include "SUMA_suma.h"
00002
00003 #undef STAND_ALONE
00004
00005 #if defined SUMA_SurfClust_STANDALONE
00006 #define STAND_ALONE
00007 #endif
00008
00009 #ifdef STAND_ALONE
00010
00011 SUMA_SurfaceViewer *SUMAg_cSV = NULL;
00012 SUMA_SurfaceViewer *SUMAg_SVv = NULL;
00013
00014 int SUMAg_N_SVv = 0;
00015 SUMA_DO *SUMAg_DOv = NULL;
00016 int SUMAg_N_DOv = 0;
00017 SUMA_CommonFields *SUMAg_CF = NULL;
00018 #else
00019 extern SUMA_CommonFields *SUMAg_CF;
00020 extern SUMA_DO *SUMAg_DOv;
00021 extern SUMA_SurfaceViewer *SUMAg_SVv;
00022 extern int SUMAg_N_SVv;
00023 extern int SUMAg_N_DOv;
00024 #endif
00025
00026 static int BuildMethod;
00027 typedef enum { SUMA_NO_BUILD_METHOD, SUMA_OFFSETS2, SUMA_OFFSETS_LL, SUMA_OFFSETS2_NO_REC } SUMA_CLUST_BUILD_METHODS;
00028
00029 void SUMA_FreeClustDatum (void * data)
00030 {
00031 static char FuncName[]={"SUMA_FreeClustDatum"};
00032 SUMA_CLUST_DATUM *Clust = NULL;
00033 SUMA_ENTRY;
00034
00035 if (!data) SUMA_RETURNe;
00036 Clust = (SUMA_CLUST_DATUM *)data;
00037 if (Clust->NodeList) SUMA_free(Clust->NodeList);
00038 if (Clust->ValueList) SUMA_free(Clust->ValueList);
00039 SUMA_free(Clust);
00040
00041 SUMA_RETURNe;
00042 }
00043
00044
00045
00046
00047
00048
00049
00050
00051 float *SUMA_CalculateNodeAreas(SUMA_SurfaceObject *SO)
00052 {
00053 static char FuncName[]={"SUMA_CalculateNodeAreas"};
00054 float *NodeAreas=NULL;
00055 int *flist = NULL, i, c;
00056
00057 SUMA_ENTRY;
00058
00059 if (!SO) { SUMA_RETURN(NodeAreas); }
00060 if (!SO->PolyArea) {
00061 if (!SUMA_SurfaceMetrics_eng(SO, "PolyArea", NULL, 0, SUMAg_CF->DsetList)) {
00062 fprintf (SUMA_STDERR,"Error %s: Failed in SUMA_SurfaceMetrics.\n", FuncName);
00063 SUMA_RETURN(NodeAreas);
00064 }
00065 }
00066
00067 NodeAreas = (float *)SUMA_malloc(SO->N_Node*sizeof(float));
00068 if (!NodeAreas) { SUMA_SL_Crit ("Failed to allocate for NodeAreas"); SUMA_RETURN(NodeAreas); }
00069
00070 for (i=0; i<SO->N_Node; ++i) {
00071 flist = SO->MF->NodeMemberOfFaceSet[i];
00072 NodeAreas[i] = 0.0;
00073 for (c = 0; c < SO->MF->N_Memb[i]; c++) {
00074 NodeAreas[i] += SO->PolyArea[flist[c]];
00075 }
00076 NodeAreas[i] /= 3.0;
00077 }
00078
00079 SUMA_RETURN(NodeAreas);
00080 }
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 SUMA_CLUST_DATUM * SUMA_Build_Cluster_From_Node(int dothisnode, SUMA_CLUST_DATUM *AddToThisClust,
00106 float *ToBeAssigned, int *N_TobeAssigned, float *NodeArea,
00107 SUMA_SurfaceObject *SO, SUMA_SURFCLUST_OPTIONS *Opt)
00108 {
00109 static char FuncName[]={"SUMA_Build_Cluster_From_Node"};
00110 SUMA_CLUST_DATUM *Clust = NULL;
00111 static int ncall;
00112 int il, jl, neighb;
00113 SUMA_GET_OFFSET_STRUCT *OffS = NULL;
00114 DList *offlist = NULL;
00115 DListElmt *elm = NULL;
00116 SUMA_OFFSET_LL_DATUM *dat=NULL;
00117 int NewClust = 0;
00118 SUMA_Boolean LocalHead = NOPE;
00119
00120 SUMA_ENTRY;
00121 ++ncall;
00122 if (dothisnode < 0) {
00123 SUMA_SL_Err("Unexpected negative index.");
00124 SUMA_RETURN(NULL);
00125 }
00126 if (!AddToThisClust) {
00127 Clust = (SUMA_CLUST_DATUM *)SUMA_malloc(sizeof(SUMA_CLUST_DATUM));
00128 Clust->N_Node = 0; Clust->totalarea = 0.0;
00129 Clust->totalvalue = 0.0; Clust->totalabsvalue = 0.0;
00130 Clust->minvalue = ToBeAssigned[dothisnode]; Clust->minnode = dothisnode;
00131 Clust->maxvalue = ToBeAssigned[dothisnode]; Clust->maxnode = dothisnode;
00132 Clust->varvalue = 0.0; Clust->centralnode = 0; Clust->weightedcentralnode = 0;
00133 Clust->NodeList = (int *)SUMA_malloc((*N_TobeAssigned) * sizeof(int));
00134 Clust->ValueList = (float *)SUMA_malloc((*N_TobeAssigned) * sizeof(float));
00135 if (!Clust->NodeList || !Clust->ValueList) {
00136 SUMA_SL_Crit("Failed to allocate for NodeList or ValueList");
00137 SUMA_free(Clust); Clust = NULL;
00138 SUMA_RETURN(NULL);
00139 }
00140 NewClust = 1;
00141 if (LocalHead) fprintf (SUMA_STDERR,"%s: New Cluster %p, with node %d\n", FuncName, Clust, dothisnode);
00142 } else {
00143 NewClust = 0;
00144 Clust = AddToThisClust;
00145 if (LocalHead) fprintf (SUMA_STDERR,"%s: Reusing Cluster %p, with node %d\n", FuncName, Clust, dothisnode);
00146 }
00147
00148
00149 if (LocalHead) fprintf(SUMA_STDERR,"%s: Adding node %d to cluster %p of %d nodes\n", FuncName, dothisnode, Clust, Clust->N_Node);
00150 Clust->NodeList[Clust->N_Node] = dothisnode;
00151 Clust->totalarea += NodeArea[dothisnode];
00152 Clust->totalvalue += ToBeAssigned[dothisnode];
00153 Clust->totalabsvalue += (float)fabs((float)ToBeAssigned[dothisnode]);
00154 if (ToBeAssigned[dothisnode] < Clust->minvalue) { Clust->minvalue = ToBeAssigned[dothisnode]; Clust->minnode = dothisnode; }
00155 if (ToBeAssigned[dothisnode] > Clust->maxvalue) { Clust->maxvalue = ToBeAssigned[dothisnode]; Clust->maxnode = dothisnode; }
00156 Clust->ValueList[Clust->N_Node] = ToBeAssigned[dothisnode];
00157 ++Clust->N_Node;
00158
00159
00160 ToBeAssigned[dothisnode] = 0; --(*N_TobeAssigned);
00161
00162 if (BuildMethod == SUMA_OFFSETS2) {
00163
00164 if (*N_TobeAssigned) {
00165
00166 OffS = SUMA_Initialize_getoffsets (SO->N_Node);
00167 SUMA_getoffsets2 (dothisnode, SO, Opt->DistLim, OffS, NULL, 0);
00168 #if 0
00169 if (NewClust) {
00170 FILE *fid=NULL;
00171 char *s=NULL, tmp[50];
00172 fid = fopen("offsets2.1D", "w");
00173 if (!fid) {
00174 SUMA_SL_Err("Could not open file for writing.\nCheck file permissions, disk space.\n");
00175 } else {
00176 s = SUMA_ShowOffset_Info(OffS, 0);
00177 if (s) { fprintf(fid,"%s\n", s); SUMA_free(s); s = NULL;}
00178 fclose(fid);
00179 }
00180 }
00181 #endif
00182
00183 for (il=1; il<OffS->N_layers; ++il) {
00184 for (jl=0; jl<OffS->layers[il].N_NodesInLayer; ++jl) {
00185 neighb = OffS->layers[il].NodesInLayer[jl];
00186 if (ToBeAssigned[neighb] && OffS->OffVect[neighb] <= Opt->DistLim) {
00187
00188 SUMA_Build_Cluster_From_Node( neighb, Clust, ToBeAssigned, N_TobeAssigned, NodeArea, SO, Opt);
00189 }
00190 }
00191 }
00192
00193
00194 if (OffS) SUMA_Free_getoffsets(OffS); OffS = NULL;
00195 }
00196 } else if (BuildMethod == SUMA_OFFSETS_LL) {
00197 if (*N_TobeAssigned) {
00198
00199 if (!(offlist = SUMA_getoffsets_ll (dothisnode, SO, Opt->DistLim, NULL, 0))) {
00200 SUMA_SL_Err("Failed to get offsets.\nNo cleanup done.");
00201 SUMA_RETURN(NULL);
00202 }
00203 #if 0
00204 if (NewClust) {
00205 FILE *fid=NULL;
00206 char *s=NULL, tmp[50];
00207 fid = fopen("offsets_ll.1D", "w");
00208 if (!fid) {
00209 SUMA_SL_Err("Could not open file for writing.\nCheck file permissions, disk space.\n");
00210 } else {
00211 s = SUMA_ShowOffset_ll_Info(offlist, 0);
00212 if (s) { fprintf(fid,"%s\n", s); SUMA_free(s); s = NULL;}
00213 fclose(fid);
00214 }
00215 }
00216 #endif
00217
00218
00219 elm = dlist_head(offlist);
00220 dat = (SUMA_OFFSET_LL_DATUM *)elm->data;
00221 if (dat->layer != 0) {
00222 SUMA_SL_Err("Unexpected non zero layer for first element.");
00223 SUMA_RETURN(NULL);
00224 }
00225 do {
00226 elm = elm->next;
00227 dat = (SUMA_OFFSET_LL_DATUM *)elm->data;
00228 neighb = dat->ni;
00229 if (ToBeAssigned[neighb] && dat->off <= Opt->DistLim) {
00230
00231 SUMA_Build_Cluster_From_Node( neighb, Clust, ToBeAssigned, N_TobeAssigned, NodeArea, SO, Opt);
00232 }
00233 } while (elm != dlist_tail(offlist));
00234 dlist_destroy(offlist);
00235 }
00236 }
00237
00238
00239
00240 SUMA_RETURN(Clust);
00241 }
00242
00243
00244
00245
00246 #define SUMA_ADD_NODE_TO_CLUST(dothisnode, Clust, NodeArea, ToBeAssigned) { \
00247 if (LocalHead) fprintf(SUMA_STDERR,"%s: Adding node %d to cluster %p of %d nodes\n", FuncName, dothisnode, Clust, Clust->N_Node); \
00248 Clust->NodeList[Clust->N_Node] = dothisnode; \
00249 Clust->totalarea += NodeArea[dothisnode]; \
00250 Clust->totalvalue += ToBeAssigned[dothisnode]; \
00251 Clust->totalabsvalue += (float)fabs((float)ToBeAssigned[dothisnode]); \
00252 if (ToBeAssigned[dothisnode] < Clust->minvalue) { Clust->minvalue = ToBeAssigned[dothisnode]; Clust->minnode = dothisnode; } \
00253 if (ToBeAssigned[dothisnode] > Clust->maxvalue) { Clust->maxvalue = ToBeAssigned[dothisnode]; Clust->maxnode = dothisnode; } \
00254 Clust->ValueList[Clust->N_Node] = ToBeAssigned[dothisnode]; \
00255 ++Clust->N_Node; \
00256 }
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273 SUMA_CLUST_DATUM * SUMA_Build_Cluster_From_Node_NoRec ( int dothisnode,
00274 float *ToBeAssigned, int *N_TobeAssigned, float *NodeArea,
00275 SUMA_SurfaceObject *SO, SUMA_SURFCLUST_OPTIONS *Opt )
00276 {
00277 static char FuncName[]={"SUMA_Build_Cluster_From_Node_NoRec"};
00278 SUMA_CLUST_DATUM *Clust = NULL;
00279 static int ncall;
00280 static int N_Orig = -1;
00281 int il, jl, neighb, itmp;
00282 SUMA_GET_OFFSET_STRUCT *OffS = NULL;
00283 DList *offlist = NULL, *candlist=NULL;
00284 DListElmt *elm = NULL, *dothiselm=NULL;
00285 SUMA_OFFSET_LL_DATUM *dat=NULL;
00286 int NewClust = 0;
00287 byte *visited=NULL;
00288 SUMA_Boolean LocalHead = NOPE;
00289
00290 SUMA_ENTRY;
00291 ++ncall;
00292
00293 if (Opt->update < 0) { N_Orig = *N_TobeAssigned; Opt->update = -Opt->update; }
00294
00295 if (dothisnode < 0) {
00296 SUMA_SL_Err("Unexpected negative index.");
00297 SUMA_RETURN(NULL);
00298 }
00299
00300 OffS = SUMA_Initialize_getoffsets (SO->N_Node);
00301 Clust = (SUMA_CLUST_DATUM *)SUMA_malloc(sizeof(SUMA_CLUST_DATUM));
00302 Clust->N_Node = 0; Clust->totalarea = 0.0;
00303 Clust->totalvalue = 0.0; Clust->totalabsvalue = 0.0;
00304 Clust->minvalue = ToBeAssigned[dothisnode]; Clust->minnode = dothisnode;
00305 Clust->maxvalue = ToBeAssigned[dothisnode]; Clust->maxnode = dothisnode;
00306 Clust->varvalue = 0.0; Clust->centralnode = 0; Clust->weightedcentralnode = 0;
00307 Clust->NodeList = (int *)SUMA_malloc((*N_TobeAssigned) * sizeof(int));
00308 Clust->ValueList = (float *)SUMA_malloc((*N_TobeAssigned) * sizeof(float));
00309 if (!Clust->NodeList || !Clust->ValueList || !OffS) {
00310 SUMA_SL_Crit("Failed to allocate for NodeList or ValueList");
00311 SUMA_free(Clust); Clust = NULL;
00312 SUMA_RETURN(NULL);
00313 }
00314 candlist = (DList*)SUMA_malloc(sizeof(DList));
00315 visited = (byte *)SUMA_calloc(SO->N_Node, sizeof(byte));
00316 if (!visited || !candlist) {
00317 SUMA_SL_Crit("Failed to allocate for visited or candlist");
00318 SUMA_free(Clust); Clust = NULL;
00319 SUMA_RETURN(NULL);
00320 }
00321 dlist_init(candlist, NULL);
00322
00323 SUMA_ADD_NODE_TO_CLUST(dothisnode, Clust, NodeArea, ToBeAssigned);
00324
00325 ToBeAssigned[dothisnode] = 0; --(*N_TobeAssigned);
00326 visited[dothisnode] = YUP;
00327 dlist_ins_next(candlist, dlist_tail(candlist), (void *)dothisnode);
00328 while (*N_TobeAssigned && dlist_size(candlist)) {
00329
00330 dothiselm = dlist_head(candlist); dothisnode = (int) dothiselm->data;
00331 SUMA_getoffsets2 (dothisnode, SO, Opt->DistLim, OffS, NULL, 0);
00332
00333 dlist_remove(candlist, dothiselm, (void*)&itmp);
00334
00335 for (il=1; il<OffS->N_layers; ++il) {
00336 for (jl=0; jl<OffS->layers[il].N_NodesInLayer; ++jl) {
00337 neighb = OffS->layers[il].NodesInLayer[jl];
00338 if (ToBeAssigned[neighb] && OffS->OffVect[neighb] <= Opt->DistLim) {
00339
00340 SUMA_ADD_NODE_TO_CLUST(neighb, Clust, NodeArea, ToBeAssigned);
00341
00342 ToBeAssigned[neighb] = 0; --(*N_TobeAssigned);
00343 if (Opt->update) {
00344 if (N_Orig - *N_TobeAssigned >= Opt->update) {
00345 if (LocalHead) fprintf(SUMA_STDERR,"%s: tick (%d nodes processed)\n", FuncName, N_Orig - *N_TobeAssigned);
00346 else fprintf(SUMA_STDERR,".");
00347
00348 N_Orig = *N_TobeAssigned;
00349 }
00350 }
00351
00352 if (!visited[neighb]) {
00353 dlist_ins_next(candlist, dlist_tail(candlist), (void *)neighb);
00354 visited[neighb] = YUP;
00355 }
00356
00357 }
00358 }
00359 }
00360
00361 SUMA_Recycle_getoffsets (OffS);
00362 }
00363
00364 if (OffS) SUMA_Free_getoffsets(OffS); OffS = NULL;
00365 if (Opt->update) fprintf(SUMA_STDERR,"\n");
00366
00367 if (candlist) { dlist_destroy(candlist); SUMA_free(candlist); candlist = NULL; }
00368 SUMA_RETURN(Clust);
00369 }
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389 DList *SUMA_FindClusters ( SUMA_SurfaceObject *SO, int *ni, float *nv, int N_ni,
00390 int dothisnode, SUMA_SURFCLUST_OPTIONS *Opt,
00391 float *NodeArea)
00392 {
00393 static char FuncName[]={"SUMA_FindClusters"};
00394 DList *list=NULL;
00395 DListElmt *elm=NULL;
00396 float *ToBeAssigned=NULL;
00397 float mean;
00398 int N_n, nc, i, kk, PureNothing=0;
00399 SUMA_CLUST_DATUM *Clust = NULL;
00400 SUMA_Boolean LocalHead = NOPE;
00401
00402 SUMA_ENTRY;
00403
00404 if (!SO || !nv || !ni) {
00405 SUMA_S_Err("Bad parameters");
00406 }
00407 if (dothisnode == -1) {
00408 SUMA_LH("Initializing");
00409 nc = 0;
00410
00411 list = (DList *)SUMA_malloc(sizeof(DList));
00412 dlist_init(list, SUMA_FreeClustDatum);
00413 ToBeAssigned = (float *)SUMA_calloc(SO->N_Node, sizeof(float));
00414 N_n = N_ni;
00415 for (i=0; i<N_n; ++i) {
00416 if (!nv[i]) {
00417 ++PureNothing;
00418 }
00419 ToBeAssigned[ni[i]] = nv[i];
00420 }
00421 if (Opt->update) {
00422 fprintf(SUMA_STDERR,"%s: Have %d nodes to work with. %d nodes have 0 value.\n", FuncName, N_n, PureNothing);
00423 }
00424 }
00425
00426 while (N_n - PureNothing > 0) {
00427 dothisnode = -1;
00428 for (i=0;i<SO->N_Node; ++i) {
00429 if (ToBeAssigned[i]) {
00430 dothisnode = i; continue;
00431 }
00432 }
00433 if (dothisnode < 0) {
00434 SUMA_S_Err("Not expected here. dothisnode < 0");
00435 SUMA_RETURN(list);
00436 }
00437
00438 if (BuildMethod == SUMA_OFFSETS2_NO_REC) {
00439 Clust = SUMA_Build_Cluster_From_Node_NoRec(dothisnode, ToBeAssigned, &N_n, NodeArea, SO, Opt);
00440 } else if (BuildMethod == SUMA_OFFSETS2 || BuildMethod == SUMA_OFFSETS_LL) {
00441 Clust = SUMA_Build_Cluster_From_Node(dothisnode, NULL, ToBeAssigned, &N_n, NodeArea, SO, Opt);
00442 } else {
00443 SUMA_SL_Err("No Such Method!");
00444 SUMA_RETURN(list);
00445 }
00446 if (!Clust) {
00447 SUMA_SL_Err("Failed in SUMA_Build_Cluster_From_Node*");
00448 SUMA_RETURN(list);
00449 }
00450 if (LocalHead) fprintf(SUMA_STDERR,"%s: Cluster %p is finished, %d nodes\n", FuncName, Clust, Clust->N_Node);
00451
00452 if (Opt->AreaLim > 0 && Clust->totalarea < Opt->AreaLim) {
00453 SUMA_LH("Cluster less than area limit");
00454 SUMA_FreeClustDatum((void *)Clust); Clust = NULL;
00455 } else {
00456 mean = Clust->totalvalue/((float)Clust->N_Node);
00457 for (kk=0; kk < Clust->N_Node; ++kk) {
00458 Clust->varvalue += (Clust->ValueList[kk] - mean) * (Clust->ValueList[kk] - mean);
00459 }
00460 if (Clust->N_Node > 1) Clust->varvalue /= (Clust->N_Node - 1);
00461 else Clust->varvalue = 0.0;
00462
00463 Clust->NodeList = (int *)SUMA_realloc(Clust->NodeList, sizeof(int)*Clust->N_Node);
00464 Clust->ValueList = (float *)SUMA_realloc(Clust->ValueList, sizeof(float)*Clust->N_Node);
00465 if (!Clust->NodeList || !Clust->ValueList) {
00466 SUMA_SL_Crit("Failed to reallocate for NodeList or ValueList");
00467 SUMA_RETURN(NULL);
00468 }
00469
00470 if (Opt->DoCentrality) {
00471 if (Opt->update) {
00472 SUMA_SL_Note("Looking for central nodes...(use -no_cent to skip this slow step)");
00473 }
00474 if (!SUMA_ClusterCenterofMass (SO, Clust, 1)) {
00475 SUMA_SL_Err("Failed to find central node");
00476 SUMA_RETURN(list);
00477 }
00478 }
00479
00480 dlist_ins_next(list, dlist_tail(list), (void *)Clust);
00481 ++nc;
00482 }
00483 }
00484
00485 if (N_n == 0) {
00486 if (LocalHead) fprintf(SUMA_STDERR,"%s: No more nodes to consider, cleaning up.\n", FuncName);
00487 if (ToBeAssigned) SUMA_free(ToBeAssigned); ToBeAssigned = NULL; \
00488 }
00489
00490 SUMA_RETURN(list);
00491 }
00492
00493
00494 SUMA_Boolean SUMA_Show_SurfClust_list(DList *list, FILE *Out, int detail, char *params)
00495 {
00496 static char FuncName[]={"SUMA_Show_SurfClust_list"};
00497 char *s = NULL;
00498
00499 SUMA_ENTRY;
00500
00501 if (Out == NULL) Out = stdout;
00502
00503 s = SUMA_Show_SurfClust_list_Info(list, detail, params);
00504 if (!s) {
00505 SUMA_SL_Err("Failed in SUMA_Show_SurfClust_list_Info");
00506 SUMA_RETURN(NOPE);
00507 } else {
00508 fprintf(Out, "%s", s);
00509 SUMA_free(s); s = NULL;
00510 }
00511
00512 SUMA_RETURN(YUP);
00513 }
00514
00515
00516 char *SUMA_Show_SurfClust_list_Info(DList *list, int detail, char *params)
00517 {
00518 static char FuncName[]={"SUMA_Show_SurfClust_list_Info"};
00519 int i, ic, max;
00520 SUMA_STRING *SS = NULL;
00521 DListElmt *elmt=NULL;
00522 SUMA_CLUST_DATUM *cd=NULL;
00523 char *s=NULL, *pad_str, str[20];
00524 int lc[]= { 6, 6, 9, 9, 9, 6, 6, 9, 6, 9, 6, 9, 9 };
00525 char Col[][12] = { {"# Rank"}, {"num Nd"}, {"Area"}, {"Mean"}, {"|Mean|"},{"Cent"}, {"W Cent"},{"Min V"}, {"Min Nd"}, {"Max V"}, {"Max Nd"} , {"Var"}, {"SEM"} };
00526 SUMA_ENTRY;
00527
00528 SS = SUMA_StringAppend (NULL, NULL);
00529
00530 if (!list) {
00531 SS = SUMA_StringAppend (SS,"NULL list.\n");
00532 SUMA_SS2S(SS,s);
00533 SUMA_RETURN(s);
00534 }
00535
00536 if (!list->size) {
00537 SS = SUMA_StringAppend (SS,"Empty list.\n");
00538 SUMA_SS2S(SS,s);
00539 SUMA_RETURN(s);
00540 }else{
00541 SS = SUMA_StringAppend_va (SS,"#Col. 0 = Rank \n"
00542 "#Col. 1 = Number of nodes\n"
00543 "#Col. 2 = Total area (units^2)\n"
00544 "#Col. 3 = Mean value\n"
00545 "#Col. 4 = Mean absolute value\n"
00546 "#Col. 5 = Central node\n"
00547 "#Col. 6 = Weighted central node\n"
00548 "#Col. 7 = Minimum value\n"
00549 "#Col. 8 = Minimum node\n"
00550 "#Col. 9 = Maximum value\n"
00551 "#Col. 10 = Maximum node\n"
00552 "#Col. 11 = Variance\n"
00553 "#Col. 12 = Standard error of the mean\n"
00554 "#Command history:\n"
00555 "#%s\n", params);
00556
00557 for (ic=0; ic<13; ++ic) {
00558 if (ic == 0) sprintf(str,"%s", Col[ic]);
00559 else sprintf(str,"%s", Col[ic]);
00560 pad_str = SUMA_pad_string(str, ' ', lc[ic], 0);
00561 SS = SUMA_StringAppend_va (SS,"%s ", pad_str);
00562 SUMA_free(pad_str);
00563 }
00564 SS = SUMA_StringAppend_va (SS,"\n");
00565 if (detail == 1) {
00566 SS = SUMA_StringAppend_va (SS,"#Other columns: list of 5 first nodes in ROI.\n");
00567 }
00568 if (detail == 2) {
00569 SS = SUMA_StringAppend_va (SS,"#Other columns: list all nodes in ROI.\n");
00570 }
00571 if (detail > 0) {
00572 SS = SUMA_StringAppend_va (SS,"#A total of %d cluster%s were found.\n", list->size, SUMA_COUNTER_PLURAL(list->size));
00573 }
00574 }
00575
00576 elmt = NULL;
00577 ic = 1;
00578 do {
00579 if (!elmt) elmt = dlist_head(list); else elmt = elmt->next;
00580 if (!elmt) SS = SUMA_StringAppend_va (SS,"#%d%s cluster element is NULL!\n", ic, SUMA_COUNTER_SUFFIX(ic));
00581 else {
00582 cd = (SUMA_CLUST_DATUM *)elmt->data;
00583 if (detail > 0) SS = SUMA_StringAppend_va (SS,"#%d%s cluster\n", ic, SUMA_COUNTER_SUFFIX(ic));
00584 SS = SUMA_StringAppend_va (SS,"%6d %6d %9.2f"
00585 " %9.3f %9.3f"
00586 " %6d %6d"
00587 " %9.3f %6d"
00588 " %9.3f %6d"
00589 " %9.3f %9.3f"
00590 , ic, cd->N_Node, cd->totalarea
00591 , cd->totalvalue/((float)cd->N_Node), cd->totalabsvalue/((float)cd->N_Node)
00592 , cd->centralnode, cd->weightedcentralnode
00593 , cd->minvalue, cd->minnode
00594 , cd->maxvalue, cd->maxnode
00595 , cd->varvalue, sqrt(cd->varvalue/cd->N_Node));
00596 if (detail > 0) {
00597 if (detail == 1) {
00598 if (cd->N_Node < 5) max = cd->N_Node; else max = 5;
00599 } else max = cd->N_Node;
00600 for (i=0;i<max; ++i) SS = SUMA_StringAppend_va (SS,"%d\t", cd->NodeList[i]);
00601 }
00602 SS = SUMA_StringAppend(SS,"\n");
00603 }
00604 ++ic;
00605 } while (elmt != dlist_tail(list));
00606
00607 SUMA_SS2S(SS,s);
00608
00609 SUMA_RETURN (s);
00610 }
00611
00612
00613 SUMA_DSET *SUMA_MaskDsetByClustList(SUMA_DSET *idset, SUMA_SurfaceObject *SO, DList *list, SUMA_Boolean FullList, char *leName)
00614 {
00615 static char FuncName[]={"SUMA_MaskDsetByClustList"};
00616 int i, j;
00617 DListElmt *elmt=NULL;
00618 SUMA_DSET *dset = NULL;
00619 int *ni=NULL, N_ni, cnt;
00620 SUMA_CLUST_DATUM *cd=NULL;
00621 byte *ismask=NULL, *rowmask=NULL, *colmask = NULL;
00622 SUMA_Boolean LocalHead = NOPE;
00623
00624 SUMA_ENTRY;
00625
00626 if (!list) {
00627 SUMA_SL_Err("NULL list");
00628 SUMA_RETURN(dset);
00629 }
00630
00631 ismask = (byte *)SUMA_calloc(SO->N_Node, sizeof(byte));
00632
00633 elmt = NULL; cnt = 0;
00634 do {
00635 if (!elmt) elmt = dlist_head(list);
00636 else elmt = elmt->next;
00637 cd = (SUMA_CLUST_DATUM *)elmt->data;
00638 for (j=0; j<cd->N_Node; ++j) {
00639 if(LocalHead) fprintf (SUMA_STDERR,"nic=%d\t", cd->NodeList[j]);
00640 ismask[cd->NodeList[j]] = 1;
00641 ++cnt;
00642 }
00643 } while (elmt != dlist_tail(list));
00644 if (LocalHead) fprintf(SUMA_STDERR,"%s:\n%d nodes in cluster list.\n", FuncName, cnt);
00645
00646
00647 rowmask = (byte *)SUMA_calloc(idset->dnel->vec_len, sizeof(byte));
00648 colmask = (byte *)SUMA_calloc(idset->dnel->vec_num, sizeof(byte));
00649
00650
00651 ni = SUMA_GetNodeDef(idset);
00652 N_ni = SDEST_VECLEN(idset);
00653 if (!ni) {
00654 SUMA_SL_Err("Failed to find node index column");
00655 SUMA_RETURN(NULL);
00656 }
00657
00658 for (i=0; i<idset->dnel->vec_len; ++i) { if (ismask[ni[i]]) { rowmask[i]=1; if(LocalHead) fprintf (SUMA_STDERR,"%d,%d\t", ni[i], i); } }
00659
00660 for (i=0; i<idset->dnel->vec_num; ++i) {
00661 if (SUMA_isDsetColumn_inferred(idset, i)) {
00662 colmask[i]=0;
00663 if (LocalHead) fprintf(SUMA_STDERR,"%s: Column %d will not be written because it is inferred.\n", FuncName, i);
00664 } else colmask[i]=1;
00665 }
00666
00667 dset = SUMA_MaskedCopyofDset(idset, rowmask, colmask, !FullList, 1);
00668 if (!dset) {
00669 SUMA_SL_Err("Failed to create masked copy of input");
00670 SUMA_RETURN(NULL);
00671 }
00672
00673 if (rowmask) SUMA_free(rowmask); rowmask = NULL;
00674 if (colmask) SUMA_free(colmask); colmask = NULL;
00675 if (ismask) SUMA_free(ismask); ismask = NULL;
00676
00677 SUMA_RETURN (dset);
00678 }
00679
00680 SUMA_DSET *SUMA_SurfClust_list_2_DsetMask(SUMA_SurfaceObject *SO, DList *list, SUMA_Boolean FullList, char *leName)
00681 {
00682 static char FuncName[]={"SUMA_SurfClust_list_2_DsetMask"};
00683 int i, ic, max, j, rank;
00684 DListElmt *elmt=NULL;
00685 SUMA_DSET *dset = NULL;
00686 int *NodeIndex=NULL, N_Node, *Val = NULL;
00687 SUMA_CLUST_DATUM *cd=NULL;
00688 SUMA_Boolean LocalHead = NOPE;
00689
00690 SUMA_ENTRY;
00691
00692 if (!list) {
00693 SUMA_SL_Err("NULL list");
00694 SUMA_RETURN(dset);
00695 }
00696 if (FullList) N_Node = SO->N_Node;
00697 else {
00698 elmt = NULL; N_Node = 0;
00699 do {
00700 if (!elmt) elmt = dlist_head(list);
00701 else elmt = elmt->next;
00702 cd = (SUMA_CLUST_DATUM *)elmt->data;
00703 N_Node += cd->N_Node;
00704 } while (elmt != dlist_tail(list));
00705 }
00706 NodeIndex = (int *)SUMA_malloc(N_Node*sizeof(int));
00707 Val = (int *)SUMA_malloc(N_Node * sizeof(int));
00708 if (!NodeIndex || !Val){
00709 SUMA_SL_Crit("Failed to allocate NodeIndex and or Val");
00710 SUMA_RETURN(dset);
00711 }
00712 if (FullList) {
00713 for (i=0; i<N_Node; ++i) {
00714 NodeIndex[i] = i; Val[i] = 0;
00715 }
00716 elmt = NULL; rank = 1;
00717 do {
00718 if (!elmt) elmt = dlist_head(list);
00719 else elmt = elmt->next;
00720 cd = (SUMA_CLUST_DATUM *)elmt->data;
00721 for (j=0; j<cd->N_Node; ++j) {
00722 Val[cd->NodeList[j]] = rank;
00723 }
00724 ++rank;
00725 } while (elmt != dlist_tail(list));
00726 } else {
00727 elmt = NULL; rank = 1; i = 0;
00728 do {
00729 if (!elmt) elmt = dlist_head(list);
00730 else elmt = elmt->next;
00731 cd = (SUMA_CLUST_DATUM *)elmt->data;
00732 for (j=0; j<cd->N_Node; ++j) {
00733 Val[i] = rank;
00734 NodeIndex[i] = cd->NodeList[j];
00735 ++i;
00736 }
00737 ++rank;
00738 } while (elmt != dlist_tail(list));
00739 }
00740
00741 SUMA_LH("Creating dset pointer");
00742 dset = SUMA_CreateDsetPointer(
00743 leName,
00744 SUMA_NODE_ROI,
00745 NULL,
00746 NULL,
00747 N_Node
00748 );
00749
00750
00751 SUMA_LH("Adding NodeDef column ...");
00752 if (!SUMA_AddDsetNelCol ( dset,
00753 "le Node Def",
00754 SUMA_NODE_INDEX,
00755 (void *)NodeIndex,
00756 NULL,
00757 1
00758 )) {
00759 fprintf (stderr,"Error %s:\nFailed in SUMA_AddNelCol", FuncName);
00760 SUMA_RETURN(NULL);
00761 }
00762
00763 if (!SUMA_AddDsetNelCol (dset, "Cluster Rank", SUMA_NODE_INT, (void *)Val, NULL ,1)) {
00764 fprintf (stderr,"Error %s:\nFailed in SUMA_AddNelCol", FuncName);
00765 SUMA_RETURN (NULL);
00766 }
00767
00768 if (NodeIndex) SUMA_free(NodeIndex); NodeIndex = NULL;
00769 if (Val) SUMA_free(Val); Val = NULL;
00770
00771 SUMA_RETURN (dset);
00772 }
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789 int SUMA_ClusterCenterofMass (SUMA_SurfaceObject *SO, SUMA_CLUST_DATUM *cd, int UseSurfDist)
00790 {
00791 static char FuncName[]={"SUMA_ClusterCenterofMass"};
00792 SUMA_GET_OFFSET_STRUCT *OffS = NULL;
00793 SUMA_OFFSET_STRUCT OS;
00794 int *CoverThisNode = NULL, i, c, ni, nc, centralnode, weightedcentralnode,
00795 nanch, k, WeightByValue = 1;
00796 float Uia[3], dia, Uca[3], dca, s[3], s_w[3], *DistVec=NULL,
00797 *weightedcentrality=NULL, minweightedcentrality = 0.0, *centrality=NULL,
00798 mincentrality = 0.0, mtotal_w, mi_w, fac, mtotal, mi;
00799 static int ncall;
00800 SUMA_Boolean LocalHead = NOPE;
00801
00802 SUMA_ENTRY;
00803
00804 centralnode = -1;
00805 weightedcentralnode = -1;
00806 if (!SO || !cd) { SUMA_RETURN(NOPE); }
00807 if (!cd->N_Node || !cd->NodeList) { SUMA_RETURN(NOPE); }
00808 OffS = SUMA_Initialize_getoffsets (SO->N_Node);
00809 if (!OffS) {
00810 SUMA_SL_Err("Failed to allocate for OffS");
00811 SUMA_RETURN(NOPE);
00812 }
00813
00814 CoverThisNode = (int *)SUMA_calloc(SO->N_Node, sizeof(int));
00815 if (!CoverThisNode) {
00816 SUMA_SL_Crit("Failed to allocate for CoverThisNode");
00817 SUMA_RETURN(NOPE);
00818 }
00819 if (cd->N_Node == SO->N_Node) {
00820
00821 SUMA_SL_Note("Cluster spans entire surface.\nNo central node.\n");
00822 cd->centralnode = 0;
00823 cd->weightedcentralnode = 0;
00824 SUMA_RETURN(YUP);
00825 }
00826
00827 for (i=0; i<cd->N_Node; ++i) { CoverThisNode[cd->NodeList[i]] = 1; }
00828 nanch = cd->NodeList[0];
00829 SUMA_getoffsets2 (cd->NodeList[0], SO, 0, OffS, CoverThisNode, cd->N_Node);
00830
00831 if (!SUMA_GetOffset2Offset (OffS, &OS)) {
00832 SUMA_SL_Err("Failed in SUMA_GetOffset2Offset");
00833 SUMA_RETURN(NOPE);
00834 }
00835 if (OffS) SUMA_Free_getoffsets(OffS); OffS = NULL;
00836
00837 DistVec = (float *) SUMA_calloc(SO->N_Node, sizeof(float));
00838 centrality = (float *)SUMA_malloc(OS.N_Neighb * sizeof(float));
00839 weightedcentrality = (float *)SUMA_malloc(OS.N_Neighb * sizeof(float));
00840 if (!DistVec || !centrality || !weightedcentrality) {
00841 SUMA_SL_Crit("Failed to allocate.");
00842 SUMA_RETURN(NOPE);
00843 }
00844 for (c=0; c < OS.N_Neighb; ++c) { DistVec[OS.Neighb_ind[c]] = OS.Neighb_dist[c]; }
00845
00846
00847
00848
00849 #if 1
00850 if (cd->N_Node == 1) {
00851 centralnode = cd->NodeList[0];
00852 weightedcentralnode = cd->NodeList[0];
00853 } else {
00854 for (c=0; c < OS.N_Neighb; ++c) {
00855 nc = OS.Neighb_ind[c];
00856 s[0] = s[1] = s[2] = 0.0; s_w[0] = s_w[1] = s_w[2] = 0.0;
00857 centrality[c] = 0.0; weightedcentrality[c] = 0.0; mtotal_w = 0.0;
00858 for (i=0; i<cd->N_Node; ++i) {
00859 mi_w = cd->ValueList[i];
00860 mtotal_w += mi_w;
00861 ni = cd->NodeList[i];
00862 SUMA_UNIT_VEC((&(SO->NodeList[3*ni])), (&(SO->NodeList[3*nanch])), Uia, dia);
00863 if (UseSurfDist) { for (k=0; k<3;++k) { fac = Uia[k] * DistVec[ni]; s_w[k] += mi_w * fac; s[k] += fac;} }
00864 else { for (k=0; k<3;++k) { fac = Uia[k] * dia; s_w[k] += mi_w * fac; s[k] += fac; } }
00865 }
00866 SUMA_UNIT_VEC((&(SO->NodeList[3*nc])), (&(SO->NodeList[3*nanch])), Uca, dca);
00867 if (UseSurfDist) { for (k=0; k<3;++k) { fac = Uca[k] * DistVec[nc]; s_w[k] -= mtotal_w * fac; s[k] -= cd->N_Node * fac; } }
00868 else { for (k=0; k<3;++k) { fac = Uca[k] * dca; s_w[k] -= mtotal_w * fac; s[k] -= cd->N_Node * fac;} }
00869
00870 SUMA_NORM_VEC(s, 3, centrality[c]); SUMA_NORM_VEC(s_w, 3, weightedcentrality[c]);
00871 if (LocalHead) fprintf(SUMA_STDERR,"%s: call %d, Centrality of node %d / %d = %f , weighted %f\n", FuncName, ncall, nc, nanch, centrality[c], weightedcentrality[c]);
00872 if (c == 0) {
00873 mincentrality = centrality[c]; centralnode = nc;
00874 minweightedcentrality = weightedcentrality[c]; weightedcentralnode = nc;
00875 } else {
00876 if (centrality[c] < mincentrality) { mincentrality = centrality[c]; centralnode = nc; }
00877 if (weightedcentrality[c] < minweightedcentrality) { minweightedcentrality = weightedcentrality[c]; weightedcentralnode = nc; }
00878 }
00879
00880 }
00881 }
00882 cd->centralnode = centralnode;
00883 cd->weightedcentralnode = weightedcentralnode;
00884 if (LocalHead) fprintf(SUMA_STDERR,"%s: Central node chosen %d, weighted %d\n", FuncName, cd->centralnode, cd->weightedcentralnode);
00885 #else
00886 if (cd->N_Node == 1) {
00887 centralnode = cd->NodeList[0];
00888 } else {
00889 for (c=0; c < OS.N_Neighb; ++c) {
00890 nc = OS.Neighb_ind[c];
00891 s[0] = s[1] = s[2] = 0.0;
00892 centrality[c] = 0.0; mtotal = 0.0;
00893 for (i=0; i<cd->N_Node; ++i) {
00894 if (WeightByValue) mi = cd->ValueList[i];
00895 else mi = 1.0;
00896 mtotal += mi;
00897 ni = cd->NodeList[i];
00898 SUMA_UNIT_VEC((&(SO->NodeList[3*ni])), (&(SO->NodeList[3*nanch])), Uia, dia);
00899 if (UseSurfDist) { for (k=0; k<3;++k) s[k] += mi * Uia[k] * DistVec[ni]; }
00900 else { for (k=0; k<3;++k) s[k] += mi * Uia[k] * dia; }
00901 }
00902 SUMA_UNIT_VEC((&(SO->NodeList[3*nc])), (&(SO->NodeList[3*nanch])), Uca, dca);
00903 if (UseSurfDist) { for (k=0; k<3;++k) s[k] -= mtotal * Uca[k] * DistVec[nc]; }
00904 else { for (k=0; k<3;++k) s[k] -= mtotal * Uca[k] * dca; }
00905
00906 SUMA_NORM_VEC(s, 3, centrality[c]);
00907 if (LocalHead) fprintf(SUMA_STDERR,"%s: call %d, Centrality of node %d / %d = %f\n", FuncName, ncall, nc, nanch, centrality[c]);
00908 if (c == 0) { mincentrality = centrality[c]; centralnode = nc; }
00909 else if (centrality[c] < mincentrality) { mincentrality = centrality[c]; centralnode = nc; }
00910 }
00911 }
00912 cd->centralnode = centralnode;
00913 if (LocalHead) fprintf(SUMA_STDERR,"%s: Central node chosen %d\n", FuncName, cd->centralnode);
00914 cd->weightedcentralnode = -1;
00915 #endif
00916
00917 if (CoverThisNode) SUMA_free(CoverThisNode); CoverThisNode = NULL;
00918 if (OS.Neighb_dist) SUMA_free(OS.Neighb_dist); OS.Neighb_dist = NULL;
00919 if (OS.Neighb_ind) SUMA_free(OS.Neighb_ind); OS.Neighb_ind = NULL;
00920 if (DistVec) SUMA_free(DistVec); DistVec = NULL;
00921 if (centrality) SUMA_free(centrality); centrality = NULL;
00922 if (weightedcentrality) SUMA_free(weightedcentrality); weightedcentrality = NULL;
00923
00924 ++ncall;
00925 SUMA_RETURN(YUP);
00926 }
00927
00928 SUMA_Boolean SUMA_Sort_ClustersList (DList *list, SUMA_SURF_CLUST_SORT_MODES SortMode)
00929 {
00930 static char FuncName[]={"SUMA_Sort_ClustersList"};
00931 DListElmt *elmt=NULL, *max_elmt=NULL, *elmt_comp=NULL;
00932 SUMA_CLUST_DATUM *cd=NULL, *cd_comp=NULL, *cd_max=NULL;
00933 int r = 0;
00934 SUMA_Boolean LocalHead = NOPE;
00935
00936 SUMA_ENTRY;
00937
00938 if (!list) { SUMA_RETURN(NOPE); }
00939 switch (SortMode) {
00940 case SUMA_SORT_CLUST_NOT_SET:
00941 SUMA_S_Err("Why are you calling me?");
00942 SUMA_RETURN(NOPE);
00943 break;
00944 case SUMA_SORT_CLUST_NO_SORT:
00945 SUMA_RETURN(YUP);
00946 break;
00947 case SUMA_SORT_CLUST_BY_NUMBER_NODES:
00948 case SUMA_SORT_CLUST_BY_AREA:
00949 elmt = NULL;
00950 do {
00951 if (!elmt) elmt = dlist_head(list);
00952 else elmt = elmt->next;
00953 cd = (SUMA_CLUST_DATUM *)elmt->data;
00954
00955 if (elmt != dlist_tail(list)) {
00956 max_elmt = elmt; cd_max = (SUMA_CLUST_DATUM *)max_elmt->data; elmt_comp = NULL;
00957 do {
00958 if (!elmt_comp) elmt_comp = elmt->next;
00959 else elmt_comp = elmt_comp->next;
00960 cd_comp = (SUMA_CLUST_DATUM *)elmt_comp->data;
00961 if (SortMode == SUMA_SORT_CLUST_BY_NUMBER_NODES) {
00962 if (cd_comp->N_Node > cd_max->N_Node) { max_elmt = elmt_comp; cd_max = (SUMA_CLUST_DATUM *)max_elmt->data; }
00963 }else if (SortMode == SUMA_SORT_CLUST_BY_AREA) {
00964 if (cd_comp->totalarea > cd_max->totalarea) { max_elmt = elmt_comp; cd_max = (SUMA_CLUST_DATUM *)max_elmt->data; }
00965 }
00966 } while (elmt_comp != dlist_tail(list));
00967 if (max_elmt != elmt) {
00968 dlist_remove(list, max_elmt, (void *)(&cd_max));
00969 dlist_ins_prev(list, elmt, (void *)cd_max);
00970 elmt = elmt->prev;
00971 }
00972 }
00973 } while (elmt != dlist_tail(list));
00974 SUMA_RETURN(YUP);
00975 break;
00976 default:
00977 break;
00978 }
00979
00980
00981
00982 SUMA_RETURN(YUP);
00983 }
00984 #ifdef SUMA_SurfClust_STANDALONE
00985 void usage_SUMA_SurfClust ()
00986 {
00987 static char FuncName[]={"usage_SUMA_SurfClust"};
00988 char * s = NULL;
00989 s = SUMA_help_basics();
00990 printf ( "\nUsage: A program to perform clustering analysis surfaces.\n"
00991 " SurfClust <-spec SpecFile> \n"
00992 " <-surf_A insurf> \n"
00993 " <-input inData.1D dcol_index> \n"
00994 " <-rmm rad>\n"
00995 " [-amm2 minarea]\n"
00996 " [-prefix OUTPREF] \n"
00997 " [-out_clusterdset] [-out_roidset] \n"
00998 " [-out_fulllist]\n"
00999 " [-sort_none | -sort_n_nodes | -sort_area]\n"
01000 "\n"
01001 " The program can outputs a table of the clusters on the surface,\n"
01002 " a mask dataset formed by the different clusters and a clustered\n"
01003 " version of the input dataset.\n"
01004 "\n"
01005 " Mandatory parameters:\n"
01006 " -spec SpecFile: The surface spec file.\n"
01007 " -surf_A insurf: The input surface name.\n"
01008 " -input inData.1D dcol_index: The input 1D dataset\n"
01009 " and the index of the\n"
01010 " datacolumn to use\n"
01011 " (index 0 for 1st column).\n"
01012 " Values of 0 indicate \n"
01013 " inactive nodes.\n"
01014 " -rmm rad: Maximum distance between an activated node\n"
01015 " and the cluster to which it belongs.\n"
01016 " Distance is measured on the surface's graph (mesh).\n"
01017 "\n"
01018 " Optional Parameters:\n"
01019 " -thresh_col tcolind: Index of thresholding column.\n"
01020 " Default is column 0.\n "
01021 " -thresh tval: Apply thresholding prior to clustering.\n"
01022 " A node n is considered if thresh_col[n] > tval.\n"
01023 " -athresh tval: Apply absolute thresholding prior to clustering.\n"
01024 " A node n is considered if | thresh_col[n] | > tval.\n"
01025 " -amm2 minarea: Do not output resutls for clusters having\n"
01026 " an area less than minarea.\n"
01027 " -prefix OUTPREF: Prefix for output.\n"
01028 " Default is the prefix of \n"
01029 " the input dataset.\n"
01030 " If this option is used, the\n"
01031 " cluster table is written to a file called\n"
01032 " OUTPREF_ClstTable_rXX_aXX.1D. Otherwise the\n"
01033 " table is written to stdout. \n"
01034 " -out_clusterdset: Output a clustered version of inData.1D \n"
01035 " preserving only the values of nodes that \n"
01036 " belong to clusters that passed the rmm and amm2\n"
01037 " conditions above.\n"
01038 " The clustered dset's prefix has\n"
01039 " _Clustered_rXX_aXX affixed to the OUTPREF\n"
01040 " -out_roidset: Output an ROI dataset with the value\n"
01041 " at each node being the rank of its\n"
01042 " cluster. The ROI dataset's prefix has\n"
01043 " _ClstMsk_rXX_aXX affixed to the OUTPREF\n"
01044 " where XX represent the values for the\n"
01045 " the -rmm and -amm2 options respectively.\n"
01046 " The program will not overwrite pre-existing\n"
01047 " dsets.\n"
01048 " -out_fulllist: Output a value for all nodes of insurf.\n"
01049 " This option must be used in conjuction with\n"
01050 " -out_roidset and/or out_clusterdset.\n"
01051 " With this option, the output files might\n"
01052 " be mostly 0, if you have small clusters.\n"
01053 " However, you should use it if you are to \n"
01054 " maintain the same row-to-node correspondence\n"
01055 " across multiple datasets.\n"
01056 " -sort_none: No sorting of ROI clusters.\n"
01057 " -sort_n_nodes: Sorting based on number of nodes\n"
01058 " in cluster.\n"
01059 " -sort_area: Sorting based on area of clusters \n"
01060 " (default).\n"
01061 " -update perc: Pacify me when perc of the data have been\n"
01062 " processed. perc is between 1%% and 50%%.\n"
01063 " Default is no update.\n"
01064 " -no_cent: Do not find the central nodes.\n"
01065 " Finding the central node is a \n"
01066 " relatively slow operation. Use\n"
01067 " this option to skip it.\n"
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087 "\n"
01088 " The cluster table output:\n"
01089 " A table where ach row shows results from one cluster.\n"
01090 " Each row contains 13 columns: \n"
01091 " Col. 0 Rank of cluster (sorting order).\n"
01092 " Col. 1 Number of nodes in cluster.\n"
01093 " Col. 2 Total area of cluster. Units are the \n"
01094 " the surface coordinates' units^2.\n"
01095 " Col. 3 Mean data value in cluster.\n"
01096 " Col. 4 Mean of absolute data value in cluster.\n"
01097 " Col. 5 Central node of cluster (see below).\n"
01098 " Col. 6 Weighted central node (see below).\n"
01099 " Col. 7 Minimum value in cluster.\n"
01100 " Col. 8 Node where minimum value occurred.\n"
01101 " Col. 9 Maximum value in cluster.\n"
01102 " Col. 10 Node where maximum value occurred.\n"
01103 " Col. 11 Variance of values in cluster.\n"
01104 " Col. 12 Standard error of the mean ( sqrt(variance / number of nodes) ).\n"
01105 " The CenterNode n is such that: \n"
01106 " ( sum (Uia * dia * wi) ) - ( Uca * dca * sum (wi) ) is minimal\n"
01107 " where i is a node in the cluster\n"
01108 " a is an anchor node on the surface\n"
01109 " sum is carried over all nodes i in a cluster\n"
01110 " w. is the weight of a node \n"
01111 " = 1.0 for central node \n"
01112 " = value at node for the weighted central node\n"
01113 " U.. is the unit vector between two nodes\n"
01114 " d.. is the distance between two nodes on the graph\n"
01115 " (an approximation of the geodesic distance)\n"
01116 " If -no_cent is used, CenterNode columns are set to 0.\n"
01117 "\n"
01118 "%s"
01119 "\n", s);
01120
01121
01122
01123
01124 SUMA_free(s); s = NULL;
01125 s = SUMA_New_Additions(0, 1); printf("%s\n", s);SUMA_free(s); s = NULL;
01126 printf(" Ziad S. Saad SSCC/NIMH/NIH ziad@nih.gov \n");
01127 exit (0);
01128 }
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140 SUMA_SURFCLUST_OPTIONS *SUMA_SurfClust_ParseInput (char *argv[], int argc)
01141 {
01142 static char FuncName[]={"SUMA_SurfClust_ParseInput"};
01143 SUMA_SURFCLUST_OPTIONS *Opt=NULL;
01144 int kar, i, ind;
01145 char *outname;
01146 SUMA_Boolean brk = NOPE;
01147 SUMA_Boolean LocalHead = NOPE;
01148
01149 SUMA_ENTRY;
01150
01151 Opt = (SUMA_SURFCLUST_OPTIONS *)SUMA_malloc(sizeof(SUMA_SURFCLUST_OPTIONS));
01152
01153 kar = 1;
01154 Opt->spec_file = NULL;
01155 Opt->out_prefix = NULL;
01156 Opt->sv_name = NULL;
01157 Opt->N_surf = -1;
01158 Opt->DistLim = -1.0;
01159 Opt->AreaLim = -1.0;
01160 Opt->in_name = NULL;
01161 Opt->nodecol = -1;
01162 Opt->labelcol = -1;
01163 Opt->OutROI = NOPE;
01164 Opt->OutClustDset = NOPE;
01165 Opt->FullROIList = NOPE;
01166 Opt->WriteFile = NOPE;
01167 Opt->DoThreshold = 0;
01168 Opt->Thresh = 0.0;
01169 Opt->tind = 0;
01170 Opt->update = 0;
01171 Opt->SortMode = SUMA_SORT_CLUST_NOT_SET;
01172 Opt->DoCentrality = 1;
01173 for (i=0; i<SURFCLUST_MAX_SURF; ++i) { Opt->surf_names[i] = NULL; }
01174 outname = NULL;
01175 BuildMethod = SUMA_OFFSETS2_NO_REC;
01176 brk = NOPE;
01177 while (kar < argc) {
01178
01179 if (strcmp(argv[kar], "-h") == 0 || strcmp(argv[kar], "-help") == 0) {
01180 usage_SUMA_SurfClust();
01181 exit (0);
01182 }
01183
01184 SUMA_SKIP_COMMON_OPTIONS(brk, kar);
01185
01186 if (!brk && (strcmp(argv[kar], "-no_cent") == 0)) {
01187 Opt->DoCentrality = 0;
01188 brk = YUP;
01189 }
01190
01191 if (!brk && (strcmp(argv[kar], "-O2") == 0)) {
01192 BuildMethod = SUMA_OFFSETS2;
01193 brk = YUP;
01194 }
01195 if (!brk && (strcmp(argv[kar], "-O2_NR") == 0)) {
01196 BuildMethod = SUMA_OFFSETS2_NO_REC;
01197 brk = YUP;
01198 }
01199
01200 if (!brk && (strcmp(argv[kar], "-Oll") == 0)) {
01201 BuildMethod = SUMA_OFFSETS_LL;
01202 brk = YUP;
01203 }
01204
01205 if (!brk && (strcmp(argv[kar], "-spec") == 0)) {
01206 kar ++;
01207 if (kar >= argc) {
01208 fprintf (SUMA_STDERR, "need argument after -spec \n");
01209 exit (1);
01210 }
01211 Opt->spec_file = argv[kar];
01212 brk = YUP;
01213 }
01214
01215 if (!brk && (strcmp(argv[kar], "-sv") == 0)) {
01216 kar ++;
01217 if (kar >= argc) {
01218 fprintf (SUMA_STDERR, "need argument after -sv \n");
01219 exit (1);
01220 }
01221 Opt->sv_name = argv[kar];
01222 brk = YUP;
01223 }
01224
01225
01226 if (!brk && (strcmp(argv[kar], "-update") == 0)) {
01227 kar ++;
01228 if (kar >= argc) {
01229 fprintf (SUMA_STDERR, "need argument after -update \n");
01230 exit (1);
01231 }
01232 Opt->update = atof(argv[kar]);
01233 if (Opt->update < 1 || Opt->update > 100) {
01234 fprintf (SUMA_STDERR, "-update needs a parameter between 1 and 50 (I have %.1f)\n", Opt->update);
01235 }
01236 brk = YUP;
01237 }
01238
01239 if (!brk && (strcmp(argv[kar], "-prefix") == 0)) {
01240 kar ++;
01241 if (kar >= argc) {
01242 fprintf (SUMA_STDERR, "need argument after -prefix \n");
01243 exit (1);
01244 }
01245 Opt->out_prefix = SUMA_copy_string(argv[kar]);
01246 Opt->WriteFile = YUP;
01247 brk = YUP;
01248 }
01249
01250 if (!brk && (strncmp(argv[kar], "-surf_", 6) == 0)) {
01251 if (kar + 1>= argc) {
01252 fprintf (SUMA_STDERR, "need argument after -surf_X SURF_NAME \n");
01253 exit (1);
01254 }
01255 ind = argv[kar][6] - 'A';
01256 if (ind < 0 || ind >= SURFCLUST_MAX_SURF) {
01257 fprintf (SUMA_STDERR, "-surf_X SURF_NAME option is out of range.\n"
01258 "Only %d surfaces are allowed. \n"
01259 "Must start with surf_A for first surface.\n", SURFCLUST_MAX_SURF);
01260 exit (1);
01261 }
01262 kar ++;
01263 Opt->surf_names[ind] = argv[kar];
01264 Opt->N_surf = ind+1;
01265 brk = YUP;
01266 }
01267
01268 if (!brk && (strcmp(argv[kar], "-input") == 0)) {
01269 kar ++;
01270 if (kar+1 >= argc) {
01271 fprintf (SUMA_STDERR, "need 2 arguments after -input \n");
01272 exit (1);
01273 }
01274 Opt->in_name = argv[kar]; kar ++;
01275
01276 Opt->labelcol = atoi(argv[kar]);
01277 brk = YUP;
01278 }
01279
01280 if (!brk && (strcmp(argv[kar], "-rmm") == 0)) {
01281 kar ++;
01282 if (kar >= argc) {
01283 fprintf (SUMA_STDERR, "need argument after -rmm \n");
01284 exit (1);
01285 }
01286 Opt->DistLim = atof(argv[kar]);
01287 brk = YUP;
01288 }
01289
01290 if (!brk && (strcmp(argv[kar], "-thresh") == 0)) {
01291 kar ++;
01292 if (kar >= argc) {
01293 fprintf (SUMA_STDERR, "need argument after -thresh \n");
01294 exit (1);
01295 }
01296 Opt->DoThreshold = 1;
01297 Opt->Thresh = atof(argv[kar]);
01298 brk = YUP;
01299 }
01300
01301 if (!brk && (strcmp(argv[kar], "-athresh") == 0)) {
01302 kar ++;
01303 if (kar >= argc) {
01304 fprintf (SUMA_STDERR, "need argument after -athresh \n");
01305 exit (1);
01306 }
01307 Opt->DoThreshold = 2;
01308 Opt->Thresh = atof(argv[kar]);
01309 brk = YUP;
01310 }
01311
01312 if (!brk && (strcmp(argv[kar], "-thresh_col") == 0)) {
01313 kar ++;
01314 if (kar >= argc) {
01315 fprintf (SUMA_STDERR, "need argument after -thresh_col \n");
01316 exit (1);
01317 }
01318 Opt->tind = atoi(argv[kar]);
01319 brk = YUP;
01320 }
01321
01322 if (!brk && (strcmp(argv[kar], "-amm2") == 0)) {
01323 kar ++;
01324 if (kar >= argc) {
01325 fprintf (SUMA_STDERR, "need argument after -amm2 \n");
01326 exit (1);
01327 }
01328 Opt->AreaLim = atof(argv[kar]);
01329 brk = YUP;
01330 }
01331
01332 if (!brk && (strcmp(argv[kar], "-out_roidset") == 0)) {
01333 Opt->OutROI = YUP;
01334 brk = YUP;
01335 }
01336
01337 if (!brk && (strcmp(argv[kar], "-out_clusterdset") == 0)) {
01338 Opt->OutClustDset = YUP;
01339 brk = YUP;
01340 }
01341
01342 if (!brk && (strcmp(argv[kar], "-out_fulllist") == 0)) {
01343 Opt->FullROIList = YUP;
01344 brk = YUP;
01345 }
01346
01347 if (!brk && (strcmp(argv[kar], "-sort_none") == 0)) {
01348 Opt->SortMode = SUMA_SORT_CLUST_NO_SORT;
01349 brk = YUP;
01350 }
01351
01352 if (!brk && (strcmp(argv[kar], "-sort_n_nodes") == 0)) {
01353 Opt->SortMode = SUMA_SORT_CLUST_BY_NUMBER_NODES;
01354 brk = YUP;
01355 }
01356
01357 if (!brk && (strcmp(argv[kar], "-sort_area") == 0)) {
01358 Opt->SortMode = SUMA_SORT_CLUST_BY_AREA;
01359 brk = YUP;
01360 }
01361
01362 if (!brk) {
01363 fprintf (SUMA_STDERR,"Error %s:\nOption %s not understood. Try -help for usage\n", FuncName, argv[kar]);
01364 exit (1);
01365 } else {
01366 brk = NOPE;
01367 kar ++;
01368 }
01369 }
01370
01371
01372 if (Opt->DistLim < 0) {
01373 fprintf (SUMA_STDERR, "must use options -rmm \n");
01374 exit(1);
01375 }
01376 if (!Opt->out_prefix) {
01377 Opt->out_prefix = SUMA_RemoveDsetExtension(Opt->in_name, SUMA_NO_DSET_FORMAT);
01378 }
01379
01380 if (Opt->SortMode == SUMA_SORT_CLUST_NOT_SET) { Opt->SortMode = SUMA_SORT_CLUST_BY_AREA; }
01381
01382 if (BuildMethod == SUMA_OFFSETS2) { SUMA_S_Note("Using Offsets2"); }
01383 else if (BuildMethod == SUMA_OFFSETS_LL) { SUMA_S_Note("Using Offsets_ll"); }
01384 else if (BuildMethod == SUMA_OFFSETS2_NO_REC) { if (LocalHead) SUMA_S_Note("Using no recursion"); }
01385 else {
01386 SUMA_SL_Err("Bad BuildMethod");
01387 exit(1);
01388 }
01389
01390 if (Opt->FullROIList && !(Opt->OutROI || Opt->OutClustDset)) {
01391 SUMA_SL_Err("-out_fulllist must be used in conjunction with -out_ROIdset or -out_clusterdset");
01392 exit(1);
01393 }
01394 SUMA_RETURN(Opt);
01395 }
01396
01397 int main (int argc,char *argv[])
01398 {
01399 static char FuncName[]={"SurfClust"};
01400 int kar, SO_read, *ni=NULL, N_ni, cnt, i, *nip=NULL;
01401 float *data_old = NULL, *far = NULL, *nv=NULL, *nt = NULL;
01402 void *SO_name = NULL;
01403 SUMA_SurfaceObject *SO = NULL, *SOnew = NULL;
01404 MRI_IMAGE *im = NULL;
01405 SUMA_DSET_FORMAT iform;
01406 SUMA_SURFCLUST_OPTIONS *Opt;
01407 SUMA_SurfSpecFile Spec;
01408 DList *list = NULL;
01409 SUMA_DSET *dset = NULL;
01410 float *NodeArea = NULL;
01411 FILE *clustout=NULL;
01412 char *ClustOutName = NULL, *params=NULL, stmp[200];
01413 SUMA_Boolean LocalHead = NOPE;
01414
01415 SUMA_mainENTRY;
01416 SUMA_STANDALONE_INIT;
01417
01418
01419
01420 SUMAg_DOv = SUMA_Alloc_DisplayObject_Struct (SUMA_MAX_DISPLAYABLE_OBJECTS);
01421
01422 if (argc < 6)
01423 {
01424 usage_SUMA_SurfClust();
01425 exit (1);
01426 }
01427
01428 Opt = SUMA_SurfClust_ParseInput (argv, argc);
01429
01430 if (Opt->WriteFile) {
01431 if (Opt->AreaLim < 0) sprintf(stmp,"_ClstTable_r%.1f.1D", Opt->DistLim);
01432 else sprintf(stmp,"_ClstTable_r%.1f_a%.1f.1D", Opt->DistLim, Opt->AreaLim);
01433 ClustOutName = SUMA_append_string(Opt->out_prefix, stmp);
01434 if (SUMA_filexists(ClustOutName)) {
01435 fprintf (SUMA_STDERR,"Error %s:\nOutput file %s exists, will not overwrite.\n", FuncName, ClustOutName);
01436 exit(1);
01437 }
01438 }
01439
01440
01441 if (!SUMA_Read_SpecFile (Opt->spec_file, &Spec)) {
01442 fprintf(SUMA_STDERR,"Error %s: Error in SUMA_Read_SpecFile\n", FuncName);
01443 exit(1);
01444 }
01445 SO_read = SUMA_spec_select_surfs(&Spec, Opt->surf_names, SURFCLUST_MAX_SURF, 0);
01446 if ( SO_read != Opt->N_surf )
01447 {
01448 if (SO_read >=0 )
01449 fprintf(SUMA_STDERR,"Error %s:\nFound %d surfaces, expected %d.\n", FuncName, SO_read, Opt->N_surf);
01450 exit(1);
01451 }
01452
01453 if (!SUMA_LoadSpec_eng(&Spec, SUMAg_DOv, &SUMAg_N_DOv, Opt->sv_name, 0, SUMAg_CF->DsetList) ) {
01454 fprintf(SUMA_STDERR,"Error %s: Failed in SUMA_LoadSpec_eng\n", FuncName);
01455 exit(1);
01456 }
01457 SO = SUMA_find_named_SOp_inDOv(Opt->surf_names[0], SUMAg_DOv, SUMAg_N_DOv);
01458 NodeArea = SUMA_CalculateNodeAreas(SO);
01459 if (!NodeArea) {
01460 SUMA_S_Err("Failed to calculate Node Areas.\n");
01461 exit(1);
01462 }
01463
01464 iform = SUMA_NO_DSET_FORMAT;
01465 dset = SUMA_LoadDset (Opt->in_name, &iform, 0);
01466 if (LocalHead) SUMA_ShowDset(dset, 0, NULL);
01467 if (!dset) { SUMA_S_Err( "Failed to load dataset.\n"
01468 "Make sure file exists\n"
01469 "and is of the specified\n"
01470 "format.");
01471 exit(1); }
01472 if (!SUMA_OKassign(dset, SO)) {
01473 SUMA_SL_Err("Failed to assign data set to surface.");
01474 exit(1);
01475 }
01476
01477 nip = SUMA_GetNodeDef(dset);
01478 N_ni = SDEST_VECLEN(dset);
01479 if (!nip) {
01480 SUMA_S_Err("Failed to find node index column");
01481 exit(1);
01482 }
01483
01484 ni = (int *)SUMA_malloc(N_ni*sizeof(int));
01485 memcpy (ni, nip, N_ni*sizeof(int));
01486 nv = SUMA_DsetCol2Float(dset, Opt->labelcol, 0);
01487 if (!nv) {
01488 SUMA_S_Err("Failed to find node value column");
01489 exit(1);
01490 }
01491
01492
01493 if (Opt->DoThreshold) {
01494 nt = SUMA_DsetCol2Float(dset, Opt->tind, 0);
01495 if (!nt) {
01496 SUMA_S_Err("Failed to find threshold column");
01497 exit(1);
01498 }
01499 cnt = 0;
01500 if (Opt->DoThreshold == 1) {
01501 if (Opt->update) fprintf(SUMA_STDERR,"%s: Thresholding at %f...\n", FuncName, Opt->Thresh);
01502 for (i=0;i<N_ni; ++i) {
01503 if (nt[i] >= Opt->Thresh) {
01504 ni[cnt] = ni[i];
01505 nv[cnt] = nv[i];
01506 ++cnt;
01507 }
01508 }
01509 } else {
01510 SUMA_LH("ABS Thresholding...");
01511 for (i=0;i<N_ni; ++i) {
01512 if (fabs(nt[i]) >= Opt->Thresh) {
01513 ni[cnt] = ni[i];
01514 nv[cnt] = nv[i];
01515 ++cnt;
01516 }
01517 }
01518 }
01519 N_ni = cnt;
01520 }
01521 if (Opt->update) {
01522 Opt->update = -(N_ni * Opt->update / 100);
01523 if (LocalHead) {
01524 fprintf(SUMA_STDERR,"Update parameter, once every %d nodes\n%d nodes to work with.\n", -(int)Opt->update, N_ni);
01525 }
01526 }
01527
01528
01529 list = SUMA_FindClusters (SO, ni, nv, N_ni, -1, Opt, NodeArea);
01530
01531
01532 if (!SUMA_Sort_ClustersList (list, Opt->SortMode)) {
01533 SUMA_S_Err("Failed to sort cluster list");
01534 exit(1);
01535 }
01536
01537
01538 params = SUMA_HistString(FuncName, argc, argv, NULL);
01539 if (Opt->WriteFile) {
01540 clustout = fopen(ClustOutName, "w");
01541 if (!clustout) {
01542 fprintf (SUMA_STDERR,"Error %s:\nFailed to open %s for writing.\nCheck permissions.\n", FuncName, ClustOutName);
01543 exit(1);
01544 }
01545 SUMA_Show_SurfClust_list(list, clustout, 0, params);
01546 fclose(clustout);clustout = NULL;
01547 } else SUMA_Show_SurfClust_list(list, NULL, 0, params);
01548
01549
01550 if (Opt->OutROI) {
01551 SUMA_DSET *dset_roi = NULL;
01552 char *ROIprefix = NULL;
01553 char *NameOut = NULL;
01554 if (Opt->AreaLim < 0) sprintf(stmp,"_ClstMsk_r%.1f", Opt->DistLim);
01555 else sprintf(stmp,"_ClstMsk_r%.1f_a%.1f", Opt->DistLim, Opt->AreaLim);
01556 ROIprefix = SUMA_append_string(Opt->out_prefix, stmp);
01557
01558 dset_roi = SUMA_SurfClust_list_2_DsetMask(SO, list, Opt->FullROIList, ROIprefix);
01559
01560 if (!dset_roi) {
01561 SUMA_S_Err("NULL dset_roi");
01562 exit(1);
01563 }
01564 NameOut = SUMA_WriteDset (ROIprefix, dset_roi, SUMA_ASCII_NIML, 0, 0);
01565 if (!NameOut) { SUMA_SL_Err("Failed to write dataset."); exit(1); }
01566 SUMA_FreeDset((void *)dset_roi); dset_roi = NULL;
01567 if (NameOut) SUMA_free(NameOut); NameOut = NULL;
01568 if (ROIprefix) SUMA_free(ROIprefix); ROIprefix = NULL;
01569 }
01570
01571 if (Opt->OutClustDset) {
01572 SUMA_DSET *dset_clust = NULL;
01573 char *Clustprefix = NULL;
01574 char *NameOut = NULL;
01575 if (Opt->AreaLim < 0) sprintf(stmp,"_Clustered_r%.1f", Opt->DistLim);
01576 else sprintf(stmp,"_Clustered_r%.1f_a%.1f", Opt->DistLim, Opt->AreaLim);
01577 Clustprefix = SUMA_append_string(Opt->out_prefix, stmp);
01578
01579
01580 dset_clust = SUMA_MaskDsetByClustList(dset, SO, list, Opt->FullROIList, Clustprefix);
01581 if (!dset_clust) {
01582 SUMA_S_Err("NULL dset_clust");
01583 exit(1);
01584 }
01585 NameOut = SUMA_WriteDset (Clustprefix, dset_clust, SUMA_ASCII_NIML, 0, 0);
01586 if (!NameOut) { SUMA_SL_Err("Failed to write dataset."); exit(1); }
01587 SUMA_FreeDset((void *)dset_clust); dset_clust = NULL;
01588 if (NameOut) SUMA_free(NameOut); NameOut = NULL;
01589 if (Clustprefix) SUMA_free(Clustprefix); Clustprefix = NULL;
01590 }
01591
01592 if (ClustOutName) SUMA_free(ClustOutName); ClustOutName = NULL;
01593 if (list) dlist_destroy(list); SUMA_free(list); list = NULL;
01594 if (ni) SUMA_free(ni); ni = NULL;
01595 if (nv) SUMA_free(nv); nv = NULL;
01596 if (nt) SUMA_free(nt); nt = NULL;
01597 if (Opt->out_prefix) SUMA_free(Opt->out_prefix); Opt->out_prefix = NULL;
01598 if (Opt) SUMA_free(Opt);
01599 if (dset) SUMA_FreeDset((void *)dset); dset = NULL;
01600 if (!SUMA_Free_Displayable_Object_Vect (SUMAg_DOv, SUMAg_N_DOv)) {
01601 SUMA_SL_Err("DO Cleanup Failed!");
01602 }
01603 exit(0);
01604 }
01605 #endif