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 /* these global variables must be declared even if they will not be used by this main */
00011 SUMA_SurfaceViewer *SUMAg_cSV = NULL; /*!< Global pointer to current Surface Viewer structure*/
00012 SUMA_SurfaceViewer *SUMAg_SVv = NULL; /*!< Global pointer to the vector containing the various Surface Viewer Structures 
00013                                     SUMAg_SVv contains SUMA_MAX_SURF_VIEWERS structures */
00014 int SUMAg_N_SVv = 0; /*!< Number of SVs realized by X */
00015 SUMA_DO *SUMAg_DOv = NULL;   /*!< Global pointer to Displayable Object structure vector*/
00016 int SUMAg_N_DOv = 0; /*!< Number of DOs stored in DOv */
00017 SUMA_CommonFields *SUMAg_CF = NULL; /*!< Global pointer to structure containing info common to all viewers */
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    \brief Calculate area of each node as one third of the sum of
00046    the areas of incident triangles. 
00047    If you change this function make sure changes are also done
00048    on RickR's compute_node_areas since the two functions use
00049    the same principle.
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    \brief builds a cluster starting from some node 
00085    
00086    \param dothisnode (int) start building from this node
00087    \param AddToThisClust (SUMA_CLUST_DATUM *)pointer to cluster to add to. Send NULL
00088                     when function is first called. This is a non NULL when the 
00089                     function recurses.
00090    \param ToBeAssigned (float *) if ToBeAssigned[i] then node i (index into SO's nodelist)
00091                                  is considered in the clustering and the value at that
00092                                  node is ToBeAssigned[i]. Vector is SO->N_Node elements
00093                                  long. Gets modified during recursion.
00094    \param N_ToBeAssigned (int *) pointer to number of value values in ToBeAssigned. 
00095                                  Changes with recursive calls. 
00096    \param NodeArea(float *) Vector containing area of each node of surface
00097    \param SO (SUMA_SurfaceObject *) the usual deal
00098    \param Opt (SUMA_SURFCLUST_OPTIONS *) options for cluster building.     
00099    
00100    \sa recursive calls SUCK a lot of memory and are a nightmare to track with 
00101    function I/O. They are also slow and not necessary! 
00102    
00103    \sa SUMA_Build_Cluster_From_Node_NoRec       
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; /* Clust->rank = -1; */  
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    /* Add node to cluster */
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    /* mark it as assigned, an reduce the number of nodes left to assign*/
00160    ToBeAssigned[dothisnode] = 0; --(*N_TobeAssigned);
00161    
00162    if (BuildMethod == SUMA_OFFSETS2) {
00163       /* Tres bad memory utilization due to recursive calls */
00164       if (*N_TobeAssigned) {
00165          /* look in its vicinity - bad memory usage due to recursive calls*/
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          /* search to see if any are to be assigned */
00183          for (il=1; il<OffS->N_layers; ++il) { /* starting at layer 1, layer 0 is the node itself */
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                      /* take that node into the cluster */
00188                      SUMA_Build_Cluster_From_Node( neighb, Clust, ToBeAssigned, N_TobeAssigned, NodeArea, SO, Opt); 
00189                }
00190             }
00191          }
00192          /* free this OffS structure (Note you can't recycle the same structure because you are using many
00193          OffS at one because of recursive calls */
00194          if (OffS) SUMA_Free_getoffsets(OffS); OffS = NULL;
00195       }
00196    } else if (BuildMethod == SUMA_OFFSETS_LL) { 
00197       if (*N_TobeAssigned) {
00198          /* look in its vicinity */
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          /* search to see if any are to be assigned, start at layer 1*/
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                /* take that node into the cluster */
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    /* If you ever use this function again, you should probably SUMA_free(offlist); dlist_destroy(candlist); SUMA_free(candlist); */
00239 
00240    SUMA_RETURN(Clust);
00241 }
00242 
00243 /*! 
00244    A macro for SUMA_Build_Cluster_From_Node_NoRec
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    \brief builds a cluster starting from some node 
00259    
00260    \param dothisnode (int) start building from this node
00261    \param ToBeAssigned (float *) if ToBeAssigned[i] then node i (index into SO's nodelist)
00262                                  is considered in the clustering and the value at that
00263                                  node is ToBeAssigned[i]. Vector is SO->N_Node elements
00264                                  long. Gets modified in function.
00265    \param N_ToBeAssigned (int *) pointer to number of value values in ToBeAssigned. 
00266                                  Gets modified in function. 
00267    \param NodeArea(float *) Vector containing area of each node of surface
00268    \param SO (SUMA_SurfaceObject *) the usual deal
00269    \param Opt (SUMA_SURFCLUST_OPTIONS *) options for cluster building.     
00270       
00271    \sa Based on recursive SUMA_Build_Cluster_From_Node       
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    /* a trick to know when to initialize */
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; /* Clust->rank = -1; */  
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    /* Add node to cluster */
00323    SUMA_ADD_NODE_TO_CLUST(dothisnode, Clust, NodeArea, ToBeAssigned);
00324    /* mark it as assigned, an reduce the number of nodes left to assign*/
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          /* look in its vicinity */
00330          dothiselm = dlist_head(candlist); dothisnode = (int) dothiselm->data;
00331          SUMA_getoffsets2 (dothisnode, SO, Opt->DistLim, OffS, NULL, 0);
00332          /* remove node from candidate list */
00333          dlist_remove(candlist, dothiselm, (void*)&itmp);
00334          /* search to see if any are to be assigned */
00335          for (il=1; il<OffS->N_layers; ++il) { /* starting at layer 1, layer 0 is the node itself */
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                      /* take that node into the cluster */
00340                      SUMA_ADD_NODE_TO_CLUST(neighb, Clust, NodeArea, ToBeAssigned);
00341                      /* mark it as assigned, an reduce the number of nodes left to assign*/
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                      /* mark it as a candidate if it has not been visited as a candidate before */
00352                      if (!visited[neighb]) {
00353                         dlist_ins_next(candlist, dlist_tail(candlist), (void *)neighb);
00354                         visited[neighb] = YUP;   
00355                      }
00356                      
00357                }
00358             }
00359          }
00360          /* recycle */
00361          SUMA_Recycle_getoffsets (OffS);
00362       }
00363    /* free this OffS structure  */
00364    if (OffS) SUMA_Free_getoffsets(OffS); OffS = NULL;
00365    if (Opt->update) fprintf(SUMA_STDERR,"\n");
00366    /* destroy the list */
00367    if (candlist) { dlist_destroy(candlist); SUMA_free(candlist); candlist  = NULL; }
00368    SUMA_RETURN(Clust);
00369 }
00370 
00371 /*!
00372    \brief Finds clusters of data on the surface.
00373    
00374    \param SO (SUMA_SurfaceObject *) pointer to surface in question
00375    \param ni (int *) vector of N_Ni node indices that are to be considered
00376                      for clustering
00377    \param nv (float *) vector of N_Ni node values (corresponding to ni) 
00378                       Only non 0 values are considered for clustering
00379    \param N_ni (int) number of values in ni and nv
00380    \param dothisnode (int) index of node (into SO's nodelist, not ni) to start/proceed from
00381    \param Opt (SUMA_SURFCLUST_OPTIONS *) structure containing clustering options
00382    \param AddToThisClust (SUMA_CLUST_DATUM *) add to this cluster
00383    \param NodeArea (float *) SO->N_Node vector of node areas. 
00384    \return ClustList (DList *) list of clusters, in no particular order. Processing of the list is to
00385                               be done later.
00386    
00387    \sa SUMA_Build_Cluster_From_Node
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) { /* initialize */
00408       SUMA_LH("Initializing");
00409       nc = 0;
00410       /* initialize the list */
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          /* reallocate to save space */
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          /* find the central node */
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 /*! Show the ViewState structure */
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 /*! Show the SurfClust list */
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 /*! Masks a data set by a clustering list*/
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    /* which nodes in mask ? */
00631    ismask = (byte *)SUMA_calloc(SO->N_Node, sizeof(byte)); /* you need to allocate for SO->N_Node, to be safe, otherwise you will need
00632                                                             to search for the highest node index.. */
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    /* now form a rowmask vector to parallel rows in idset */
00647    rowmask = (byte *)SUMA_calloc(idset->dnel->vec_len, sizeof(byte));
00648    colmask = (byte *)SUMA_calloc(idset->dnel->vec_num, sizeof(byte));
00649    
00650    /* get the node index column in idset */
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    /* now, fill rowmask */
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    /* fill colmask*/
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 /*! Turn the clusters to a cluster dataset mask*/
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,         /* usually the filename */
00744                                  SUMA_NODE_ROI,                /* mix and match */
00745                                  NULL,    /* no idcode, let the function create one from the filename*/
00746                                  NULL,       /* no domain str specified */
00747                                  N_Node    /* Number of nodes allocated for */
00748                                  ); /* DO NOT free dset, it is store in DsetList */
00749                            
00750         /* form the dataset */
00751    SUMA_LH("Adding NodeDef column ...");
00752    if (!SUMA_AddDsetNelCol (   dset, /* the famed nel */ 
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    \brief Finds a node that best approximates the property of the center of mass on the surface.
00776    Note that the real center of mass of a curved surface is rarely on the surface, so this is 
00777    an attempt at localizing a node that is central to an ROI and that can reflect the weight, 
00778    or activity, distribution in that ROI.
00779    
00780    \param SO (SUMA_SurfaceObject *)
00781    \param cd (SUMA_CLUST_DATUM *) basically the output of SUMA_SUMA_Build_Cluster_From_Node or 
00782                                   an element of the list returned by SUMA_FindClusters
00783                                   Function will fill centralnode and weightedcentralnode
00784    \param UseSurfDist (int) 0: use distances along the surface (approximated by distances on the graph)
00785                             1: use Euclidian distances. 
00786    \return ans (int) : NOPE failed, YUP succeeded.
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       /* not much to do, return 0 */
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]; /* anchor node */
00829    SUMA_getoffsets2 (cd->NodeList[0], SO, 0, OffS, CoverThisNode, cd->N_Node);
00830    /* help me with a nicer structure (for sanity)*/
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    /* put the distance in an easy to search array */
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    /* Now calculate the center of massity of each node 
00847    This is no center of mass in the proper definition of the term*/
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]; /* Node index into SO of center node */
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; /* recalculated each time, not a big deal ... */
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]; /* Node index into SO of other node in cluster */
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]; /* Node index into SO of center node */
00891          s[0] = s[1] = s[2] = 0.0;
00892          centrality[c] = 0.0; mtotal = 0.0; /* recalculated each time, not a big deal ... */
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]; /* Node index into SO of other node in cluster */
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             /* compare to all ahead of this element */
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) { /* max is the real deal */
00968                   dlist_remove(list, max_elmt, (void *)(&cd_max));
00969                   dlist_ins_prev(list, elmt, (void *)cd_max);
00970                   elmt = elmt->prev; /* continue from that one */
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"/* [<-sv sv_name>] */
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                /*"     -sv SurfaceVolume \n"
01069                "        If you supply a surface volume, the coordinates of the input surface.\n"
01070                "        are modified to SUMA's convention and aligned with SurfaceVolume.\n"
01071                "     -acpc: Apply acpc transform (which must be in acpc version of \n"
01072                "        SurfaceVolume) to the surface vertex coordinates. \n"
01073                "        This option must be used with the -sv option.\n"
01074                "     -tlrc: Apply Talairach transform (which must be a talairach version of \n"
01075                "        SurfaceVolume) to the surface vertex coordinates. \n"
01076                "        This option must be used with the -sv option.\n"
01077                "     -MNI_rai/-MNI_lpi: Apply Andreas Meyer Lindenberg's transform to turn \n"
01078                "        AFNI tlrc coordinates (RAI) into MNI coord space \n"
01079                "        in RAI (with -MNI_rai) or LPI (with -MNI_lpi)).\n"
01080                "        NOTE: -MNI_lpi option has not been tested yet (I have no data\n"
01081                "        to test it on. Verify alignment with AFNI and please report\n"
01082                "        any bugs.\n" 
01083                "        This option can be used without the -tlrc option.\n"
01084                "        But that assumes that surface nodes are already in\n"
01085                "        AFNI RAI tlrc coordinates .\n"    
01086                "\n" */
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                /* Can use -O2_NR, -O2 and -Oll to try different implementations 
01121                of the cluster building function. -O2 and -Oll are recursive and
01122                work well for small clusters, a catastrophy for large clusters.
01123                -O2_NR is fast and reliable. */
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    \brief parse the arguments for SurfSmooth program
01132    
01133    \param argv (char *)
01134    \param argc (int)
01135    \return Opt (SUMA_SURFCLUST_OPTIONS *) options structure.
01136                To free it, use 
01137                SUMA_free(Opt->out_name); 
01138                SUMA_free(Opt);
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) { /* loop accross command ine options */
01178                 /*fprintf(stdout, "%s verbose: Parsing command line...\n", FuncName);*/
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          /* no need for that one Opt->nodecol = atoi(argv[kar]); kar ++; */
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    /* sanitorium */
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 {/* Main */    
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    /* Allocate space for DO structure */
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    /* read  surfaces */
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    /* now read into SUMAg_DOv */
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    /* load the data */   
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    /* get the node index column */
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    /* copy nip's contents because you will be modifying in the thresholding below */
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    /* any thresholding ? */
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); /* make it negative before you begin a clustering operation */
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    /* make the call */
01529    list = SUMA_FindClusters (SO, ni, nv, N_ni, -1, Opt, NodeArea);
01530    
01531    /* sort the list */
01532    if (!SUMA_Sort_ClustersList (list, Opt->SortMode)) {
01533       SUMA_S_Err("Failed to sort cluster list");
01534       exit(1);
01535    }
01536        
01537    /* Show the results */
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       /* Call this function, write out the resultant dset to disk then cleanup */
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       /* Call this function, write out the resultant dset to disk then cleanup */
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