# Gaussian Mixture Models + EM algorithm - Demo

OpenCV has an EM reference implementation:

## Demo 1a (general covariance matrices for Gaussians)

Here I generated some samples on the letters of the word “hello” and tried to learn a Gaussian Mixture Model (GMM) with n=1,2,…,100 components.

I.e. describing the sample distribution with 1,2,…,100 Gaussians where each Gaussian is described by

• a mean 2D position
• a (full) 2×2 covariance matrix
• and a weight (remember: GMM := sum of weighted Gaussians)

Conclusions:

• it roughly works
• the more gaussian you use in your GMM model, the longer the EM training takes
• but it is still “fast”, i.e. in the order of less or just a few seconds for all n (1,2,…,100)
• note! if the weight of a gaussian is 0.0 after learning:
• you get an invalid cluster center with location cx=-1.#INF, cy=-1.#INF for that Gaussian
• bug in OpenCV's implementation: sometimes the resulting GMM has one component that covers all the sample area and has weight 1.0 while all other GMM components have weight 0.0 (see animation)
• practical solution for that bug: re-start GMM training with different start conditions? (random seed)
• note (from 2012/05/09): this bug has recently be fixed in the trunk version of OpenCV!

BTW: the weight of each Gaussians is visualized by the “white-ness” of the corresponding ellipse's line or fill-color, respectively.

## Demo 1b (spherical covariance matrices)

As demo 1a, but now we constrain the covariance matrix of the Gaussians to be spherical, i.e.

```    | var   0 |
M = | 0   var |```

## Demo 1c (diagonal covariance matrices)

As demo 1a, but now we constrain the covariance matrix of the Gaussians to be diagonal, i.e.

```    | var_x     0 |
M = | 0      var_y|```

## Demo 2

Similar to demo 1 but with another sample distribution: see the red crosses.

## Sample source code

```bool IsNumber(double x)
{
// This looks like it should always be true,
// but it's false if x is a NaN.
return (x == x);
}

GMM::GMM()
{
// 1. Prepare training vectors

// 1.1 Create a list of sample points from a RGB image (some image structure, e.g. text on white background)
int SampleImgHeight = SampleImg.rows;
int SampleImgWidth  = SampleImg.cols;
QList<QPoint> ListSamplePoints;
for (int y=0; y<SampleImgHeight; y++)
{
for (int x=0; x<SampleImgWidth; x++)
{
// Get pixel color at that position
Vec3b bgrPixel = SampleImg.at<Vec3b>(y, x);

uchar b = bgrPixel.val[0];
//uchar g = bgrPixel.val[1];
//uchar r = bgrPixel.val[2];

if (b==0)
{
if (rand() % 25 == 0)
{
ListSamplePoints.append( QPoint(x,y) );

} // (add to list of sample points?)

} // if (color is not white)

} // for (x)
} // for (y)

// 1.2 Write the just generated sample list points in <ListSamplePoints> to a matrix
Mat labels;
int NrSamples = ListSamplePoints.count();
Mat samples( NrSamples, 2, CV_32FC1 );
for (int s=0; s<NrSamples; s++)
{
QPoint p = ListSamplePoints.at(s);

samples.at<float>(s,0) = (float) p.x();
samples.at<float>(s,1) = (float) p.y();
}
Logger::Access()->Log( QString("Sample matrix dimensions: %0 x %1").arg( samples.rows ).arg( samples.cols ) );

// Try out different nrs of GMM components to represent the sample distribution
for (int NrGMMComponents = 1; NrGMMComponents<=100; NrGMMComponents++)
{
Logger::Access()->Log( QString("\nLearning to represent the sample distributions with %0 gaussians.").arg( NrGMMComponents ) );

// 2. Prepare GMM training

// 2.1 Set GMM parameters
CvEMParams params;
params.covs      = NULL;
params.means     = NULL;
params.weights   = NULL;
params.probs     = NULL;
params.nclusters = NrGMMComponents;
params.cov_mat_type       = CvEM::COV_MAT_GENERIC; // DIAGONAL, GENERIC, SPHERICAL
params.start_step         = CvEM::START_AUTO_STEP;
params.term_crit.max_iter = 1500;
params.term_crit.epsilon  = 0.001;
params.term_crit.type     = CV_TERMCRIT_ITER|CV_TERMCRIT_EPS;
//params.term_crit.type     = CV_TERMCRIT_ITER;

// 2.2 Estimate GMM params for all <NrGMMComponents> Gaussian Mixture Components
Logger::Access()->Log( "Started GMM training" );
CvEM em_model;
em_model.train( samples, Mat(), params, &labels );
Logger::Access()->Log( "Finished GMM training" );

// 3. Visualize results

// 3.1 Prepare visualization images
Mat img  = Mat::zeros( Size( SampleImgWidth, SampleImgHeight ), CV_8UC3 );
Mat img2 = Mat::zeros( Size( SampleImgWidth, SampleImgHeight ), CV_8UC3 );

char txt[500];
sprintf_s(txt, "nr of gaussians: %d", NrGMMComponents);
putText(img,  txt, Point(10,30), CV_FONT_HERSHEY_SIMPLEX, .75, CV_RGB(0,255,0) );
putText(img2, txt, Point(10,30), CV_FONT_HERSHEY_SIMPLEX, .75, CV_RGB(0,255,0) );

// 3.2 Show classification for each pixel to which GMM component it belongs?
if (false)
{
// Classify every image pixel
// A single sample is a 1x2 float matrix
Mat sample( 1, 2, CV_32FC1 );
for(int i = 0; i < img.rows; i++ )
{
for(int j = 0; j < img.cols; j++ )
{
sample.at<float>(0) = (float)j;
sample.at<float>(1) = (float)i;
int response = cvRound(em_model.predict( sample ));
//Scalar c = colors[response];
//circle( img, Point(j, i), 1, c*0.75, CV_FILLED );
}
}
}

// 3.3 Draw the samples
for (int i = 0; i < NrSamples; i++)
{
float x = samples.at<float>(i,0);
float y = samples.at<float>(i,1);
Point pt = Point(x,y);
circle( img, pt, 1, CV_RGB(255,0,0), CV_FILLED );
}

// 3.4 Draw the cluster centers
cv::Mat Means   = em_model.getMeans();
cv::Mat Weights = em_model.getWeights();
std::vector<cv::Mat> covs;
em_model.getCovs(covs);
for (int i=0; i<params.nclusters; i++)
{
double cx = Means.at<double>(i,0);
double cy = Means.at<double>(i,1);
double wi = Weights.at<double>(0,i);
Logger::Access()->Log( QString("Weight for cluster #%0 = %1").arg(i).arg(wi,0,'f',4) );

// Ups! Invalid cluster center!
if (!IsNumber(cx))
{
Logger::Access()->Log( "Invalid cluster center! This should not happen! But it happens! :( => So it seems there is an error in OpenCV's EM algorithm... " );
continue;
}

circle( img, cv::Point(cx,cy), 5, CV_RGB(255,255,255), CV_FILLED );

cv::Mat CurrCovMat = covs.at(i);
double m11 = CurrCovMat.at<double>(0,0);
double m12 = CurrCovMat.at<double>(0,1);

double m21 = CurrCovMat.at<double>(1,0);
double m22 = CurrCovMat.at<double>(1,1);

// Get Eigenvectors + Eigenvalues of that covariance matrix
Mat eigenvalues;
Mat eigenvectors;
eigen( CurrCovMat, eigenvalues, eigenvectors );

if (false)
{
Logger::Access()->Log( QString("eigenvectors is a %0 x %1 matrix.").arg( eigenvectors.rows ).arg( eigenvectors.cols ) );
Logger::Access()->Log( QString("eigenvalues  is a %0 x %1 matrix.").arg( eigenvalues.rows ).arg( eigenvalues.cols ) );
for (int r=0; r<=1; r++)
for (int c=0; c<=1; c++)
Logger::Access()->Log( QString("   m%0%1 = %2").arg(r).arg(c).arg( eigenvectors.at<double>(r,c) ) );
for (int i=0; i<=1; i++)
Logger::Access()->Log( QString("Length of eigenvec #%0 = %1").arg(i).arg( eigenvalues.at<double>(i,0) ) );
}

double eigenvec1_len = eigenvalues.at<double>(0,0);
double len1          = sqrt(eigenvec1_len)*3;
double eigenvec1_x   = eigenvectors.at<double>(0,0) * len1;
double eigenvec1_y   = eigenvectors.at<double>(1,0) * len1;

double eigenvec2_len = eigenvalues.at<double>(1,0);
double len2          = sqrt(eigenvec2_len)*3;
double eigenvec2_x   = eigenvectors.at<double>(0,1) * len2;
double eigenvec2_y   = eigenvectors.at<double>(1,1) * len2;

/*
if ( (NrGMMComponents!=1) && ((len1>SampleImgHeight/2) || (len1>SampleImgWidth/2)) )
{
Logger::Access()->Log( "Yet another OpenCV EM algorithm error..." );
continue;
}
*/

// Show axes of ellipse
line( img, cv::Point(cx,cy), cv::Point(cx+eigenvec1_x, cy+eigenvec1_y), CV_RGB(255,255,0) );
line( img, cv::Point(cx,cy), cv::Point(cx+eigenvec2_x, cy+eigenvec2_y), CV_RGB(0,255,0) );

// Show ellipse rotated into direction of eigenvec1
double dx = eigenvec1_x;
double dy = eigenvec1_y;
double angle_rad = atan2( dy, dx );
double angle_deg = angle_rad * (180.0/M_PI); // convert radians (0,2PI) to degree (0°,360°)
cv::RotatedRect* myRotatedRect = new cv::RotatedRect( cvPoint(cx,cy), cvSize(len1, len2), angle_deg );
if (myRotatedRect != nullptr)
{
int g = wi*255.0 * (i+1);
cv::Scalar s = CV_RGB(g,g,g);
ellipse( img, *myRotatedRect, s );
ellipse( img2, *myRotatedRect, s, -1 );

delete myRotatedRect;
}
}

// 3.5 Save image
char filename1[500];
char filename2[500];
sprintf_s(filename1, "W:\\tmp\\visu1_sample_distribution_represented_with_%03d_gaussians.png", NrGMMComponents);
sprintf_s(filename2, "W:\\tmp\\visu2_sample_distribution_represented_with_%03d_gaussians.png", NrGMMComponents);
cv::imwrite(filename1, img);
cv::imwrite(filename2, img2);

} // for (NrGMMComponents)

}```