Expectation-Maximization in CSharp

This example requires Emgu CV 1.5.0.0

What is an Expectaion-Maximization Classifier

According to wikipedia,

An expectation-maximization (EM) algorithm is used in statistics for finding maximum likelihood estimates of parameters in probabilistic models, where the model depends on unobserved latent variables. EM is an iterative method which alternates between performing an expectation (E) step, which computes an expectation of the log likelihood with respect to the current estimate of the distribution for the latent variables, and a maximization (M) step, which computes the parameters which maximize the expected log likelihood found on the E step. These parameters are then used to determine the distribution of the latent variables in the next E step.

Source Code

using System.Drawing;
using Emgu.CV.Structure;
using Emgu.CV.ML;
using Emgu.CV.ML.Structure;

...

int N = 4; //number of clusters
int N1 = (int)Math.Sqrt((double)4);

Bgr[] colors = new Bgr[] { 
   new Bgr(0, 0, 255), 
   new Bgr(0, 255, 0),
   new Bgr(0, 255, 255),
   new Bgr(255, 255, 0)};

int nSamples = 100;

Matrix<float> samples = new Matrix<float>(nSamples, 2);
Matrix<Int32> labels = new Matrix<int>(nSamples, 1);
Image<Bgr, Byte> img = new Image<Bgr,byte>(500, 500);
Matrix<float> sample = new Matrix<float>(1, 2);

CvInvoke.cvReshape(samples.Ptr, samples.Ptr, 2, 0);
for (int i = 0; i < N; i++)
{
   Matrix<float> rows = samples.GetRows(i * nSamples / N, (i + 1) * nSamples / N, 1);
   double scale = ((i % N1) + 1.0) / (N1 + 1);
   MCvScalar mean = new MCvScalar(scale * img.Width, scale * img.Height);
   MCvScalar sigma = new MCvScalar(30, 30);
   ulong seed = (ulong)DateTime.Now.Ticks;
   CvInvoke.cvRandArr(ref seed, rows.Ptr, Emgu.CV.CvEnum.RAND_TYPE.CV_RAND_NORMAL, mean, sigma);
}
CvInvoke.cvReshape(samples.Ptr, samples.Ptr, 1, 0);

using (EM emModel1 = new EM())
using (EM emModel2 = new EM())
{
   EMParams parameters1 = new EMParams();
   parameters1.Nclusters = N;
   parameters1.CovMatType = Emgu.CV.ML.MlEnum.EM_COVARIAN_MATRIX_TYPE.COV_MAT_DIAGONAL;
   parameters1.StartStep = Emgu.CV.ML.MlEnum.EM_INIT_STEP_TYPE.START_AUTO_STEP;
   parameters1.TermCrit = new MCvTermCriteria(10, 0.01);
   emModel1.Train(samples, null, parameters1, labels);

   EMParams parameters2 = new EMParams();
   parameters2.Nclusters = N;
   parameters2.CovMatType = Emgu.CV.ML.MlEnum.EM_COVARIAN_MATRIX_TYPE.COV_MAT_GENERIC;
   parameters2.StartStep = Emgu.CV.ML.MlEnum.EM_INIT_STEP_TYPE.START_E_STEP;
   parameters2.TermCrit = new MCvTermCriteria(100, 1.0e-6);
   parameters2.Means = emModel1.GetMeans();
   parameters2.Covs = emModel1.GetCovariances();
   parameters2.Weights = emModel1.GetWeights();
               
   emModel2.Train(samples, null, parameters2, labels);
   
   #region Classify every image pixel
   for (int i = 0; i < img.Height; i++)
      for (int j = 0; j < img.Width; j++)
      {
         sample.Data[0, 0] = i;
         sample.Data[0, 1] = j;
         int response = (int) emModel2.Predict(sample, null);

         Bgr color = colors[response];
         
         img.Draw(
            new CircleF(new PointF(i, j), 1), 
            new Bgr(color.Blue*0.5, color.Green * 0.5, color.Red * 0.5 ), 
            0);
      }
   #endregion 

   #region draw the clustered samples
   for (int i = 0; i < nSamples; i++)
   {
      img.Draw(new CircleF(new PointF(samples.Data[i, 0], samples.Data[i, 1]), 1), colors[labels.Data[i, 0]], 0);
   }
   #endregion 

   Emgu.CV.UI.ImageViewer.Show(img);
}

Result