Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
StatClust.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
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
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
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
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
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
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
00150
00151
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
00165
00166
00167 vector_destroy (&v);
00168 vector_destroy (&sv);
00169 }
00170
00171
00172
00173
00174
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
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
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
00238
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
00262
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
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
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
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
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
00367
00368
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
00386
00387
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
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
00477
00478
00479 cluster * consolidate_clusters (cluster * aclust, cluster * bclust,
00480 cluster * head_clust)
00481 {
00482 cluster * abclust = NULL;
00483
00484
00485
00486 abclust = merge_clusters (aclust, bclust);
00487
00488
00489
00490 head_clust = remove_clusters (aclust, bclust, head_clust);
00491
00492
00493
00494 add_cluster (abclust, head_clust);
00495
00496
00497
00498 return (abclust);
00499 }
00500
00501
00502
00503
00504
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
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
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
00553 head_clust = consolidate_clusters (aclust, bclust, head_clust);
00554
00555
00556 return (head_clust);
00557 }
00558
00559
00560
00561
00562
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
00577 guard = initialize_cluster();
00578 guard->next_cluster = head_clust;
00579 ip = guard;
00580
00581 while (ip->next_cluster != NULL)
00582 {
00583
00584 i = ip->next_cluster;
00585 mp = ip;
00586 m = i;
00587 jp = i;
00588
00589
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
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
00611 ip = i;
00612
00613 }
00614
00615
00616
00617 head_clust = guard->next_cluster;
00618 delete_cluster (guard);
00619
00620 return (head_clust);
00621 }
00622
00623
00624
00625
00626
00627
00628
00629
00630 void calc_covariance
00631 (
00632 matrix * s,
00633 matrix * sinv
00634 )
00635
00636 {
00637 int ivox;
00638 int ip, jp;
00639 int ok;
00640
00641 vector mean;
00642 matrix covar;
00643 matrix cinv;
00644
00645 char message[80];
00646
00647
00648
00649 vector_initialize (&mean);
00650 matrix_initialize (&covar);
00651 matrix_initialize (&cinv);
00652
00653
00654
00655 vector_create (SC_dimension, &mean);
00656 matrix_create (SC_dimension, SC_dimension, &covar);
00657
00658
00659
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
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
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
00692
00693
00694
00695
00696 ok = matrix_inverse (covar, &cinv);
00697 if (! ok)
00698 SC_error
00699 ("Unable to calculate inverse of covariance matrix");
00700
00701
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
00708 ok = matrix_inverse (*sinv, s);
00709 if (! ok)
00710 SC_error
00711 ("Unable to calculate square root of covariance matrix");
00712
00713
00714
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