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

时间：2023-03-09 14:41:11
[转]运动检测（前景检测）之（二）混合高斯模型GMM

转自：http://blog.****.net/zouxy09/article/details/9622401

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

帧差、背景减除（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.****.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

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

MOG_BGS.cpp

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

Main.cpp

  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.****.net/zouxy09
  6. // Email  : zouxy09@qq.com
  7. #include "opencv2/opencv.hpp"
  8. #include "MOG_BGS.h"
  9. #include <iostream>
  10. #include <cstdio>
  11. using namespace cv;
  12. using namespace std;
  13. int main(int argc, char* argv[])
  14. {
  15. Mat frame, gray, mask;
  16. VideoCapture capture;
  17. capture.open("video.avi");
  18. if (!capture.isOpened())
  19. {
  20. cout<<"No camera or video input!\n"<<endl;
  21. return -1;
  22. }
  23. MOG_BGS Mog_Bgs;
  24. int count = 0;
  25. while (1)
  26. {
  27. count++;
  28. capture >> frame;
  29. if (frame.empty())
  30. break;
  31. cvtColor(frame, gray, CV_RGB2GRAY);
  32. if (count == 1)
  33. {
  34. Mog_Bgs.init(gray);
  35. Mog_Bgs.processFirstFrame(gray);
  36. cout<<" Using "<<TRAIN_FRAMES<<" frames to training GMM..."<<endl;
  37. }
  38. else if (count < TRAIN_FRAMES)
  39. {
  40. Mog_Bgs.trainGMM(gray);
  41. }
  42. else if (count == TRAIN_FRAMES)
  43. {
  44. Mog_Bgs.getFitNum(gray);
  45. cout<<" Training GMM complete!"<<endl;
  46. }
  47. else
  48. {
  49. Mog_Bgs.testGMM(gray);
  50. mask = Mog_Bgs.getMask();
  51. morphologyEx(mask, mask, MORPH_OPEN, Mat());
  52. erode(mask, mask, Mat(7, 7, CV_8UC1), Point(-1, -1));  // You can use Mat(5, 5, CV_8UC1) here for less distortion
  53. dilate(mask, mask, Mat(7, 7, CV_8UC1), Point(-1, -1));
  54. imshow("mask", mask);
  55. }
  56. imshow("input", frame);
  57. if ( cvWaitKey(10) == 'q' )
  58. break;
  59. }
  60. return 0;
  61. }

相关文章