1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install,
7 // copy or use the software.
8 //
9 //
10 // Intel License Agreement
11 // For Open Source Computer Vision Library
12 //
13 // Copyright( C) 2000, Intel Corporation, all rights reserved.
14 // Third party copyrights are property of their respective owners.
15 //
16 // Redistribution and use in source and binary forms, with or without modification,
17 // are permitted provided that the following conditions are met:
18 //
19 // * Redistribution's of source code must retain the above copyright notice,
20 // this list of conditions and the following disclaimer.
21 //
22 // * Redistribution's in binary form must reproduce the above copyright notice,
23 // this list of conditions and the following disclaimer in the documentation
24 // and/or other materials provided with the distribution.
25 //
26 // * The name of Intel Corporation may not be used to endorse or promote products
27 // derived from this software without specific prior written permission.
28 //
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the Intel Corporation or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 //(including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort(including negligence or otherwise) arising in any way out of
38 // the use of this software, even ifadvised of the possibility of such damage.
39 //
40 //M*/
41
42 #include "_ml.h"
43
44
45 /*
46 CvEM:
47 * params.nclusters - number of clusters to cluster samples to.
48 * means - calculated by the EM algorithm set of gaussians' means.
49 * log_weight_div_det - auxilary vector that k-th component is equal to
50 (-2)*ln(weights_k/det(Sigma_k)^0.5),
51 where <weights_k> is the weight,
52 <Sigma_k> is the covariation matrice of k-th cluster.
53 * inv_eigen_values - set of 1*dims matrices, <inv_eigen_values>[k] contains
54 inversed eigen values of covariation matrice of the k-th cluster.
55 In the case of <cov_mat_type> == COV_MAT_DIAGONAL,
56 inv_eigen_values[k] = Sigma_k^(-1).
57 * covs_rotate_mats - used only if cov_mat_type == COV_MAT_GENERIC, in all the
58 other cases it is NULL. <covs_rotate_mats>[k] is the orthogonal
59 matrice, obtained by the SVD-decomposition of Sigma_k.
60 Both <inv_eigen_values> and <covs_rotate_mats> fields are used for representation of
61 covariation matrices and simplifying EM calculations.
62 For fixed k denote
63 u = covs_rotate_mats[k],
64 v = inv_eigen_values[k],
65 w = v^(-1);
66 if <cov_mat_type> == COV_MAT_GENERIC, then Sigma_k = u w u',
67 else Sigma_k = w.
68 Symbol ' means transposition.
69 */
70
71
CvEM()72 CvEM::CvEM()
73 {
74 means = weights = probs = inv_eigen_values = log_weight_div_det = 0;
75 covs = cov_rotate_mats = 0;
76 }
77
CvEM(const CvMat * samples,const CvMat * sample_idx,CvEMParams params,CvMat * labels)78 CvEM::CvEM( const CvMat* samples, const CvMat* sample_idx,
79 CvEMParams params, CvMat* labels )
80 {
81 means = weights = probs = inv_eigen_values = log_weight_div_det = 0;
82 covs = cov_rotate_mats = 0;
83
84 // just invoke the train() method
85 train(samples, sample_idx, params, labels);
86 }
87
~CvEM()88 CvEM::~CvEM()
89 {
90 clear();
91 }
92
93
clear()94 void CvEM::clear()
95 {
96 int i;
97
98 cvReleaseMat( &means );
99 cvReleaseMat( &weights );
100 cvReleaseMat( &probs );
101 cvReleaseMat( &inv_eigen_values );
102 cvReleaseMat( &log_weight_div_det );
103
104 if( covs || cov_rotate_mats )
105 {
106 for( i = 0; i < params.nclusters; i++ )
107 {
108 if( covs )
109 cvReleaseMat( &covs[i] );
110 if( cov_rotate_mats )
111 cvReleaseMat( &cov_rotate_mats[i] );
112 }
113 cvFree( &covs );
114 cvFree( &cov_rotate_mats );
115 }
116 }
117
118
set_params(const CvEMParams & _params,const CvVectors & train_data)119 void CvEM::set_params( const CvEMParams& _params, const CvVectors& train_data )
120 {
121 CV_FUNCNAME( "CvEM::set_params" );
122
123 __BEGIN__;
124
125 int k;
126
127 params = _params;
128 params.term_crit = cvCheckTermCriteria( params.term_crit, 1e-6, 10000 );
129
130 if( params.cov_mat_type != COV_MAT_SPHERICAL &&
131 params.cov_mat_type != COV_MAT_DIAGONAL &&
132 params.cov_mat_type != COV_MAT_GENERIC )
133 CV_ERROR( CV_StsBadArg, "Unknown covariation matrix type" );
134
135 switch( params.start_step )
136 {
137 case START_M_STEP:
138 if( !params.probs )
139 CV_ERROR( CV_StsNullPtr, "Probabilities must be specified when EM algorithm starts with M-step" );
140 break;
141 case START_E_STEP:
142 if( !params.means )
143 CV_ERROR( CV_StsNullPtr, "Mean's must be specified when EM algorithm starts with E-step" );
144 break;
145 case START_AUTO_STEP:
146 break;
147 default:
148 CV_ERROR( CV_StsBadArg, "Unknown start_step" );
149 }
150
151 if( params.nclusters < 1 )
152 CV_ERROR( CV_StsOutOfRange, "The number of clusters (mixtures) should be > 0" );
153
154 if( params.probs )
155 {
156 const CvMat* p = params.weights;
157 if( !CV_IS_MAT(p) ||
158 CV_MAT_TYPE(p->type) != CV_32FC1 &&
159 CV_MAT_TYPE(p->type) != CV_64FC1 ||
160 p->rows != train_data.count ||
161 p->cols != params.nclusters )
162 CV_ERROR( CV_StsBadArg, "The array of probabilities must be a valid "
163 "floating-point matrix (CvMat) of 'nsamples' x 'nclusters' size" );
164 }
165
166 if( params.means )
167 {
168 const CvMat* m = params.means;
169 if( !CV_IS_MAT(m) ||
170 CV_MAT_TYPE(m->type) != CV_32FC1 &&
171 CV_MAT_TYPE(m->type) != CV_64FC1 ||
172 m->rows != params.nclusters ||
173 m->cols != train_data.dims )
174 CV_ERROR( CV_StsBadArg, "The array of mean's must be a valid "
175 "floating-point matrix (CvMat) of 'nsamples' x 'dims' size" );
176 }
177
178 if( params.weights )
179 {
180 const CvMat* w = params.weights;
181 if( !CV_IS_MAT(w) ||
182 CV_MAT_TYPE(w->type) != CV_32FC1 &&
183 CV_MAT_TYPE(w->type) != CV_64FC1 ||
184 w->rows != 1 && w->cols != 1 ||
185 w->rows + w->cols - 1 != params.nclusters )
186 CV_ERROR( CV_StsBadArg, "The array of weights must be a valid "
187 "1d floating-point vector (CvMat) of 'nclusters' elements" );
188 }
189
190 if( params.covs )
191 for( k = 0; k < params.nclusters; k++ )
192 {
193 const CvMat* cov = params.covs[k];
194 if( !CV_IS_MAT(cov) ||
195 CV_MAT_TYPE(cov->type) != CV_32FC1 &&
196 CV_MAT_TYPE(cov->type) != CV_64FC1 ||
197 cov->rows != cov->cols || cov->cols != train_data.dims )
198 CV_ERROR( CV_StsBadArg,
199 "Each of covariation matrices must be a valid square "
200 "floating-point matrix (CvMat) of 'dims' x 'dims'" );
201 }
202
203 __END__;
204 }
205
206
207 /****************************************************************************************/
208 float
predict(const CvMat * _sample,CvMat * _probs) const209 CvEM::predict( const CvMat* _sample, CvMat* _probs ) const
210 {
211 float* sample_data = 0;
212 void* buffer = 0;
213 int allocated_buffer = 0;
214 int cls = 0;
215
216 CV_FUNCNAME( "CvEM::predict" );
217 __BEGIN__;
218
219 int i, k, dims;
220 int nclusters;
221 int cov_mat_type = params.cov_mat_type;
222 double opt = FLT_MAX;
223 size_t size;
224 CvMat diff, expo;
225
226 dims = means->cols;
227 nclusters = params.nclusters;
228
229 CV_CALL( cvPreparePredictData( _sample, dims, 0, params.nclusters, _probs, &sample_data ));
230
231 // allocate memory and initializing headers for calculating
232 size = sizeof(double) * (nclusters + dims);
233 if( size <= CV_MAX_LOCAL_SIZE )
234 buffer = cvStackAlloc( size );
235 else
236 {
237 CV_CALL( buffer = cvAlloc( size ));
238 allocated_buffer = 1;
239 }
240 expo = cvMat( 1, nclusters, CV_64FC1, buffer );
241 diff = cvMat( 1, dims, CV_64FC1, (double*)buffer + nclusters );
242
243 // calculate the probabilities
244 for( k = 0; k < nclusters; k++ )
245 {
246 const double* mean_k = (const double*)(means->data.ptr + means->step*k);
247 const double* w = (const double*)(inv_eigen_values->data.ptr + inv_eigen_values->step*k);
248 double cur = log_weight_div_det->data.db[k];
249 CvMat* u = cov_rotate_mats[k];
250 // cov = u w u' --> cov^(-1) = u w^(-1) u'
251 if( cov_mat_type == COV_MAT_SPHERICAL )
252 {
253 double w0 = w[0];
254 for( i = 0; i < dims; i++ )
255 {
256 double val = sample_data[i] - mean_k[i];
257 cur += val*val*w0;
258 }
259 }
260 else
261 {
262 for( i = 0; i < dims; i++ )
263 diff.data.db[i] = sample_data[i] - mean_k[i];
264 if( cov_mat_type == COV_MAT_GENERIC )
265 cvGEMM( &diff, u, 1, 0, 0, &diff, CV_GEMM_B_T );
266 for( i = 0; i < dims; i++ )
267 {
268 double val = diff.data.db[i];
269 cur += val*val*w[i];
270 }
271 }
272
273 expo.data.db[k] = cur;
274 if( cur < opt )
275 {
276 cls = k;
277 opt = cur;
278 }
279 /* probability = (2*pi)^(-dims/2)*exp( -0.5 * cur ) */
280 }
281
282 if( _probs )
283 {
284 CV_CALL( cvConvertScale( &expo, &expo, -0.5 ));
285 CV_CALL( cvExp( &expo, &expo ));
286 if( _probs->cols == 1 )
287 CV_CALL( cvReshape( &expo, &expo, 0, nclusters ));
288 CV_CALL( cvConvertScale( &expo, _probs, 1./cvSum( &expo ).val[0] ));
289 }
290
291 __END__;
292
293 if( sample_data != _sample->data.fl )
294 cvFree( &sample_data );
295 if( allocated_buffer )
296 cvFree( &buffer );
297
298 return (float)cls;
299 }
300
301
302
train(const CvMat * _samples,const CvMat * _sample_idx,CvEMParams _params,CvMat * labels)303 bool CvEM::train( const CvMat* _samples, const CvMat* _sample_idx,
304 CvEMParams _params, CvMat* labels )
305 {
306 bool result = false;
307 CvVectors train_data;
308 CvMat* sample_idx = 0;
309
310 train_data.data.fl = 0;
311 train_data.count = 0;
312
313 CV_FUNCNAME("cvEM");
314
315 __BEGIN__;
316
317 int i, nsamples, nclusters, dims;
318
319 clear();
320
321 CV_CALL( cvPrepareTrainData( "cvEM",
322 _samples, CV_ROW_SAMPLE, 0, CV_VAR_CATEGORICAL,
323 0, _sample_idx, false, (const float***)&train_data.data.fl,
324 &train_data.count, &train_data.dims, &train_data.dims,
325 0, 0, 0, &sample_idx ));
326
327 CV_CALL( set_params( _params, train_data ));
328 nsamples = train_data.count;
329 nclusters = params.nclusters;
330 dims = train_data.dims;
331
332 if( labels && (!CV_IS_MAT(labels) || CV_MAT_TYPE(labels->type) != CV_32SC1 ||
333 labels->cols != 1 && labels->rows != 1 || labels->cols + labels->rows - 1 != nsamples ))
334 CV_ERROR( CV_StsBadArg,
335 "labels array (when passed) must be a valid 1d integer vector of <sample_count> elements" );
336
337 if( nsamples <= nclusters )
338 CV_ERROR( CV_StsOutOfRange,
339 "The number of samples should be greater than the number of clusters" );
340
341 CV_CALL( log_weight_div_det = cvCreateMat( 1, nclusters, CV_64FC1 ));
342 CV_CALL( probs = cvCreateMat( nsamples, nclusters, CV_64FC1 ));
343 CV_CALL( means = cvCreateMat( nclusters, dims, CV_64FC1 ));
344 CV_CALL( weights = cvCreateMat( 1, nclusters, CV_64FC1 ));
345 CV_CALL( inv_eigen_values = cvCreateMat( nclusters,
346 params.cov_mat_type == COV_MAT_SPHERICAL ? 1 : dims, CV_64FC1 ));
347 CV_CALL( covs = (CvMat**)cvAlloc( nclusters * sizeof(*covs) ));
348 CV_CALL( cov_rotate_mats = (CvMat**)cvAlloc( nclusters * sizeof(cov_rotate_mats[0]) ));
349
350 for( i = 0; i < nclusters; i++ )
351 {
352 CV_CALL( covs[i] = cvCreateMat( dims, dims, CV_64FC1 ));
353 CV_CALL( cov_rotate_mats[i] = cvCreateMat( dims, dims, CV_64FC1 ));
354 cvZero( cov_rotate_mats[i] );
355 }
356
357 init_em( train_data );
358 log_likelihood = run_em( train_data );
359 if( log_likelihood <= -DBL_MAX/10000. )
360 EXIT;
361
362 if( labels )
363 {
364 if( nclusters == 1 )
365 cvZero( labels );
366 else
367 {
368 CvMat sample = cvMat( 1, dims, CV_32F );
369 CvMat prob = cvMat( 1, nclusters, CV_64F );
370 int lstep = CV_IS_MAT_CONT(labels->type) ? 1 : labels->step/sizeof(int);
371
372 for( i = 0; i < nsamples; i++ )
373 {
374 int idx = sample_idx ? sample_idx->data.i[i] : i;
375 sample.data.ptr = _samples->data.ptr + _samples->step*idx;
376 prob.data.ptr = probs->data.ptr + probs->step*i;
377
378 labels->data.i[i*lstep] = cvRound(predict(&sample, &prob));
379 }
380 }
381 }
382
383 result = true;
384
385 __END__;
386
387 if( sample_idx != _sample_idx )
388 cvReleaseMat( &sample_idx );
389
390 cvFree( &train_data.data.ptr );
391
392 return result;
393 }
394
395
init_em(const CvVectors & train_data)396 void CvEM::init_em( const CvVectors& train_data )
397 {
398 CvMat *w = 0, *u = 0, *tcov = 0;
399
400 CV_FUNCNAME( "CvEM::init_em" );
401
402 __BEGIN__;
403
404 double maxval = 0;
405 int i, force_symm_plus = 0;
406 int nclusters = params.nclusters, nsamples = train_data.count, dims = train_data.dims;
407
408 if( params.start_step == START_AUTO_STEP || nclusters == 1 || nclusters == nsamples )
409 init_auto( train_data );
410 else if( params.start_step == START_M_STEP )
411 {
412 for( i = 0; i < nsamples; i++ )
413 {
414 CvMat prob;
415 cvGetRow( params.probs, &prob, i );
416 cvMaxS( &prob, 0., &prob );
417 cvMinMaxLoc( &prob, 0, &maxval );
418 if( maxval < FLT_EPSILON )
419 cvSet( &prob, cvScalar(1./nclusters) );
420 else
421 cvNormalize( &prob, &prob, 1., 0, CV_L1 );
422 }
423 EXIT; // do not preprocess covariation matrices,
424 // as in this case they are initialized at the first iteration of EM
425 }
426 else
427 {
428 CV_ASSERT( params.start_step == START_E_STEP && params.means );
429 if( params.weights && params.covs )
430 {
431 cvConvert( params.means, means );
432 cvReshape( weights, weights, 1, params.weights->rows );
433 cvConvert( params.weights, weights );
434 cvReshape( weights, weights, 1, 1 );
435 cvMaxS( weights, 0., weights );
436 cvMinMaxLoc( weights, 0, &maxval );
437 if( maxval < FLT_EPSILON )
438 cvSet( &weights, cvScalar(1./nclusters) );
439 cvNormalize( weights, weights, 1., 0, CV_L1 );
440 for( i = 0; i < nclusters; i++ )
441 CV_CALL( cvConvert( params.covs[i], covs[i] ));
442 force_symm_plus = 1;
443 }
444 else
445 init_auto( train_data );
446 }
447
448 CV_CALL( tcov = cvCreateMat( dims, dims, CV_64FC1 ));
449 CV_CALL( w = cvCreateMat( dims, dims, CV_64FC1 ));
450 if( params.cov_mat_type == COV_MAT_GENERIC )
451 CV_CALL( u = cvCreateMat( dims, dims, CV_64FC1 ));
452
453 for( i = 0; i < nclusters; i++ )
454 {
455 if( force_symm_plus )
456 {
457 cvTranspose( covs[i], tcov );
458 cvAddWeighted( covs[i], 0.5, tcov, 0.5, 0, tcov );
459 }
460 else
461 cvCopy( covs[i], tcov );
462 cvSVD( tcov, w, u, 0, CV_SVD_MODIFY_A + CV_SVD_U_T + CV_SVD_V_T );
463 if( params.cov_mat_type == COV_MAT_SPHERICAL )
464 cvSetIdentity( covs[i], cvScalar(cvTrace(w).val[0]/dims) );
465 else if( params.cov_mat_type == COV_MAT_DIAGONAL )
466 cvCopy( w, covs[i] );
467 else
468 {
469 // generic case: covs[i] = (u')'*max(w,0)*u'
470 cvGEMM( u, w, 1, 0, 0, tcov, CV_GEMM_A_T );
471 cvGEMM( tcov, u, 1, 0, 0, covs[i], 0 );
472 }
473 }
474
475 __END__;
476
477 cvReleaseMat( &w );
478 cvReleaseMat( &u );
479 cvReleaseMat( &tcov );
480 }
481
482
init_auto(const CvVectors & train_data)483 void CvEM::init_auto( const CvVectors& train_data )
484 {
485 CvMat* hdr = 0;
486 const void** vec = 0;
487 CvMat* class_ranges = 0;
488 CvMat* labels = 0;
489
490 CV_FUNCNAME( "CvEM::init_auto" );
491
492 __BEGIN__;
493
494 int nclusters = params.nclusters, nsamples = train_data.count, dims = train_data.dims;
495 int i, j;
496
497 if( nclusters == nsamples )
498 {
499 CvMat src = cvMat( 1, dims, CV_32F );
500 CvMat dst = cvMat( 1, dims, CV_64F );
501 for( i = 0; i < nsamples; i++ )
502 {
503 src.data.ptr = train_data.data.ptr[i];
504 dst.data.ptr = means->data.ptr + means->step*i;
505 cvConvert( &src, &dst );
506 cvZero( covs[i] );
507 cvSetIdentity( cov_rotate_mats[i] );
508 }
509 cvSetIdentity( probs );
510 cvSet( weights, cvScalar(1./nclusters) );
511 }
512 else
513 {
514 int max_count = 0;
515
516 CV_CALL( class_ranges = cvCreateMat( 1, nclusters+1, CV_32SC1 ));
517 if( nclusters > 1 )
518 {
519 CV_CALL( labels = cvCreateMat( 1, nsamples, CV_32SC1 ));
520 kmeans( train_data, nclusters, labels, cvTermCriteria( CV_TERMCRIT_ITER,
521 params.means ? 1 : 10, 0.5 ), params.means );
522 CV_CALL( cvSortSamplesByClasses( (const float**)train_data.data.fl,
523 labels, class_ranges->data.i ));
524 }
525 else
526 {
527 class_ranges->data.i[0] = 0;
528 class_ranges->data.i[1] = nsamples;
529 }
530
531 for( i = 0; i < nclusters; i++ )
532 {
533 int left = class_ranges->data.i[i], right = class_ranges->data.i[i+1];
534 max_count = MAX( max_count, right - left );
535 }
536 CV_CALL( hdr = (CvMat*)cvAlloc( max_count*sizeof(hdr[0]) ));
537 CV_CALL( vec = (const void**)cvAlloc( max_count*sizeof(vec[0]) ));
538 hdr[0] = cvMat( 1, dims, CV_32F );
539 for( i = 0; i < max_count; i++ )
540 {
541 vec[i] = hdr + i;
542 hdr[i] = hdr[0];
543 }
544
545 for( i = 0; i < nclusters; i++ )
546 {
547 int left = class_ranges->data.i[i], right = class_ranges->data.i[i+1];
548 int cluster_size = right - left;
549 CvMat avg;
550
551 if( cluster_size <= 0 )
552 continue;
553
554 for( j = left; j < right; j++ )
555 hdr[j - left].data.fl = train_data.data.fl[j];
556
557 CV_CALL( cvGetRow( means, &avg, i ));
558 CV_CALL( cvCalcCovarMatrix( vec, cluster_size, covs[i],
559 &avg, CV_COVAR_NORMAL | CV_COVAR_SCALE ));
560 weights->data.db[i] = (double)cluster_size/(double)nsamples;
561 }
562 }
563
564 __END__;
565
566 cvReleaseMat( &class_ranges );
567 cvReleaseMat( &labels );
568 cvFree( &hdr );
569 cvFree( &vec );
570 }
571
572
kmeans(const CvVectors & train_data,int nclusters,CvMat * labels,CvTermCriteria termcrit,const CvMat * centers0)573 void CvEM::kmeans( const CvVectors& train_data, int nclusters, CvMat* labels,
574 CvTermCriteria termcrit, const CvMat* centers0 )
575 {
576 CvMat* centers = 0;
577 CvMat* old_centers = 0;
578 CvMat* counters = 0;
579
580 CV_FUNCNAME( "CvEM::kmeans" );
581
582 __BEGIN__;
583
584 CvRNG rng = cvRNG(-1);
585 int i, j, k, nsamples, dims;
586 int iter = 0;
587 double max_dist = DBL_MAX;
588
589 termcrit = cvCheckTermCriteria( termcrit, 1e-6, 100 );
590 termcrit.epsilon *= termcrit.epsilon;
591 nsamples = train_data.count;
592 dims = train_data.dims;
593 nclusters = MIN( nclusters, nsamples );
594
595 CV_CALL( centers = cvCreateMat( nclusters, dims, CV_64FC1 ));
596 CV_CALL( old_centers = cvCreateMat( nclusters, dims, CV_64FC1 ));
597 CV_CALL( counters = cvCreateMat( 1, nclusters, CV_32SC1 ));
598 cvZero( old_centers );
599
600 if( centers0 )
601 {
602 CV_CALL( cvConvert( centers0, centers ));
603 }
604 else
605 {
606 for( i = 0; i < nsamples; i++ )
607 labels->data.i[i] = i*nclusters/nsamples;
608 cvRandShuffle( labels, &rng );
609 }
610
611 for( ;; )
612 {
613 CvMat* temp;
614
615 if( iter > 0 || centers0 )
616 {
617 for( i = 0; i < nsamples; i++ )
618 {
619 const float* s = train_data.data.fl[i];
620 int k_best = 0;
621 double min_dist = DBL_MAX;
622
623 for( k = 0; k < nclusters; k++ )
624 {
625 const double* c = (double*)(centers->data.ptr + k*centers->step);
626 double dist = 0;
627
628 for( j = 0; j <= dims - 4; j += 4 )
629 {
630 double t0 = c[j] - s[j];
631 double t1 = c[j+1] - s[j+1];
632 dist += t0*t0 + t1*t1;
633 t0 = c[j+2] - s[j+2];
634 t1 = c[j+3] - s[j+3];
635 dist += t0*t0 + t1*t1;
636 }
637
638 for( ; j < dims; j++ )
639 {
640 double t = c[j] - s[j];
641 dist += t*t;
642 }
643
644 if( min_dist > dist )
645 {
646 min_dist = dist;
647 k_best = k;
648 }
649 }
650
651 labels->data.i[i] = k_best;
652 }
653 }
654
655 if( ++iter > termcrit.max_iter )
656 break;
657
658 CV_SWAP( centers, old_centers, temp );
659 cvZero( centers );
660 cvZero( counters );
661
662 // update centers
663 for( i = 0; i < nsamples; i++ )
664 {
665 const float* s = train_data.data.fl[i];
666 k = labels->data.i[i];
667 double* c = (double*)(centers->data.ptr + k*centers->step);
668
669 for( j = 0; j <= dims - 4; j += 4 )
670 {
671 double t0 = c[j] + s[j];
672 double t1 = c[j+1] + s[j+1];
673
674 c[j] = t0;
675 c[j+1] = t1;
676
677 t0 = c[j+2] + s[j+2];
678 t1 = c[j+3] + s[j+3];
679
680 c[j+2] = t0;
681 c[j+3] = t1;
682 }
683 for( ; j < dims; j++ )
684 c[j] += s[j];
685 counters->data.i[k]++;
686 }
687
688 if( iter > 1 )
689 max_dist = 0;
690
691 for( k = 0; k < nclusters; k++ )
692 {
693 double* c = (double*)(centers->data.ptr + k*centers->step);
694 if( counters->data.i[k] != 0 )
695 {
696 double scale = 1./counters->data.i[k];
697 for( j = 0; j < dims; j++ )
698 c[j] *= scale;
699 }
700 else
701 {
702 const float* s;
703 for( j = 0; j < 10; j++ )
704 {
705 i = cvRandInt( &rng ) % nsamples;
706 if( counters->data.i[labels->data.i[i]] > 1 )
707 break;
708 }
709 s = train_data.data.fl[i];
710 for( j = 0; j < dims; j++ )
711 c[j] = s[j];
712 }
713
714 if( iter > 1 )
715 {
716 double dist = 0;
717 const double* c_o = (double*)(old_centers->data.ptr + k*old_centers->step);
718 for( j = 0; j < dims; j++ )
719 {
720 double t = c[j] - c_o[j];
721 dist += t*t;
722 }
723 if( max_dist < dist )
724 max_dist = dist;
725 }
726 }
727
728 if( max_dist < termcrit.epsilon )
729 break;
730 }
731
732 cvZero( counters );
733 for( i = 0; i < nsamples; i++ )
734 counters->data.i[labels->data.i[i]]++;
735
736 // ensure that we do not have empty clusters
737 for( k = 0; k < nclusters; k++ )
738 if( counters->data.i[k] == 0 )
739 for(;;)
740 {
741 i = cvRandInt(&rng) % nsamples;
742 j = labels->data.i[i];
743 if( counters->data.i[j] > 1 )
744 {
745 labels->data.i[i] = k;
746 counters->data.i[j]--;
747 counters->data.i[k]++;
748 break;
749 }
750 }
751
752 __END__;
753
754 cvReleaseMat( ¢ers );
755 cvReleaseMat( &old_centers );
756 cvReleaseMat( &counters );
757 }
758
759
760 /****************************************************************************************/
761 /* log_weight_div_det[k] = -2*log(weights_k) + log(det(Sigma_k)))
762
763 covs[k] = cov_rotate_mats[k] * cov_eigen_values[k] * (cov_rotate_mats[k])'
764 cov_rotate_mats[k] are orthogonal matrices of eigenvectors and
765 cov_eigen_values[k] are diagonal matrices (represented by 1D vectors) of eigen values.
766
767 The <alpha_ik> is the probability of the vector x_i to belong to the k-th cluster:
768 <alpha_ik> ~ weights_k * exp{ -0.5[ln(det(Sigma_k)) + (x_i - mu_k)' Sigma_k^(-1) (x_i - mu_k)] }
769 We calculate these probabilities here by the equivalent formulae:
770 Denote
771 S_ik = -0.5(log(det(Sigma_k)) + (x_i - mu_k)' Sigma_k^(-1) (x_i - mu_k)) + log(weights_k),
772 M_i = max_k S_ik = S_qi, so that the q-th class is the one where maximum reaches. Then
773 alpha_ik = exp{ S_ik - M_i } / ( 1 + sum_j!=q exp{ S_ji - M_i })
774 */
run_em(const CvVectors & train_data)775 double CvEM::run_em( const CvVectors& train_data )
776 {
777 CvMat* centered_sample = 0;
778 CvMat* covs_item = 0;
779 CvMat* log_det = 0;
780 CvMat* log_weights = 0;
781 CvMat* cov_eigen_values = 0;
782 CvMat* samples = 0;
783 CvMat* sum_probs = 0;
784 log_likelihood = -DBL_MAX;
785
786 CV_FUNCNAME( "CvEM::run_em" );
787 __BEGIN__;
788
789 int nsamples = train_data.count, dims = train_data.dims, nclusters = params.nclusters;
790 double min_variation = FLT_EPSILON;
791 double min_det_value = MAX( DBL_MIN, pow( min_variation, dims ));
792 double likelihood_bias = -CV_LOG2PI * (double)nsamples * (double)dims / 2., _log_likelihood = -DBL_MAX;
793 int start_step = params.start_step;
794
795 int i, j, k, n;
796 int is_general = 0, is_diagonal = 0, is_spherical = 0;
797 double prev_log_likelihood = -DBL_MAX / 1000., det, d;
798 CvMat whdr, iwhdr, diag, *w, *iw;
799 double* w_data;
800 double* sp_data;
801
802 if( nclusters == 1 )
803 {
804 double log_weight;
805 CV_CALL( cvSet( probs, cvScalar(1.)) );
806
807 if( params.cov_mat_type == COV_MAT_SPHERICAL )
808 {
809 d = cvTrace(*covs).val[0]/dims;
810 d = MAX( d, FLT_EPSILON );
811 inv_eigen_values->data.db[0] = 1./d;
812 log_weight = pow( d, dims*0.5 );
813 }
814 else
815 {
816 w_data = inv_eigen_values->data.db;
817
818 if( params.cov_mat_type == COV_MAT_GENERIC )
819 cvSVD( *covs, inv_eigen_values, *cov_rotate_mats, 0, CV_SVD_U_T );
820 else
821 cvTranspose( cvGetDiag(*covs, &diag), inv_eigen_values );
822
823 cvMaxS( inv_eigen_values, FLT_EPSILON, inv_eigen_values );
824 for( j = 0, det = 1.; j < dims; j++ )
825 det *= w_data[j];
826 log_weight = sqrt(det);
827 cvDiv( 0, inv_eigen_values, inv_eigen_values );
828 }
829
830 log_weight_div_det->data.db[0] = -2*log(weights->data.db[0]/log_weight);
831 log_likelihood = DBL_MAX/1000.;
832 EXIT;
833 }
834
835 if( params.cov_mat_type == COV_MAT_GENERIC )
836 is_general = 1;
837 else if( params.cov_mat_type == COV_MAT_DIAGONAL )
838 is_diagonal = 1;
839 else if( params.cov_mat_type == COV_MAT_SPHERICAL )
840 is_spherical = 1;
841 /* In the case of <cov_mat_type> == COV_MAT_DIAGONAL, the k-th row of cov_eigen_values
842 contains the diagonal elements (variations). In the case of
843 <cov_mat_type> == COV_MAT_SPHERICAL - the 0-ths elements of the vectors cov_eigen_values[k]
844 are to be equal to the mean of the variations over all the dimensions. */
845
846 CV_CALL( log_det = cvCreateMat( 1, nclusters, CV_64FC1 ));
847 CV_CALL( log_weights = cvCreateMat( 1, nclusters, CV_64FC1 ));
848 CV_CALL( covs_item = cvCreateMat( dims, dims, CV_64FC1 ));
849 CV_CALL( centered_sample = cvCreateMat( 1, dims, CV_64FC1 ));
850 CV_CALL( cov_eigen_values = cvCreateMat( inv_eigen_values->rows, inv_eigen_values->cols, CV_64FC1 ));
851 CV_CALL( samples = cvCreateMat( nsamples, dims, CV_64FC1 ));
852 CV_CALL( sum_probs = cvCreateMat( 1, nclusters, CV_64FC1 ));
853 sp_data = sum_probs->data.db;
854
855 // copy the training data into double-precision matrix
856 for( i = 0; i < nsamples; i++ )
857 {
858 const float* src = train_data.data.fl[i];
859 double* dst = (double*)(samples->data.ptr + samples->step*i);
860
861 for( j = 0; j < dims; j++ )
862 dst[j] = src[j];
863 }
864
865 if( start_step != START_M_STEP )
866 {
867 for( k = 0; k < nclusters; k++ )
868 {
869 if( is_general || is_diagonal )
870 {
871 w = cvGetRow( cov_eigen_values, &whdr, k );
872 if( is_general )
873 cvSVD( covs[k], w, cov_rotate_mats[k], 0, CV_SVD_U_T );
874 else
875 cvTranspose( cvGetDiag( covs[k], &diag ), w );
876 w_data = w->data.db;
877 for( j = 0, det = 1.; j < dims; j++ )
878 det *= w_data[j];
879 if( det < min_det_value )
880 {
881 if( start_step == START_AUTO_STEP )
882 det = min_det_value;
883 else
884 EXIT;
885 }
886 log_det->data.db[k] = det;
887 }
888 else
889 {
890 d = cvTrace(covs[k]).val[0]/(double)dims;
891 if( d < min_variation )
892 {
893 if( start_step == START_AUTO_STEP )
894 d = min_variation;
895 else
896 EXIT;
897 }
898 cov_eigen_values->data.db[k] = d;
899 log_det->data.db[k] = d;
900 }
901 }
902
903 cvLog( log_det, log_det );
904 if( is_spherical )
905 cvScale( log_det, log_det, dims );
906 }
907
908 for( n = 0; n < params.term_crit.max_iter; n++ )
909 {
910 if( n > 0 || start_step != START_M_STEP )
911 {
912 // e-step: compute probs_ik from means_k, covs_k and weights_k.
913 CV_CALL(cvLog( weights, log_weights ));
914
915 // S_ik = -0.5[log(det(Sigma_k)) + (x_i - mu_k)' Sigma_k^(-1) (x_i - mu_k)] + log(weights_k)
916 for( k = 0; k < nclusters; k++ )
917 {
918 CvMat* u = cov_rotate_mats[k];
919 const double* mean = (double*)(means->data.ptr + means->step*k);
920 w = cvGetRow( cov_eigen_values, &whdr, k );
921 iw = cvGetRow( inv_eigen_values, &iwhdr, k );
922 cvDiv( 0, w, iw );
923
924 w_data = (double*)(inv_eigen_values->data.ptr + inv_eigen_values->step*k);
925
926 for( i = 0; i < nsamples; i++ )
927 {
928 double *csample = centered_sample->data.db, p = log_det->data.db[k];
929 const double* sample = (double*)(samples->data.ptr + samples->step*i);
930 double* pp = (double*)(probs->data.ptr + probs->step*i);
931 for( j = 0; j < dims; j++ )
932 csample[j] = sample[j] - mean[j];
933 if( is_general )
934 cvGEMM( centered_sample, u, 1, 0, 0, centered_sample, CV_GEMM_B_T );
935 for( j = 0; j < dims; j++ )
936 p += csample[j]*csample[j]*w_data[is_spherical ? 0 : j];
937 pp[k] = -0.5*p + log_weights->data.db[k];
938
939 // S_ik <- S_ik - max_j S_ij
940 if( k == nclusters - 1 )
941 {
942 double max_val = 0;
943 for( j = 0; j < nclusters; j++ )
944 max_val = MAX( max_val, pp[j] );
945 for( j = 0; j < nclusters; j++ )
946 pp[j] -= max_val;
947 }
948 }
949 }
950
951 CV_CALL(cvExp( probs, probs )); // exp( S_ik )
952 cvZero( sum_probs );
953
954 // alpha_ik = exp( S_ik ) / sum_j exp( S_ij ),
955 // log_likelihood = sum_i log (sum_j exp(S_ij))
956 for( i = 0, _log_likelihood = likelihood_bias; i < nsamples; i++ )
957 {
958 double* pp = (double*)(probs->data.ptr + probs->step*i), sum = 0;
959 for( j = 0; j < nclusters; j++ )
960 sum += pp[j];
961 sum = 1./MAX( sum, DBL_EPSILON );
962 for( j = 0; j < nclusters; j++ )
963 {
964 double p = pp[j] *= sum;
965 sp_data[j] += p;
966 }
967 _log_likelihood -= log( sum );
968 }
969
970 // check termination criteria
971 if( fabs( (_log_likelihood - prev_log_likelihood) / prev_log_likelihood ) < params.term_crit.epsilon )
972 break;
973 prev_log_likelihood = _log_likelihood;
974 }
975
976 // m-step: update means_k, covs_k and weights_k from probs_ik
977 cvGEMM( probs, samples, 1, 0, 0, means, CV_GEMM_A_T );
978
979 for( k = 0; k < nclusters; k++ )
980 {
981 double sum = sp_data[k], inv_sum = 1./sum;
982 CvMat* cov = covs[k], _mean, _sample;
983
984 w = cvGetRow( cov_eigen_values, &whdr, k );
985 w_data = w->data.db;
986 cvGetRow( means, &_mean, k );
987 cvGetRow( samples, &_sample, k );
988
989 // update weights_k
990 weights->data.db[k] = sum;
991
992 // update means_k
993 cvScale( &_mean, &_mean, inv_sum );
994
995 // compute covs_k
996 cvZero( cov );
997 cvZero( w );
998
999 for( i = 0; i < nsamples; i++ )
1000 {
1001 double p = probs->data.db[i*nclusters + k]*inv_sum;
1002 _sample.data.db = (double*)(samples->data.ptr + samples->step*i);
1003
1004 if( is_general )
1005 {
1006 cvMulTransposed( &_sample, covs_item, 1, &_mean );
1007 cvScaleAdd( covs_item, cvRealScalar(p), cov, cov );
1008 }
1009 else
1010 for( j = 0; j < dims; j++ )
1011 {
1012 double val = _sample.data.db[j] - _mean.data.db[j];
1013 w_data[is_spherical ? 0 : j] += p*val*val;
1014 }
1015 }
1016
1017 if( is_spherical )
1018 {
1019 d = w_data[0]/(double)dims;
1020 d = MAX( d, min_variation );
1021 w->data.db[0] = d;
1022 log_det->data.db[k] = d;
1023 }
1024 else
1025 {
1026 if( is_general )
1027 cvSVD( cov, w, cov_rotate_mats[k], 0, CV_SVD_U_T );
1028 cvMaxS( w, min_variation, w );
1029 for( j = 0, det = 1.; j < dims; j++ )
1030 det *= w_data[j];
1031 log_det->data.db[k] = det;
1032 }
1033 }
1034
1035 cvConvertScale( weights, weights, 1./(double)nsamples, 0 );
1036 cvMaxS( weights, DBL_MIN, weights );
1037
1038 cvLog( log_det, log_det );
1039 if( is_spherical )
1040 cvScale( log_det, log_det, dims );
1041 } // end of iteration process
1042
1043 //log_weight_div_det[k] = -2*log(weights_k/det(Sigma_k))^0.5) = -2*log(weights_k) + log(det(Sigma_k)))
1044 if( log_weight_div_det )
1045 {
1046 cvScale( log_weights, log_weight_div_det, -2 );
1047 cvAdd( log_weight_div_det, log_det, log_weight_div_det );
1048 }
1049
1050 /* Now finalize all the covariation matrices:
1051 1) if <cov_mat_type> == COV_MAT_DIAGONAL we used array of <w> as diagonals.
1052 Now w[k] should be copied back to the diagonals of covs[k];
1053 2) if <cov_mat_type> == COV_MAT_SPHERICAL we used the 0-th element of w[k]
1054 as an average variation in each cluster. The value of the 0-th element of w[k]
1055 should be copied to the all of the diagonal elements of covs[k]. */
1056 if( is_spherical )
1057 {
1058 for( k = 0; k < nclusters; k++ )
1059 cvSetIdentity( covs[k], cvScalar(cov_eigen_values->data.db[k]));
1060 }
1061 else if( is_diagonal )
1062 {
1063 for( k = 0; k < nclusters; k++ )
1064 cvTranspose( cvGetRow( cov_eigen_values, &whdr, k ),
1065 cvGetDiag( covs[k], &diag ));
1066 }
1067 cvDiv( 0, cov_eigen_values, inv_eigen_values );
1068
1069 log_likelihood = _log_likelihood;
1070
1071 __END__;
1072
1073 cvReleaseMat( &log_det );
1074 cvReleaseMat( &log_weights );
1075 cvReleaseMat( &covs_item );
1076 cvReleaseMat( ¢ered_sample );
1077 cvReleaseMat( &cov_eigen_values );
1078 cvReleaseMat( &samples );
1079 cvReleaseMat( &sum_probs );
1080
1081 return log_likelihood;
1082 }
1083
1084
get_nclusters() const1085 int CvEM::get_nclusters() const
1086 {
1087 return params.nclusters;
1088 }
1089
get_means() const1090 const CvMat* CvEM::get_means() const
1091 {
1092 return means;
1093 }
1094
get_covs() const1095 const CvMat** CvEM::get_covs() const
1096 {
1097 return (const CvMat**)covs;
1098 }
1099
get_weights() const1100 const CvMat* CvEM::get_weights() const
1101 {
1102 return weights;
1103 }
1104
get_probs() const1105 const CvMat* CvEM::get_probs() const
1106 {
1107 return probs;
1108 }
1109
1110 /* End of file. */
1111