00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2001, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   This file contains routines for hierarchical clustering of user specified 
00010   parametric data.
00011 
00012 
00013   File:    StatClust.c
00014   Author:  B. Douglas Ward
00015   Date:    08 October 1999
00016 
00017   Mod:     Replaced C "pow" function, significantly improving execution speed.
00018   Date:    11 October 1999
00019 
00020   Mod:     Replaced dataset interface code with call to THD_open_dataset.
00021            Restructured code for initializing hierarchical clustering.
00022   Date:    19 October 1999
00023 
00024   Mod:     At each output cluster agglomeration step, print to the screen
00025            which clusters are to be combined, and their distance.
00026   Date:    05 September 2000
00027 
00028   Mod:     Corrected error in sort_clusters routine.
00029   Date:    30 April 2001
00030 
00031 */
00032 
00033 /*---------------------------------------------------------------------------*/
00034 /*
00035   Structure declarations 
00036 */
00037 
00038 
00039 struct cluster;
00040 struct voxel;
00041 
00042 
00043 typedef struct cluster
00044 {
00045   float * centroid;
00046   int num_voxels;
00047   struct voxel * voxel_ptr;
00048   float nearest_dist;
00049   struct cluster * nearest_cluster;
00050   struct cluster * next_cluster;
00051 } cluster;
00052 
00053 
00054 typedef struct voxel
00055 {
00056   int index;
00057   struct voxel * next_voxel;
00058 } voxel;
00059 
00060 
00061 /*---------------------------------------------------------------------------*/
00062 /*
00063   Create a new voxel having the given index number.
00064 */
00065 
00066 voxel * new_voxel (int index)
00067 {
00068   voxel * voxel_ptr = NULL;
00069   
00070   voxel_ptr = (voxel *) malloc (sizeof(voxel));
00071   MTEST (voxel_ptr);
00072 
00073   voxel_ptr->index = index;
00074   voxel_ptr->next_voxel = NULL;
00075   return (voxel_ptr);
00076 }
00077 
00078 
00079 /*---------------------------------------------------------------------------*/
00080 /*
00081   Print the index number for this voxel.
00082 */
00083 
00084 void print_voxel (voxel * voxel_ptr)
00085 {
00086   if (voxel_ptr != NULL)
00087     printf ("%d ", voxel_ptr->index);
00088 }
00089 
00090 
00091 /*---------------------------------------------------------------------------*/
00092 /*
00093   Print all voxels in the linked list.
00094 */
00095 
00096 void print_all_voxels (voxel * voxel_ptr)
00097 {
00098   while (voxel_ptr != NULL)
00099     {
00100       print_voxel (voxel_ptr);
00101       voxel_ptr = voxel_ptr->next_voxel;
00102     }
00103   printf ("\n");
00104 }
00105 
00106 
00107 /*---------------------------------------------------------------------------*/
00108 /*
00109   Create an empty cluster.
00110 */
00111   
00112 cluster * initialize_cluster ()
00113 {
00114   cluster * clust_ptr = NULL;
00115 
00116   clust_ptr = (cluster *) malloc (sizeof(cluster));
00117   MTEST (clust_ptr);
00118 
00119   clust_ptr->next_cluster = NULL;
00120   clust_ptr->num_voxels = 0;
00121   clust_ptr->voxel_ptr = NULL;
00122   clust_ptr->centroid = NULL;
00123   clust_ptr->nearest_dist = 0.0;
00124   clust_ptr->nearest_cluster = NULL;
00125   return (clust_ptr);
00126   
00127 }
00128 
00129 
00130 /*---------------------------------------------------------------------------*/
00131 /*
00132   Print the contents of one cluster.
00133 */
00134 
00135 void print_cluster (cluster * clust_ptr, char * str, matrix s)
00136 {
00137   int i;
00138   vector v, sv;
00139   vector_initialize (&v);
00140   vector_initialize (&sv);
00141 
00142   printf ("Cluster %s \n", str);
00143 
00144   if (clust_ptr->voxel_ptr != NULL)
00145     {
00146       printf ("# Voxels = %d \n", clust_ptr->num_voxels);
00147 
00148       /*
00149       printf ("Voxels: ");
00150       
00151       print_all_voxels (clust_ptr->voxel_ptr);
00152       */
00153     }
00154 
00155   if (clust_ptr->centroid != NULL)
00156     {
00157       printf ("Centroid: \n");
00158       array_to_vector (SC_dimension, clust_ptr->centroid, &v);
00159       vector_multiply (s, v, &sv);
00160       vector_print (sv);
00161     }
00162 
00163   /*
00164   printf ("Nearest cluster distance = %f \n", clust_ptr->nearest_dist);
00165   */
00166 
00167   vector_destroy (&v);
00168   vector_destroy (&sv);
00169 }
00170 
00171 
00172 /*---------------------------------------------------------------------------*/
00173 /*
00174   Print the contents of all clusters in the linked list.
00175 */
00176 
00177 void print_all_clusters (cluster * clust_ptr, matrix s)
00178 {
00179   int iclust = 0;
00180   char str[30];
00181 
00182   while (clust_ptr != NULL)
00183     {
00184       iclust++;
00185       sprintf (str, "#%d", iclust);
00186       print_cluster (clust_ptr, str, s);
00187       clust_ptr = clust_ptr->next_cluster;
00188     }
00189 
00190 }
00191 
00192 
00193 /*---------------------------------------------------------------------------*/
00194 /*
00195   Save voxel contents of this cluster into byte array (sub-brick).
00196 */
00197 
00198 void save_cluster (cluster * clust_ptr, byte iclust, byte * bar)
00199 {
00200   int i;
00201   voxel * voxel_ptr = NULL;
00202 
00203 
00204   voxel_ptr = clust_ptr->voxel_ptr;
00205 
00206 
00207   while (voxel_ptr != NULL)
00208     {
00209       bar[voxel_ptr->index] = iclust;
00210       voxel_ptr = voxel_ptr->next_voxel;
00211     }
00212 
00213 }
00214 
00215 
00216 /*---------------------------------------------------------------------------*/
00217 /*
00218   Save voxel contents of all clusters into byte array (sub-brick).
00219 */
00220 
00221 void save_all_clusters (cluster * clust_ptr, byte * bar)
00222 {
00223   byte iclust = 0;
00224 
00225   while (clust_ptr != NULL)
00226     {
00227       iclust++;
00228       save_cluster (clust_ptr, iclust, bar);
00229       clust_ptr = clust_ptr->next_cluster;
00230     }
00231 
00232 }
00233 
00234 
00235 /*---------------------------------------------------------------------------*/
00236 /*
00237   Calculate the distance between two clusters as the Euclidean distance
00238   between their centroids.
00239 */
00240 
00241 float cluster_distance (cluster * aclust, cluster * bclust)
00242 {
00243   float sumsqr;
00244   float delta;
00245   int i;
00246 
00247   sumsqr = 0.0;
00248   for (i = 0;  i < SC_dimension;  i++)
00249     {
00250       delta = aclust->centroid[i] - bclust->centroid[i];
00251       sumsqr += delta * delta;
00252     }
00253 
00254   return (sqrt(sumsqr));
00255   
00256 }
00257 
00258 
00259 /*---------------------------------------------------------------------------*/
00260 /*
00261   Find the cluster which is nearest to new_cluster.
00262   Set the nearest_dist and nearest_cluster structure elements accordingly.
00263 */
00264   
00265 void find_nearest_cluster (cluster * new_cluster, cluster * head_cluster)
00266 {
00267   const float MAX_DIST = 1.0e+30;
00268 
00269   cluster * clust_ptr = NULL;
00270   float dist;
00271 
00272 
00273   /*----- Initialize nearest cluster elements -----*/
00274   new_cluster->nearest_dist = MAX_DIST;
00275   new_cluster->nearest_cluster = NULL;
00276 
00277 
00278   clust_ptr = head_cluster;
00279 
00280   while (clust_ptr != NULL)
00281     {
00282       if (clust_ptr != new_cluster)
00283         {
00284           dist = cluster_distance (new_cluster, clust_ptr);
00285 
00286           if (dist < new_cluster->nearest_dist)
00287             {
00288               new_cluster->nearest_dist = dist;
00289               new_cluster->nearest_cluster = clust_ptr;
00290             }
00291 
00292           if (dist < clust_ptr->nearest_dist)
00293             {
00294               clust_ptr->nearest_dist = dist;
00295               clust_ptr->nearest_cluster = new_cluster;
00296             }
00297         }
00298 
00299       clust_ptr = clust_ptr->next_cluster;
00300     }
00301 }
00302 
00303 
00304 /*---------------------------------------------------------------------------*/
00305 /*
00306   Add a new cluster to the linked list of clusters.
00307 */
00308 
00309 void add_cluster (cluster * new_cluster, cluster * head_cluster)
00310 {
00311 
00312   new_cluster->next_cluster = head_cluster;
00313 
00314   find_nearest_cluster (new_cluster, head_cluster);
00315 
00316 }
00317 
00318 
00319 /*---------------------------------------------------------------------------*/
00320 /*
00321   Create a new cluster containing a single voxel, and add to list of clusters.
00322 */
00323 
00324 cluster * new_cluster (int index, float * centroid, cluster * head_clust)
00325 {
00326   cluster * clust_ptr = NULL;
00327   voxel * voxel_ptr = NULL;
00328 
00329   clust_ptr = initialize_cluster ();
00330 
00331   clust_ptr->num_voxels = 1;
00332   clust_ptr->voxel_ptr = new_voxel(index);
00333   clust_ptr->centroid = centroid;
00334 
00335   add_cluster (clust_ptr, head_clust);
00336 
00337   return (clust_ptr);
00338   
00339 }
00340 
00341 
00342 /*---------------------------------------------------------------------------*/
00343 /*
00344   Deallocate memory for this cluster (excluding list of voxels).
00345 */
00346 
00347 void delete_cluster (cluster * clust_ptr)
00348 {
00349   if (clust_ptr != NULL)
00350     {
00351       clust_ptr->voxel_ptr = NULL;
00352 
00353       if (clust_ptr->centroid != NULL)
00354         {
00355           free (clust_ptr->centroid);
00356           clust_ptr->centroid = NULL;
00357         }
00358 
00359       free (clust_ptr);
00360     }
00361 }
00362 
00363 
00364 /*---------------------------------------------------------------------------*/
00365 /*
00366   Deallocate memory for this cluster.
00367 
00368   Note:  This routine does not actually delete the list of voxels.
00369 */
00370 
00371 void destroy_cluster (cluster * clust_ptr)
00372 {
00373   if (clust_ptr != NULL)
00374     {
00375       if (clust_ptr->voxel_ptr != NULL)
00376         free (clust_ptr->voxel_ptr);
00377   
00378       delete_cluster (clust_ptr);
00379     }
00380 }
00381 
00382 
00383 /*---------------------------------------------------------------------------*/
00384 /*
00385   Remove two clusters from linked list of clusters.
00386   Reset cluster pointers, and recalculate nearest cluster distances 
00387   where needed.  Finally, delete the two clusters from memory.
00388 */
00389 
00390 cluster * remove_clusters (cluster * aclust, cluster * bclust, 
00391                            cluster * head_clust)
00392 {
00393   cluster * clust_ptr = NULL;
00394   cluster * next_clust = NULL;
00395   
00396 
00397   while ((head_clust != NULL) && 
00398          ((head_clust == aclust) || (head_clust == bclust)))
00399     head_clust = head_clust->next_cluster;
00400 
00401 
00402   if (head_clust != NULL)
00403     {
00404 
00405       clust_ptr = head_clust;
00406       next_clust = clust_ptr->next_cluster;
00407       while (next_clust != NULL)
00408         {
00409           if ((next_clust == aclust) || (next_clust == bclust))
00410             clust_ptr->next_cluster = next_clust->next_cluster;
00411           else
00412             clust_ptr = next_clust;
00413           
00414           next_clust = clust_ptr->next_cluster;
00415         }
00416 
00417 
00418       clust_ptr = head_clust;
00419       while (clust_ptr != NULL)
00420         {
00421           if ((clust_ptr->nearest_cluster == aclust) 
00422               || (clust_ptr->nearest_cluster == bclust))
00423             {
00424               find_nearest_cluster (clust_ptr, head_clust);
00425             }
00426           clust_ptr = clust_ptr->next_cluster;
00427         }
00428     }
00429 
00430 
00431   delete_cluster (aclust);
00432   delete_cluster (bclust);
00433 
00434   return (head_clust);
00435 }
00436 
00437 
00438 /*---------------------------------------------------------------------------*/
00439 /*
00440   Merge two clusters into one new cluster.
00441 */
00442 
00443 cluster * merge_clusters (cluster * aclust, cluster * bclust)
00444 {
00445   cluster * abclust = NULL;
00446   voxel * voxel_ptr = NULL;
00447   int na, nb;
00448   int i;
00449 
00450   abclust = initialize_cluster ();
00451 
00452   na = aclust->num_voxels;
00453   nb = bclust->num_voxels;
00454   abclust->num_voxels = na + nb;
00455 
00456   abclust->centroid = (float *) malloc (sizeof(float) * SC_dimension);
00457   MTEST (abclust->centroid);
00458 
00459   for (i = 0;  i < SC_dimension;  i++)
00460     abclust->centroid[i] 
00461       = (na*aclust->centroid[i] + nb*bclust->centroid[i]) / (na+nb);
00462 
00463   abclust->voxel_ptr = aclust->voxel_ptr;
00464 
00465   voxel_ptr = abclust->voxel_ptr;
00466   while (voxel_ptr->next_voxel != NULL)
00467     voxel_ptr = voxel_ptr->next_voxel;
00468   voxel_ptr->next_voxel = bclust->voxel_ptr;
00469   
00470   return (abclust);
00471 }
00472 
00473 
00474 /*---------------------------------------------------------------------------*/
00475 /*
00476   Consolidate two clusters.
00477 */
00478   
00479 cluster * consolidate_clusters (cluster * aclust, cluster * bclust, 
00480                                 cluster * head_clust)
00481 {
00482   cluster * abclust = NULL;
00483 
00484 
00485   /*----- Merge two clusters into one new cluster -----*/
00486   abclust = merge_clusters (aclust, bclust);
00487 
00488 
00489   /*----- Remove the two original clusters from the linked list -----*/
00490   head_clust = remove_clusters (aclust, bclust, head_clust);
00491 
00492 
00493   /*----- Add the merged cluster to the linked list -----*/
00494   add_cluster (abclust, head_clust);
00495   
00496 
00497   /*----- Merged cluster is now at the top of the list -----*/
00498   return (abclust);
00499 }
00500 
00501 
00502 /*---------------------------------------------------------------------------*/
00503 /*
00504   Agglomerate clusters by merging the two clusters which are closest together.
00505 */
00506 
00507 cluster * agglomerate_clusters (cluster * head_clust, int print_flag)
00508 {
00509   const float MAX_DIST = 1.0e+30;
00510 
00511   cluster * clust_ptr = NULL;
00512   cluster * aclust    = NULL;
00513   cluster * bclust    = NULL;
00514   float min_dist;
00515 
00516 
00517   /*----- Find the two clusters which are closest together -----*/
00518   min_dist = MAX_DIST;
00519   clust_ptr = head_clust;
00520   while (clust_ptr != NULL)
00521     {
00522       if (clust_ptr->nearest_dist < min_dist)
00523         {
00524           min_dist = clust_ptr->nearest_dist;
00525           aclust = clust_ptr;
00526           bclust = clust_ptr->nearest_cluster;
00527         } 
00528       clust_ptr = clust_ptr->next_cluster;
00529     }
00530 
00531 
00532   /*----- Identify clusters which are to be merged -----*/
00533   if (print_flag)
00534     {
00535       int iclust, iaclust, ibclust;
00536 
00537       clust_ptr = head_clust;
00538       iclust = 0;
00539       while (clust_ptr != NULL)
00540         {
00541           iclust++;
00542           if (aclust == clust_ptr)  iaclust = iclust;
00543           if (bclust == clust_ptr)  ibclust = iclust;
00544           clust_ptr = clust_ptr->next_cluster;
00545         }
00546       
00547       printf ("Merging cluster #%d and cluster #%d \n", iaclust, ibclust);
00548       printf ("Distance = %f \n", min_dist);
00549     }
00550 
00551 
00552   /*----- Merge these two clusters -----*/
00553   head_clust = consolidate_clusters (aclust, bclust, head_clust);
00554 
00555 
00556   return (head_clust);
00557 }
00558 
00559 
00560 /*---------------------------------------------------------------------------*/
00561 /*
00562   Sort clusters in order of size.
00563 */
00564 
00565 cluster * sort_clusters (cluster * head_clust)
00566 {
00567   cluster * i  = NULL; 
00568   cluster * ip = NULL; 
00569   cluster * m  = NULL;
00570   cluster * mp = NULL;
00571   cluster * j  = NULL;
00572   cluster * jp = NULL;
00573   cluster * guard = NULL;
00574 
00575 
00576   /*----- Create guard cluster in case head cluster must be replaced -----*/
00577   guard = initialize_cluster();
00578   guard->next_cluster = head_clust;
00579   ip = guard;
00580 
00581   while (ip->next_cluster != NULL)
00582     {
00583       /*----- Initialize search for next largest cluster -----*/
00584       i = ip->next_cluster;  /* current top of list */
00585       mp = ip;               /* cluster pointing to next largest cluster */
00586       m = i;                 /* next largest cluster */
00587       jp = i;
00588 
00589       /*----- Search through list for next largest cluster -----*/
00590       while (jp->next_cluster != NULL)
00591         {
00592           j = jp->next_cluster;
00593           if (j->num_voxels > m->num_voxels)
00594             {
00595               mp = jp;
00596               m = j;
00597             }
00598           jp = j;
00599         }
00600 
00601       /*----- Now move next largest cluster to top of list -----*/
00602       if (m != i)
00603         {
00604           ip->next_cluster = m;
00605           mp->next_cluster = m->next_cluster;
00606           m->next_cluster = i;
00607           i = m;
00608         }
00609 
00610       /*----- Move down the list -----*/
00611       ip = i;
00612         
00613     }
00614 
00615   
00616   /*----- Replace head cluster -----*/
00617   head_clust = guard->next_cluster;
00618   delete_cluster (guard);
00619 
00620   return (head_clust);
00621 }
00622 
00623 
00624 /*---------------------------------------------------------------------------*/
00625 /*
00626   Calculate covariance matrix for input parameters.
00627   Return square root of covariance matrix, and its inverse.
00628 */
00629 
00630 void calc_covariance 
00631 (
00632   matrix * s,                /* square root of covariance matrix */
00633   matrix * sinv              /* inverse of square root of covariance matrix */
00634 )
00635 
00636 {
00637   int ivox;                  /* voxel indices */
00638   int ip, jp;                /* parameter indices */
00639   int ok;                    /* Boolean for successful matrix calc. */
00640 
00641   vector mean;               /* mean parameter vector */ 
00642   matrix covar;              /* variance-covariance matrix */
00643   matrix cinv;               /* inverse of covariance matrix */
00644 
00645   char message[80];          /* error message */
00646 
00647 
00648   /*----- Initialize vectors and matrices -----*/
00649   vector_initialize (&mean);
00650   matrix_initialize (&covar);
00651   matrix_initialize (&cinv);
00652 
00653 
00654   /*----- Allocate space for mean and covariance matrices -----*/
00655   vector_create (SC_dimension, &mean);
00656   matrix_create (SC_dimension, SC_dimension, &covar);
00657 
00658 
00659   /*----- Calculate parameter sums and products  -----*/
00660   for (ivox = 0;  ivox < SC_nvox;  ivox++)
00661     for (ip = 0;  ip < SC_dimension;  ip++)
00662       {
00663         mean.elts[ip] += SC_parameters[ip][ivox];
00664         for (jp = 0;  jp < SC_dimension;  jp++)
00665           if ((ip == jp) || (SC_statdist == 2))
00666             covar.elts[ip][jp] += 
00667               SC_parameters[ip][ivox] * SC_parameters[jp][ivox];
00668       }
00669 
00670 
00671   /*----- Calculate the mean parameter vector -----*/
00672   for (ip = 0;  ip < SC_dimension;  ip++)
00673     mean.elts[ip] = mean.elts[ip] / SC_nvox;
00674   if (SC_verb)  
00675     vector_sprint ("Mean parameter vector: ", mean);
00676       
00677 
00678   /*----- Calculate the covariance matrix -----*/
00679   for (ip = 0;  ip < SC_dimension;  ip++)
00680     for (jp = 0;  jp < SC_dimension;  jp++)
00681       if ((ip == jp) || (SC_statdist == 2)) 
00682         covar.elts[ip][jp] = (covar.elts[ip][jp] 
00683                               - SC_nvox * mean.elts[ip] * mean.elts[jp])
00684           / (SC_nvox - 1);
00685   if (SC_verb)
00686     if (SC_statdist == 1)
00687       matrix_sprint ("Parameter variance (diagonal) matrix: ", covar);
00688     else
00689       matrix_sprint ("Parameter covariance matrix: ", covar);
00690 
00691   /*----- Note:  The following sequence of calculations is necessary 
00692     in order to generate an error message in the event of 
00693     perfectly correlated input parameters -----*/
00694 
00695   /*----- Calculate inverse of covariance matrix -----*/
00696   ok = matrix_inverse (covar, &cinv);
00697   if (! ok)  
00698     SC_error 
00699       ("Unable to calculate inverse of covariance matrix");
00700   
00701   /*----- Calculate square root of inverse covariance matrix -----*/
00702   ok = matrix_sqrt (cinv, sinv);
00703   if (! ok)  
00704     SC_error 
00705       ("Unable to calculate square root of inverse of covariance matrix");
00706   
00707   /*----- Calculate square root of covariance matrix -----*/
00708   ok = matrix_inverse (*sinv, s);
00709   if (! ok)  
00710     SC_error 
00711       ("Unable to calculate square root of covariance matrix");
00712   
00713 
00714   /*----- Deallocate memory -----*/
00715   vector_destroy (&mean);
00716   matrix_destroy (&covar);
00717   matrix_destroy (&cinv);
00718 
00719   
00720 }
00721 
00722 
00723 /*---------------------------------------------------------------------------*/
00724 
00725 
00726 
00727 
00728 
00729 
00730