Gaussian Mixture Models + EM algorithm - Demo

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)
    //cv::Mat SampleImg = imread("W:\\01_data\\07_testimgs\\blobs_and_rectangles.png");
    cv::Mat SampleImg = imread("W:\\01_data\\07_testimgs\\hello.png");
    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)
 
}
 
public/gaussian_mixture_models_em_algorithm_-_demo.txt · Last modified: 2012/05/09 11:21 (external edit) · []
Recent changes RSS feed Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki