00310 {
00311 
00312   new_cluster->next_cluster = head_cluster;
00313 
00314   find_nearest_cluster (new_cluster, head_cluster);
00315 
00316 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

00097 {
00098   while (voxel_ptr != NULL)
00099     {
00100       print_voxel (voxel_ptr);
00101       voxel_ptr = voxel_ptr->next_voxel;
00102     }
00103   printf ("\n");
00104 }

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 }

00085 {
00086   if (voxel_ptr != NULL)
00087     printf ("%d ", voxel_ptr->index);
00088 }

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 }

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 }

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 }

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 }