Merging requests
[mldemos:mldemos.git] / _AlgorithmsPlugins / KernelMethods / kmeans.cpp
1 /*********************************************************************\r
2 MLDemos: A User-Friendly visualization toolkit for machine learning\r
3 Copyright (C) 2010  Basilio Noris\r
4 Contact: mldemos@b4silio.com\r
5 \r
6 This library is free software; you can redistribute it and/or\r
7 modify it under the terms of the GNU Lesser General Public\r
8 License as published by the Free Software Foundation; either\r
9 version 2.1 of the License, or (at your option) any later version.\r
10 \r
11 This library is distributed in the hope that it will be useful,\r
12 but WITHOUT ANY WARRANTY; without even the implied warranty of\r
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU\r
14 Library General Public License for more details.\r
15 \r
16 You should have received a copy of the GNU Lesser General Public\r
17 License along with this library; if not, write to the Free\r
18 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.\r
19 *********************************************************************/\r
20 #include <stdio.h>\r
21 #include <vector>\r
22 #include "public.h"\r
23 #include "basicMath.h"\r
24 #include <mymaths.h>\r
25 #include "kmeans.h"\r
26 #include <QTime>\r
27 #include <QDebug>\r
28 \r
29 using namespace std;\r
30 \r
31 \r
32 KMeansCluster::KMeansCluster(u32 cnt)\r
33     : clusters(cnt), sigma(NULL), pi(NULL), bGMM(false), bSoft(false), beta(1), dim(2), power(2), plusPlus(true)\r
34 {\r
35     InitClusters();\r
36 }\r
37 \r
38 KMeansCluster::~KMeansCluster()\r
39 {\r
40     Clear();\r
41 }\r
42 \r
43 float KMeansCluster::Distance( fvec a, fvec b )\r
44 {\r
45     fvec dif = a-b;\r
46     float d = 0;\r
47     if(power == 0) // infinite distance\r
48     {\r
49         FOR(i, dif.size()) d = max(d, abs(dif[i]));\r
50     }\r
51     else if(power == 1) // manhattan distance\r
52     {\r
53         FOR(i, dif.size()) d += abs(dif[i]);\r
54     }\r
55     else if(power == 2)\r
56     {\r
57         d = dif*dif;\r
58     }\r
59     else if(power > 2)\r
60     {\r
61         FOR(i, dif.size())\r
62         {\r
63             float p = abs(dif[i]);\r
64             float p2 = 1;\r
65             FOR(j, power) p2 *= p;\r
66             d += p2;\r
67         }\r
68     }\r
69     return d;\r
70 }\r
71 \r
72 float KMeansCluster::Distance2( fvec a, fvec b )\r
73 {\r
74     float d = 0;\r
75     if(power == 0) // infinite distance\r
76     {\r
77         d = max(abs(a[0]-b[0]), abs(a[1]-b[1]));\r
78     }\r
79     else if(power == 1) // manhattan distance\r
80     {\r
81         d = abs(a[0]-b[0]) + abs(a[1]-b[1]);\r
82     }\r
83     else if(power == 2)\r
84     {\r
85         float d0 = a[0]-b[0];\r
86         float d1 = a[1]-b[1];\r
87         d = d0*d0 + d1*d1;\r
88     }\r
89     else if(power > 2)\r
90     {\r
91         float d0 = abs(a[0]-b[0]);\r
92         float d1 = abs(a[1]-b[1]);\r
93         float p0 = 1;\r
94         float p1 = 1;\r
95         FOR(j, power)\r
96         {\r
97             p0 *= d0;\r
98             p1 *= d1;\r
99         }\r
100         d = p0 + p1;\r
101     }\r
102     return d;\r
103 }\r
104 \r
105 void KMeansCluster::Update(bool bFirstIteration)\r
106 {\r
107     bool bSuperposed = false;\r
108     FOR(i, clusters)\r
109     {\r
110         FOR(j,i)\r
111         {\r
112             if(means[i] == means[j]) // we have 2 superposed clusters\r
113             {\r
114                 // we replace it with a new one\r
115                 FOR(d,means[i].size()) means[i][d] = rand()/(float)RAND_MAX;\r
116                 //bSuperposed = true;\r
117                 break;\r
118             }\r
119         }\r
120     }\r
121     if(bSuperposed) InitClusters();\r
122 \r
123     if(bGMM) GMMClustering(points, means, sigma, pi, clusters, bFirstIteration);\r
124     else if (bSoft) SoftKmeansClustering(points, means, clusters, beta, bFirstIteration);\r
125     else if(!bFirstIteration) KmeansClustering(points, means, clusters);\r
126     FOR(i, clusters)\r
127     {\r
128         float mindist = 1;\r
129         u32 closest = 0;\r
130         FOR(p, points.size())\r
131         {\r
132             float d = (points[p].point - means[i])*(points[p].point - means[i]);\r
133             if (d < mindist)\r
134             {\r
135                 mindist = d;\r
136                 closest = p;\r
137             }\r
138         }\r
139         closestIndices[i] = closest;\r
140     }\r
141 }\r
142 \r
143 void KMeansCluster::AddPoint(fvec point)\r
144 {\r
145     if(point.size() != dim) dim = point.size();\r
146     ClusterPoint cpoint;\r
147     cpoint.point = point;\r
148     points.push_back(cpoint);\r
149 }\r
150 \r
151 void KMeansCluster::AddPoints(std::vector<fvec> points)\r
152 {\r
153     FOR(i, points.size()) AddPoint(points[i]);\r
154 }\r
155 \r
156 void KMeansCluster::Clear()\r
157 {\r
158     points.clear();\r
159 }\r
160 \r
161 void KMeansCluster::SetClusters(u32 clusters)\r
162 {\r
163     this->clusters = max((u32)0,clusters);\r
164     InitClusters();\r
165 }\r
166 \r
167 void KMeansCluster::InitClusters()\r
168 {\r
169     srand(QTime::currentTime().msec());\r
170 \r
171     KILL(pi);\r
172     if(sigma) FOR(i, clusters) KILL(sigma[i]);\r
173     KILL(sigma);\r
174     if(!clusters) return;\r
175     means.resize(clusters);\r
176     pi = new double[clusters];\r
177     sigma = new double *[clusters];\r
178     closestIndices.resize(clusters);\r
179     FOR(i,clusters){\r
180         means[i].resize(dim);\r
181         pi[i] = 1.f/clusters;\r
182         sigma[i] = new double[4];\r
183         sigma[i][0] = sigma[i][3] = 0.1;\r
184         sigma[i][1] = sigma[i][2] = 0.05;\r
185     }\r
186     if(!points.size()){\r
187         // no points, just choose random centers\r
188         FOR(i,clusters)\r
189         {\r
190             FOR(d,dim)\r
191             {\r
192                 means[i][d] = float(rand())/RAND_MAX;\r
193             }\r
194             closestIndices[i] = 0;\r
195         }\r
196     }\r
197     else if (plusPlus){\r
198         InitClustersPlusPlus();\r
199     }\r
200     else // choose a point at random\r
201     {\r
202         FOR(i,clusters)\r
203         {\r
204             int index = rand()%points.size();\r
205             means[i] = points[index].point;\r
206             closestIndices[i] = index;\r
207         }\r
208     }\r
209 }\r
210 \r
211 /** Use K-means++ to set the initial cluster centers, see\r
212 * <a href="http://en.wikipedia.org/wiki/K-means%2B%2B">K-means++ (wikipedia)</a>\r
213 * This code is based on Apache Commons Maths' KMeansPlusPlusClusterer.java\r
214 */\r
215 void KMeansCluster::InitClustersPlusPlus()\r
216 {\r
217     // Set the corresponding element in this array to indicate when points are no longer available.\r
218     bvec pointTaken(points.size(),false);\r
219 \r
220     // Choose first cluster center uniformly at random from among the data points.\r
221     int firstPointIndex = rand() % points.size(); // not uniform, but fair?\r
222     means[0] = points[firstPointIndex].point;\r
223     closestIndices[0] = firstPointIndex;\r
224     pointTaken[firstPointIndex] = true; // must mark it as taken\r
225     qDebug("first point at rand = %d[%f;%f]", firstPointIndex, points[firstPointIndex].point[0], points[firstPointIndex].point[1]);\r
226 \r
227     // Stores the squared minimum distance of each point to its nearest cluster center\r
228     fvec minDistSquared(points.size(), 0.0f);\r
229 \r
230     // Initialize the distances. Easy, since the only cluster is the first point\r
231     FOR(i, points.size())\r
232     {\r
233         if (i != firstPointIndex) // first point isn't considered\r
234         {\r
235             float d = Distance(points[firstPointIndex].point, points[i].point);\r
236             qDebug("initial minDistSquared to %i[%f;%f] = %f", i, points[i].point[0], points[i].point[1], d);\r
237             minDistSquared[i] = d * d;\r
238         }\r
239     }\r
240 \r
241     for(u32 centerCount = 1; centerCount < clusters; ++centerCount) // start at 1!\r
242     {\r
243         qDebug("picking cluster center %i----------------------------", centerCount);\r
244 \r
245         // Sum up the squared distances for the points not already taken.\r
246         float distSqSum = 0.0f;\r
247         FOR(j,points.size())\r
248         {\r
249             if (!pointTaken[j]) {\r
250                 distSqSum += minDistSquared[j];\r
251             }\r
252         }\r
253         qDebug(" distSqSum=%f",distSqSum);\r
254 \r
255         // Choose one new point at random as a new center, using a weighted\r
256         // probability distribution where a point x is chosen with probability proportional to D(x)^2\r
257         float r = (rand() / float(RAND_MAX)) * distSqSum;\r
258         qDebug(" random index %f",r);\r
259         // The index of the next point to be added to the resultSet.\r
260         bool nextPointFound= false;\r
261         u32 nextPointIndex = 0;\r
262 \r
263         // Sum through the squared min distances again, stopping when sum >= r.\r
264         float sum = 0.0f;\r
265         for(size_t j=0; j < points.size() && !nextPointFound; ++j)\r
266         {\r
267             if (!pointTaken[j])\r
268             {\r
269                 sum += minDistSquared[j];\r
270                 if (sum >= r)\r
271                 {\r
272                     nextPointIndex = (u32) j;\r
273                     nextPointFound = true;\r
274                 }\r
275             }\r
276         }\r
277 \r
278         // If no point was found in the previous for loop, probably because distances are extremely small.\r
279         // Just pick the last available point.\r
280         if (!nextPointFound)\r
281         {\r
282             qDebug("loop empty, pick one at rand");\r
283             for(size_t j=0; j < points.size() && !nextPointFound; ++j)\r
284             {\r
285                 if (!pointTaken[j])\r
286                 {\r
287                     nextPointIndex = j;\r
288                     nextPointFound = true;\r
289                 }\r
290             }\r
291         }\r
292 \r
293         // assert(nextPointFound); // if wanted\r
294 \r
295         // Set the new cluster point\r
296         means[centerCount] = points[nextPointIndex].point;\r
297         closestIndices[centerCount] = nextPointIndex;\r
298         pointTaken[nextPointIndex] = true;\r
299         qDebug(" new point %d[%f;%f]", nextPointIndex, points[nextPointIndex].point[0], points[nextPointIndex].point[1]);\r
300 \r
301         // Update minDistSquared. We only have to compute the distance to the new center, and update it if it is shorter\r
302         for(size_t j=0; j < points.size(); ++j)\r
303         {\r
304             if (!pointTaken[j])\r
305             {\r
306                 float d = Distance(points[nextPointIndex].point, points[j].point);\r
307                 float dSqr = d * d;\r
308                 if (dSqr < minDistSquared[j]) {\r
309                     minDistSquared[j] = dSqr;\r
310                 }\r
311             }\r
312         }\r
313     }\r
314 }\r
315 \r
316 \r
317 \r
318 ivec KMeansCluster::GetClosestPoints()\r
319 {\r
320     return closestIndices;\r
321 }\r
322 \r
323 inline float fastExp(const float x)\r
324 {\r
325     if(-x>90) return 0;\r
326     else return expf(x);\r
327 }\r
328 \r
329 void KMeansCluster::Test( fvec sample, fvec &res )\r
330 {\r
331     if(res.size() != clusters) res.resize(clusters);\r
332     if(bSoft)\r
333     {\r
334         fvec distances;\r
335         distances.resize(clusters);\r
336         // compute the distance to each clusters\r
337         float distanceSum = 0;\r
338         if(dim==2)\r
339         {\r
340             for (int j=0;j<clusters;j++){\r
341                 float d0 = means[j][0] - sample[0];\r
342                 float d1 = means[j][1] - sample[1];\r
343 \r
344                 distances[j] = fastExp(-beta * sqrtf(d0*d0 + d1*d1));\r
345                 distanceSum += distances[j];\r
346             }\r
347         }\r
348         else\r
349         {\r
350             for (int j=0;j<clusters;j++){\r
351                 distances[j] = fastExp(-beta * sqrtf((means[j] - sample)*(means[j] - sample)));\r
352                 distanceSum += distances[j];\r
353             }\r
354         }\r
355 \r
356         // compute the weights for each cluster\r
357         for (int j=0;j<clusters;j++){\r
358             res[j] = distances[j] / float(distanceSum);\r
359         }\r
360     }\r
361     else\r
362     {\r
363         FOR(d, res.size()) res[d] = 0;\r
364         int minIndex = 0;\r
365         float minDist = FLT_MAX;\r
366         if(dim==2)\r
367         {\r
368             FOR(j, clusters)\r
369             {\r
370                 float distance = Distance2(sample, means[j]);\r
371                 if(distance < minDist)\r
372                 {\r
373                     minIndex = j;\r
374                     minDist = distance;\r
375                 }\r
376             }\r
377         }\r
378         else\r
379         {\r
380             FOR(j, clusters)\r
381             {\r
382                 float distance = Distance(sample, means[j]);\r
383                 if(distance < minDist)\r
384                 {\r
385                     minIndex = j;\r
386                     minDist = distance;\r
387                 }\r
388             }\r
389         }\r
390         res[minIndex] = 1;\r
391     }\r
392 }\r
393 \r
394 \r
395 /**\r
396 * performs the K-mean clustering algorithm\r
397 *\r
398 * @param points[]  : each element of this array contains two elements\r
399 *                    point, a fvec with the coordinates of the points\r
400 *                    cluster, the number of the cluster the point belongs to\r
401 * @param nbPoints  : number of point in the array\r
402 * @param oldMeans  : array of fvec containing the old positions of the cluster centers\r
403 * @param nbCluster : number of clusters\r
404 * @param limits    : boundaries of the values for the points coordinates (default is the image size [320x240])\r
405 *\r
406 */\r
407 \r
408 void KMeansCluster::KmeansClustering(std::vector<ClusterPoint> &points, vector<fvec> &oldMeans, int nbClusters)\r
409 {\r
410     // check that we didnt try to use zero clusters\r
411     nbClusters = !nbClusters ? 1 : nbClusters;\r
412 \r
413     if((u32)nbClusters > points.size()) nbClusters = points.size();\r
414 \r
415     // new means (centers) for the clusters\r
416     vector<fvec>means;\r
417     means.resize(nbClusters);\r
418 \r
419     // used to check all the distances from the current point to each cluster\r
420     fvec distances;\r
421     distances.resize(nbClusters);\r
422 \r
423     means = oldMeans;\r
424 \r
425     int nbPoints = points.size();\r
426 \r
427     // did at least one point move from one cluster to another ?\r
428     bool bSomethingChanged = true;\r
429 \r
430     // the kmeans loop\r
431     //while(bSomethingChanged){\r
432         bSomethingChanged = false;\r
433 \r
434         //classify the points into clusters\r
435         for (register int i=0; i<nbPoints; i++){\r
436             // compute the distance to each clusters\r
437             for (register int j=0;j<nbClusters;j++){\r
438                 distances[j] = Distance(points[i].point, means[j]);\r
439             }\r
440 \r
441             // find the closest cluster\r
442             if (points[i].cluster != FindSmallest(distances)){\r
443                 bSomethingChanged = true;\r
444                 points[i].cluster = FindSmallest(distances);\r
445             }\r
446         }\r
447 \r
448         //compute the new means for each cluster\r
449         Mean(points, means, nbClusters);\r
450     //}\r
451 \r
452     oldMeans = means;\r
453 }\r
454 \r
455 \r
456 /**\r
457 * returns the square distance between two points\r
458 *\r
459 * @param p1 : the first point\r
460 * @param p2 : the second one\r
461 * @return   : the distance\r
462 *\r
463 *       note: if this function is more than a couple of lines long,\r
464 *             you're doing something wrong!\r
465 */\r
466 float SquareNorm(fvec p1, fvec p2)\r
467 {\r
468     return p1*p2;\r
469 }\r
470 \r
471 \r
472 /**\r
473 * returns the index of the smallest value in the array\r
474 *\r
475 * @param values   : the array containing the distances\r
476 * @param nbValues : the length of the array\r
477 * @return   : the index of the smallest distance in the array\r
478 *\r
479 */\r
480 int FindSmallest(fvec values)\r
481 {\r
482 \r
483     // initialize the minimum value and index to the first values in the array\r
484     int minIndex = 0;\r
485     float minValue = values[0];\r
486 \r
487     // go through the array and get the smallest value\r
488     for (int i=0; i<values.size(); i++){\r
489         if (values[i]<minValue){\r
490             minIndex = i;\r
491             minValue = values[i];\r
492         }\r
493     }\r
494 \r
495     return minIndex;\r
496 }\r
497 \r
498 /**\r
499 * computes the means for each cluster\r
500 *\r
501 * @param points     : the array containing the cluster points\r
502 * @param nbPoints   : the length of the points array\r
503 * @param means      : the array containing the means for the clusters, will be modified by the function\r
504 * @param nbClusters : the number of clusters (also = length of the means array)\r
505 *\r
506 */\r
507 void KMeansCluster::Mean(std::vector<ClusterPoint> &points, vector<fvec> &means, int nbClusters)\r
508 {\r
509     // counters to know how many points went into each cluster\r
510     int *nbPointInCluster = new int[nbClusters];\r
511 \r
512     // reinitialize the center for each cluster\r
513     for(int k = 0; k < nbClusters; k++){\r
514         FOR(d,dim) means[k][d] = 0;\r
515         nbPointInCluster[k] = 0;\r
516     }\r
517 \r
518     // sum the points\r
519     for (register int i = 0; (u32)i<points.size(); i++){\r
520         means[points[i].cluster] += points[i].point;\r
521         nbPointInCluster[points[i].cluster] += 1;\r
522     }\r
523 \r
524     // normalize by the number of points added to each cluster\r
525     for(int k = 0; k < nbClusters; k++){\r
526         if (nbPointInCluster[k]){\r
527             means[k] /= (float)nbPointInCluster[k];\r
528         }\r
529     }\r
530 \r
531     // delete the unused memory\r
532     delete [] nbPointInCluster;\r
533 }\r
534 \r
535 \r
536 // computes the means for each cluster\r
537 void KMeansCluster::SoftMean(std::vector<ClusterPoint> &points, vector<fvec> &means, int nbClusters)\r
538 {\r
539     // counters to know how many points went into each cluster\r
540     float *weightsOfPointsInCluster = new float[nbClusters];\r
541 \r
542     // reinitialize the center for each cluster\r
543     for(int k = 0; k < nbClusters; k++){\r
544         FOR(d,dim) means[k][d] = 0;\r
545         weightsOfPointsInCluster[k] = 0;\r
546     }\r
547 \r
548     // sum the points, for each cluster use the point's weight\r
549     for (unsigned int i = 0; i<points.size(); i++){\r
550         for(int k = 0; k<nbClusters; k++){\r
551             means[k] += points[i].point * points[i].weights[k];\r
552             weightsOfPointsInCluster[k] += points[i].weights[k];\r
553         }\r
554     }\r
555 \r
556     // normalize by the weights of the points added to each cluster\r
557     for(int k = 0; k < nbClusters; k++){\r
558         if (weightsOfPointsInCluster[k] != 0){\r
559             means[k] /= weightsOfPointsInCluster[k];\r
560         }\r
561     }\r
562 \r
563     // delete the unused memory\r
564     delete [] weightsOfPointsInCluster;\r
565 }\r
566 \r
567 \r
568 /**\r
569 * performs the Soft K-mean clustering algorithm\r
570 *\r
571 * @param points[]  : each element of this array contains two elements\r
572 *                    point, a fvec with the coordinates of the points\r
573 *                    weights, the weights of influence of the point on each cluster\r
574 * @param nbPoints  : number of point in the array\r
575 * @param oldMeans  : array of fvec containing the old positions of the cluster centers\r
576 * @param nbCluster : number of clusters\r
577 * @param limits    : boundaries of the values for the points coordinates (default is the image size [320x240])\r
578 * @param beta      : soft boundary stiffness (sigma = 1 / sqrt(beta))\r
579 *\r
580 */\r
581 \r
582 void KMeansCluster::SoftKmeansClustering(std::vector<ClusterPoint> &points, vector<fvec> &oldMeans, int nbClusters, float beta, bool bEStep)\r
583 {\r
584     // check that we didnt try to use zero clusters\r
585     nbClusters = !nbClusters ? 1 : nbClusters;\r
586 \r
587     if((unsigned int)nbClusters > points.size()) nbClusters = points.size();\r
588 \r
589     // contains the means for each cluster\r
590     vector<fvec>means;\r
591     means.resize(nbClusters);\r
592 \r
593     // used to check all the distances from the current point to each cluster\r
594     fvec distances;\r
595     distances.resize(nbClusters);\r
596 \r
597     // Random number generation for initial means of clusters\r
598     // initialize the random seed with the current cpu time\r
599     srand(QTime::currentTime().msec());\r
600 \r
601     means = oldMeans;\r
602 \r
603     int nbPoints = points.size();\r
604 \r
605     // initialize the points weights\r
606     for (int i=0; i < nbPoints; i++)\r
607     {\r
608         KILL(points[i].weights);\r
609         points[i].weights = new float[nbClusters];\r
610         FOR(j, nbClusters) points[i].weights[j] = 0;\r
611     }\r
612 \r
613     //classify the points into clusters\r
614     for (int i=0; i<nbPoints; i++){\r
615         fvec point = points[i].point;\r
616         // compute the distance to each clusters\r
617         float distanceSum = 0;\r
618         for (int j=0;j<nbClusters;j++){\r
619             distances[j] = expf(-beta * sqrtf((means[j] - point)*(means[j] - point)));\r
620             distanceSum += distances[j];\r
621         }\r
622 \r
623         // compute the weights for each cluster\r
624         for (int j=0;j<nbClusters;j++){\r
625             points[i].weights[j] = distances[j] / float(distanceSum);\r
626         }\r
627     }\r
628 \r
629     if(!bEStep)\r
630     {\r
631         //compute the new means for each cluster\r
632         SoftMean(points, means, nbClusters);\r
633     }\r
634 \r
635     oldMeans = means;\r
636 }\r
637 \r
638 void KMeansCluster::GMMClustering(std::vector<ClusterPoint> &points, vector<fvec> &oldMeans, double **oldSigma, double*oldPi, int nbClusters, bool bEStep)\r
639 {\r
640     // check that we didnt try to use zero clusters\r
641     nbClusters = !nbClusters ? 1 : nbClusters;\r
642 \r
643     if((unsigned int)nbClusters > points.size()) nbClusters = points.size();\r
644 \r
645     // contains the means for each cluster\r
646     vector<fvec>means;\r
647     means.resize(nbClusters);\r
648     double *pi = new double[nbClusters];\r
649     double **sigma = new double *[nbClusters];\r
650     FOR(i, nbClusters)\r
651     {\r
652         sigma[i] = new double[4];\r
653         sigma[i][0] = sigma[i][3] = 0.1f;\r
654         sigma[i][1] = sigma[i][2] = 0;\r
655     }\r
656 \r
657     // used to check all the distances from the current point to each cluster\r
658     double *distances = new double[nbClusters];\r
659 \r
660     // Random number generation for initial means of clusters\r
661     // initialize the random seed with the current cpu time\r
662     srand(QTime::currentTime().msec());\r
663 \r
664     // initialize the means as the old values\r
665     // divide by [320x240] to avoid numerical precision problems\r
666     // WARNING: from now on the values will go from 0 to 1\r
667     for (int i=0; i<nbClusters; i++){\r
668         means[i] = oldMeans[i];\r
669     }\r
670     for (int i=0; i<nbClusters; i++){\r
671         pi[i] = oldPi[i];\r
672     }\r
673     for (int i=0; i<nbClusters; i++){\r
674         for (int j=0; j<4; j++) sigma[i][j] = oldSigma[i][j];\r
675     }\r
676 \r
677     int nbPoints = points.size();\r
678 \r
679 \r
680     if(bEStep)\r
681     {\r
682         // initialize the points weights\r
683         u32 cnt = 0;\r
684         for (int i=0; i < nbPoints; i++)\r
685         {\r
686             KILL(points[i].weights);\r
687             points[i].weights = new float[nbClusters];\r
688             FOR(j, nbClusters) points[i].weights[j] = 0;\r
689             points[i].weights[cnt++%nbClusters] = 1.f;\r
690         }\r
691     }\r
692 \r
693     if(!bEStep)\r
694     {\r
695         if(points[0].weights == NULL)\r
696         {\r
697             // initialize the points weights\r
698             for (int i=0; i < nbPoints; i++)\r
699             {\r
700                 KILL(points[i].weights);\r
701                 points[i].weights = new float[nbClusters];\r
702                 FOR(j, nbClusters) points[i].weights[j] = 0;\r
703             }\r
704         }\r
705 \r
706         //classify the points into clusters\r
707         u32 flipper=0;\r
708         for (int i=0; i<nbPoints; i++){\r
709             fvec point = points[i].point;\r
710             // compute the distance to each clusters\r
711             double distanceSum = 0;\r
712             for (int j=0;j<nbClusters;j++){\r
713                 fvec a = (point - means[j]);\r
714                 double *s = sigma[j];\r
715                 double sdet = s[0]*s[3] - s[1]*s[2];\r
716                 double sinv[4] = {s[3], -s[1], -s[2], s[0]};\r
717                 FOR(k,4) sinv[k] /= sdet;\r
718 \r
719                 double b = a[0]*a[0]*sinv[0] + a[0]*a[1]*(sinv[1]+sinv[2]) + a[1]*a[1]*sinv[3];\r
720                 b *= -0.5;\r
721                 double dist = exp(b);\r
722                 dist /= sqrt(sdet);\r
723                 distances[j] = pi[j]*dist;\r
724                 distances[j] /= 2*(double)PIf;\r
725                 distanceSum += distances[j];\r
726             }\r
727 \r
728             // compute the weights for each cluster\r
729             if(distanceSum != distanceSum)\r
730             {\r
731                 FOR(j, nbClusters) points[i].weights[j] = 0;\r
732                 points[i].weights[flipper++%nbClusters] = 1;\r
733             }\r
734             else\r
735                 for (int j=0;j<nbClusters;j++){\r
736                     points[i].weights[j] = (float)(distances[j] / distanceSum);\r
737                 }\r
738         }\r
739     }\r
740 \r
741     //compute the new means for each cluster\r
742     for (int i=0; i<nbClusters; i++)\r
743     {\r
744         fvec mean;\r
745         mean.resize(2,0);\r
746         float respTotal = 0;\r
747         for (int j=0; j<nbPoints; j++)\r
748         {\r
749             mean += points[j].point * points[j].weights[i];\r
750             respTotal += points[j].weights[i];\r
751         }\r
752         means[i] = mean / respTotal;\r
753     }\r
754 \r
755     //compute the new prior for each cluster\r
756     float respTotal = 0;\r
757     for (int i=0; i<nbClusters; i++)\r
758     {\r
759         float resp = 0;\r
760         for (int j=0; j<nbPoints; j++)\r
761         {\r
762             resp += points[j].weights[i];\r
763         }\r
764         respTotal += resp;\r
765         pi[i] = resp;\r
766     }\r
767     for (int i=0; i<nbClusters; i++)\r
768     {\r
769         pi[i] /= respTotal;\r
770     }\r
771 \r
772     //compute the new sigma for each cluster\r
773     for (int i=0; i<nbClusters; i++)\r
774     {\r
775         float sums[3];\r
776         float resp = 0;\r
777         for (int j=0; j<3; j++) sums[j] = 0;\r
778         for (int j=0; j<nbPoints; j++)\r
779         {\r
780             float r = points[j].weights[i];\r
781             if(r==0) continue;\r
782             fvec diff = points[j].point - means[i];\r
783             sums[0] += r*(diff[0]*diff[0]);\r
784             sums[1] += r*(diff[0]*diff[1]);\r
785             sums[2] += r*(diff[1]*diff[1]);\r
786             resp += r;\r
787         }\r
788         for (int j=0; j<3; j++) sums[j] /= resp;\r
789         sigma[i][0] = sums[0];\r
790         sigma[i][1] = sigma[i][2] = sums[1];\r
791         sigma[i][3] = sums[2];\r
792     }\r
793 \r
794     // we copy the new values into the old ones\r
795     for (int i=0; i<nbClusters; i++){\r
796         oldMeans[i] = means[i];\r
797     }\r
798     for (int i=0; i<nbClusters; i++){\r
799         oldPi[i] = pi[i];\r
800     }\r
801     for (int i=0; i<nbClusters; i++){\r
802         for (int j=0; j<4; j++) oldSigma[i][j] = sigma[i][j];\r
803     }\r
804 \r
805     // we delete the memory we used\r
806     delete [] distances;\r
807     delete [] pi;\r
808     for (int i=0; i<nbClusters; i++) delete [] sigma[i];\r
809     delete [] sigma;\r
810 }\r