Expectation-Maximization in CSharp

This example requires Emgu CV 1.5.0.0

using System.Drawing;
using Emgu.CV.Structure;
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);
}