运动检测（前景检测）之（二）混合高斯模型GMM

因为监控发展的需求，目前前景检测的研究还是很多的，也出现了很多新的方法和思路。个人了解的大概概括为以下一些：

帧差、背景减除（GMM、CodeBook、 SOBS、 SACON、 VIBE、 W4、多帧平均……）、光流（稀疏光流、稠密光流）、运动竞争（Motion Competition）、运动模版（运动历史图像）、时间熵……等等。如果加上他们的改进版，那就是很大的一个家族了。

        对于上一些方法的一点简单的对比分析可以参考下：

http://www.cnblogs.com/ronny/archive/2012/04/12/2444053.html

        至于哪个最好，看使用环境吧，各有千秋，有一些适用的情况更多，有一些在某些情况下表现更好。这些都需要针对自己的使用情况作测试确定的。呵呵。

        推荐一个牛逼的库：http://code.google.com/p/bgslibrary/里面包含了各种背景减除的方法，可以让自己少做很多力气活。

        还有王先荣博客上存在不少的分析：

http://www.cnblogs.com/xrwang/archive/2010/02/21/ForegroundDetection.html

        下面的博客上转载王先荣的上面几篇，然后加上自己分析了两篇：

http://blog.csdn.net/stellar0

        本文主要关注其中的一种背景减除方法：GMM。tornadomeet的博客上对ViBe进行了分析，我这里就不再啰嗦了，具体的理论分析可以参考：

http://www.cnblogs.com/tornadomeet/archive/2012/06/02/2531565.html

        里面有了GMM的代码，并有了详细的注释。我之前根据这个代码（在这里，非常感谢tornadomeet）改写了一个Mat格式的版本，现在发上来和大家交流，具体如下：（在VS2010+OpenCV2.4.2中测试通过）。（当然了，OpenCV也已经提供了MOG的背景减除方法）

MOG_BGS.h

[cpp] view plaincopy
  1. #pragma once  
  2. #include <iostream>  
  3. #include "opencv2/opencv.hpp"  
  4.   
  5. using namespace cv;  
  6. using namespace std;  
  7.   
  8. //定义gmm模型用到的变量  
  9.  #define GMM_MAX_COMPONT 6          //每个GMM最多的高斯模型个数  
  10.  #define GMM_LEARN_ALPHA 0.005    
  11.  #define GMM_THRESHOD_SUMW 0.7  
  12.  #define TRAIN_FRAMES 60    // 对前 TRAIN_FRAMES 帧建模  
  13.   
  14. class MOG_BGS  
  15. {  
  16. public:  
  17.     MOG_BGS(void);  
  18.     ~MOG_BGS(void);  
  19.   
  20.     void init(const Mat _image);  
  21.     void processFirstFrame(const Mat _image);  
  22.     void trainGMM(const Mat _image);  
  23.     void getFitNum(const Mat _image);  
  24.     void testGMM(const Mat _image);  
  25.     Mat getMask(void){return m_mask;};  
  26.    
  27. private:  
  28.     Mat m_weight[GMM_MAX_COMPONT];  //权值  
  29.     Mat m_mean[GMM_MAX_COMPONT];    //均值  
  30.     Mat m_sigma[GMM_MAX_COMPONT];   //方差  
  31.   
  32.     Mat m_mask;  
  33.     Mat m_fit_num;  
  34. };  

MOG_BGS.cpp

[cpp] view plaincopy
  1. #include "MOG_BGS.h"  
  2.   
  3. MOG_BGS::MOG_BGS(void)  
  4. {  
  5.   
  6. }  
  7.   
  8. MOG_BGS::~MOG_BGS(void)  
  9. {  
  10.   
  11. }  
  12.   
  13. // 全部初始化为0  
  14. void MOG_BGS::init(const Mat _image)  
  15. {  
  16.     /****initialization the three parameters ****/  
  17.      for(int i = 0; i < GMM_MAX_COMPONT; i++)  
  18.      {  
  19.          m_weight[i] = Mat::zeros(_image.size(), CV_32FC1);  
  20.          m_mean[i] = Mat::zeros(_image.size(), CV_8UC1);  
  21.          m_sigma[i] = Mat::zeros(_image.size(), CV_32FC1);  
  22.      }  
  23.      m_mask = Mat::zeros(_image.size(),CV_8UC1);  
  24.      m_fit_num = Mat::ones(_image.size(),CV_8UC1);  
  25. }  
  26.   
  27. //gmm第一帧初始化函数实现  
  28. //捕获到第一帧时对高斯分布进行初始化．主要包括对每个高斯分布的权值、期望和方差赋初值．  
  29. //其中第一个高斯分布的权值为1，期望为第一个像素数据．其余高斯分布权值为0，期望为0．  
  30. //每个高斯分布都被赋予适当的相等的初始方差 15  
  31. void MOG_BGS::processFirstFrame(const Mat _image)  
  32. {  
  33.     for(int i = 0; i < GMM_MAX_COMPONT; i++)  
  34.     {  
  35.         if (i == 0)  
  36.         {  
  37.             m_weight[i].setTo(1.0);  
  38.             _image.copyTo(m_mean[i]);  
  39.             m_sigma[i].setTo(15.0);  
  40.         }  
  41.         else  
  42.         {  
  43.             m_weight[i].setTo(0.0);  
  44.             m_mean[i].setTo(0);  
  45.             m_sigma[i].setTo(15.0);  
  46.         }  
  47.     }  
  48. }  
  49.    
  50. // 通过新的帧来训练GMM  
  51. void MOG_BGS::trainGMM(const Mat _image)  
  52. {  
  53.     for(int i = 0; i < _image.rows; i++)  
  54.     {  
  55.         for(int j = 0; j < _image.cols; j++)  
  56.         {  
  57.              int num_fit = 0;  
  58.   
  59.              /**************************** Update parameters Start ******************************************/  
  60.              for(int k = 0 ; k < GMM_MAX_COMPONT; k++)  
  61.              {  
  62.                  int delm = abs(_image.at<uchar>(i, j) - m_mean[k].at<uchar>(i, j));  
  63.                  long dist = delm * delm;  
  64.                  // 判断是否匹配：采样值与高斯分布的均值的距离小于3倍方差（表示匹配）  
  65.                  if( dist < 3.0 * m_sigma[k].at<float>(i, j))   
  66.                  {  
  67.                      // 如果匹配  
  68.                      /****update the weight****/  
  69.                      m_weight[k].at<float>(i, j) += GMM_LEARN_ALPHA * (1 - m_weight[k].at<float>(i, j));  
  70.    
  71.                      /****update the average****/  
  72.                      m_mean[k].at<uchar>(i, j) += (GMM_LEARN_ALPHA / m_weight[k].at<uchar>(i, j)) * delm;  
  73.    
  74.                      /****update the variance****/  
  75.                      m_sigma[k].at<float>(i, j) += (GMM_LEARN_ALPHA / m_weight[k].at<float>(i, j)) * (dist - m_sigma[k].at<float>(i, j));  
  76.                  }  
  77.                  else  
  78.                  {  
  79.                     // 如果不匹配。则该该高斯模型的权值变小  
  80.                      m_weight[k].at<float>(i, j) += GMM_LEARN_ALPHA * (0 - m_weight[k].at<float>(i, j));  
  81.                      num_fit++; // 不匹配的模型个数  
  82.                  }          
  83.              }  
  84.              /**************************** Update parameters End ******************************************/      
  85.           
  86.   
  87.              /*********************** Sort Gaussian component by 'weight / sigma' Start ****************************/  
  88.              //对gmm各个高斯进行排序,从大到小排序,排序依据为 weight / sigma  
  89.              for(int kk = 0; kk < GMM_MAX_COMPONT; kk++)  
  90.              {  
  91.                  for(int rr=kk; rr< GMM_MAX_COMPONT; rr++)  
  92.                  {  
  93.                      if(m_weight[rr].at<float>(i, j)/m_sigma[rr].at<float>(i, j) > m_weight[kk].at<float>(i, j)/m_sigma[kk].at<float>(i, j))  
  94.                      {  
  95.                          //权值交换  
  96.                          float temp_weight = m_weight[rr].at<float>(i, j);  
  97.                          m_weight[rr].at<float>(i, j) = m_weight[kk].at<float>(i, j);  
  98.                          m_weight[kk].at<float>(i, j) = temp_weight;  
  99.    
  100.                          //均值交换  
  101.                          uchar temp_mean = m_mean[rr].at<uchar>(i, j);  
  102.                          m_mean[rr].at<uchar>(i, j) = m_mean[kk].at<uchar>(i, j);  
  103.                          m_mean[kk].at<uchar>(i, j) = temp_mean;  
  104.    
  105.                          //方差交换  
  106.                          float temp_sigma = m_sigma[rr].at<float>(i, j);  
  107.                          m_sigma[rr].at<float>(i, j) = m_sigma[kk].at<float>(i, j);  
  108.                          m_sigma[kk].at<float>(i, j) = temp_sigma;  
  109.                      }  
  110.                  }  
  111.              }  
  112.              /*********************** Sort Gaussian model by 'weight / sigma' End ****************************/  
  113.    
  114.   
  115.              /*********************** Create new Gaussian component Start ****************************/  
  116.              if(num_fit == GMM_MAX_COMPONT && 0 == m_weight[GMM_MAX_COMPONT - 1].at<float>(i, j))  
  117.              {  
  118.                  //if there is no exit component fit,then start a new component  
  119.                  //当有新值出现的时候，若目前分布个数小于M，新添一个分布，以新采样值作为均值，并赋予较大方差和较小权值  
  120.                   for(int k = 0 ; k < GMM_MAX_COMPONT; k++)  
  121.                  {  
  122.                      if(0 == m_weight[k].at<float>(i, j))  
  123.                      {  
  124.                          m_weight[k].at<float>(i, j) = GMM_LEARN_ALPHA;  
  125.                          m_mean[k].at<uchar>(i, j) = _image.at<uchar>(i, j);  
  126.                          m_sigma[k].at<float>(i, j) = 15.0;  
  127.                           
  128.                          //normalization the weight,let they sum to 1  
  129.                          for(int q = 0; q < GMM_MAX_COMPONT && q != k; q++)  
  130.                          {  
  131.                             //对其他的高斯模型的权值进行更新，保持权值和为1  
  132.                              /****update the other unfit's weight,u and sigma remain unchanged****/  
  133.                              m_weight[q].at<float>(i, j) *= (1 - GMM_LEARN_ALPHA);  
  134.                          }  
  135.                          break; //找到第一个权值不为0的即可  
  136.                       }                              
  137.                   }  
  138.              }  
  139.              else if(num_fit == GMM_MAX_COMPONT && m_weight[GMM_MAX_COMPONT -1].at<float>(i, j) != 0)  
  140.              {  
  141.                  //如果GMM_MAX_COMPONT都曾经赋值过，则用新来的高斯代替权值最弱的高斯，权值不变，只更新均值和方差  
  142.                  m_mean[GMM_MAX_COMPONT-1].at<uchar>(i, j) = _image.at<uchar>(i, j);  
  143.                  m_sigma[GMM_MAX_COMPONT-1].at<float>(i, j) = 15.0;  
  144.              }  
  145.              /*********************** Create new Gaussian component End ****************************/  
  146.          }  
  147.     }  
  148. }  
  149.   
  150.  //对输入图像每个像素gmm选择合适的高斯分量个数  
  151.  //排序后最有可能是背景分布的排在最前面，较小可能的短暂的分布趋向于末端．我们将排序后的前fit_num个分布选为背景模型;  
  152.  //在排过序的分布中，累积概率超过GMM_THRESHOD_SUMW的前fit_num个分布被当作背景模型，剩余的其它分布被当作前景模型．  
  153. void MOG_BGS::getFitNum(const Mat _image)  
  154. {  
  155.     for(int i = 0; i < _image.rows; i++)  
  156.     {  
  157.         for(int j = 0; j < _image.cols; j++)  
  158.         {  
  159.             float sum_w = 0.0;  //重新赋值为0，给下一个像素做累积  
  160.             for(uchar k = 0; k < GMM_MAX_COMPONT; k++)  
  161.             {  
  162.                 sum_w += m_weight[k].at<float>(i, j);  
  163.                 if(sum_w >= GMM_THRESHOD_SUMW)   //如果这里THRESHOD_SUMW=0.6的话，那么得到的高斯数目都为1，因为每个像素都有一个权值接近1  
  164.                 {  
  165.                      m_fit_num.at<uchar>(i, j) = k + 1;  
  166.                      break;  
  167.                 }  
  168.             }  
  169.         }  
  170.     }  
  171. }  
  172.   
  173.  //gmm测试函数的实现  
  174. void MOG_BGS::testGMM(const Mat _image)  
  175. {  
  176.     for(int i = 0; i < _image.rows; i++)  
  177.     {  
  178.         for(int j = 0; j < _image.cols; j++)  
  179.         {  
  180.             int k = 0;  
  181.             for( ; k < m_fit_num.at<uchar>(i, j); k++)  
  182.             {  
  183.                 if(abs(_image.at<uchar>(i, j) - m_mean[k].at<uchar>(i, j)) < (uchar)( 2.5 * m_sigma[k].at<float>(i, j)))  
  184.                 {  
  185.                     m_mask.at<uchar>(i, j) = 0;  
  186.                     break;  
  187.                 }  
  188.             }  
  189.             if(k == m_fit_num.at<uchar>(i, j))  
  190.             {  
  191.                 m_mask.at<uchar>(i, j) = 255;  
  192.             }  
  193.         }  
  194.     }  
  195. }  

Main.cpp

[cpp] view plaincopy
  1. // This is based on the "An Improved Adaptive Background Mixture Model for  
  2. // Real-time Tracking with Shadow Detection" by P. KaewTraKulPong and R. Bowden  
  3. // Author : zouxy  
  4. // Date   : 2013-4-13  
  5. // HomePage : http://blog.csdn.net/zouxy09  
  6. // Email  : zouxy09@qq.com  
  7.   
  8. #include "opencv2/opencv.hpp"  
  9. #include "MOG_BGS.h"  
  10. #include <iostream>  
  11. #include <cstdio>  
  12.   
  13. using namespace cv;  
  14. using namespace std;  
  15.   
  16. int main(int argc, char* argv[])  
  17. {  
  18.     Mat frame, gray, mask;  
  19.     VideoCapture capture;  
  20.     capture.open("video.avi");  
  21.   
  22.     if (!capture.isOpened())  
  23.     {  
  24.         cout<<"No camera or video input! "<<endl;  
  25.         return -1;  
  26.     }  
  27.   
  28.     MOG_BGS Mog_Bgs;  
  29.     int count = 0;  
  30.   
  31.     while (1)  
  32.     {  
  33.         count++;  
  34.         capture >> frame;  
  35.         if (frame.empty())  
  36.             break;  
  37.         cvtColor(frame, gray, CV_RGB2GRAY);  
  38.       
  39.         if (count == 1)  
  40.         {  
  41.             Mog_Bgs.init(gray);  
  42.             Mog_Bgs.processFirstFrame(gray);  
  43.             cout<<" Using "<<TRAIN_FRAMES<<" frames to training GMM..."<<endl;  
  44.         }  
  45.         else if (count < TRAIN_FRAMES)  
  46.         {  
  47.             Mog_Bgs.trainGMM(gray);  
  48.         }  
  49.         else if (count == TRAIN_FRAMES)  
  50.         {  
  51.             Mog_Bgs.getFitNum(gray);  
  52.             cout<<" Training GMM complete!"<<endl;  
  53.         }  
  54.         else  
  55.         {  
  56.             Mog_Bgs.testGMM(gray);  
  57.             mask = Mog_Bgs.getMask();  
  58.             morphologyEx(mask, mask, MORPH_OPEN, Mat());  
  59.             erode(mask, mask, Mat(7, 7, CV_8UC1), Point(-1, -1));  // You can use Mat(5, 5, CV_8UC1) here for less distortion  
  60.             dilate(mask, mask, Mat(7, 7, CV_8UC1), Point(-1, -1));  
  61.             imshow("mask", mask);  
  62.         }  
  63.   
  64.         imshow("input", frame);   
  65.   
  66.         if ( cvWaitKey(10) == 'q' )  
  67.             break;  
  68.     }  
  69.   
  70.     return 0;  
  71. }  
【推广】 免费学中医，健康全家人
原文地址：https://www.cnblogs.com/ywsoftware/p/4434309.html