背景建模或前景检测之GMM

博客：背景建模或前景检测之GMM

自己的工程代码在开源中国：GMM源代码+详细注释

本文通过OpenCV来实现一种前景检测算法——GMM，算法采用的思想来自论文[1][2][4]。在进行前景检测前，先对背景进行训练，对图像中每个背景采用一个混合高斯模型进行模拟，每个背景的混合高斯的个数可以自适应。然后在测试阶段，对新来的像素进行GMM匹配，如果该像素值能够匹配其中一个高斯，则认为是背景，否则认为是前景。由于整个过程GMM模型在不断更新学习中，所以对动态背景有一定的鲁棒性。最后通过对一个有树枝摇摆的动态背景进行前景检测，取得了较好的效果。

关键字：GMM,opencv，前景检测

前言

 

　　前景检测主要分为帧差法，平均背景法，光流法，前景建模法，背景非参数估计，背景建模法等。而本文要实现的方法属于背景建模法中的一种——GMM，也称混合高斯模型。

　　混合高斯模型最早在计算机视觉中的应用是Stauffer et al.[1],作者将其用来做前景检测，主要是用于视频监控领域，这个系统和稳定且有自学能力，能在户外环境跑16个多月。KaewTraKulPong et al.[2]将GMM的训练过程做了改进，将训练过程分为2步进行，前L帧采用EM算法进行权值，均值，方差更新，后面的过程就采用[1]中的方法进行更新，取得了更好的检测效果。Zivkovic et al.在[3]中对GMM理论做了全面的论述，使得GMM理论的使用不仅金限于计算机视觉领域。并且该作者在[4]将该理论进一步具体到背景减图的前景检测中来，即加入了参数估计的先验知识，取得了很好的效果和稳定性。

　　最近几年陆续有学者对GMM的背景见图中的应用做了更深一步的研究，其代表性贡献见论文[5][6]。

 

  实现过程

 

　　本文中主要是根据[2][4]中提出的算法，采用其中的更新方差来根性模型中的3个参数，最后结合用opencv基本底层函数实现该算法。其主要过程如下：

1. 首先将每个高斯的均值,方差，权值都设置为0，即初始化个模型矩阵参数。
2. 采用视频中的T帧用来训练GMM模型。对每一个像素而言，建立其模型个数最大GMM_MAX_COMPONT个高斯的GMM模型。当第一个像素来，单独为其在程序中设置好其固定的初始均值，方差，并且权值设置为1。
3. 非第一帧训练过程中，当后面来的像素值时，与前面已有的高斯的均值比较，如果该像素点的值与其模型均值差在3倍的方差内，则任务属于该高斯。此时用如下方程进行更新：

　　　　

　　4. 当到达训练的帧数T后，进行不同像素点GMM个数自适应的选择。首先用权值除以方差对各个高斯进行从大到小排序，然后选取最前面B个高斯，使

　　

　　这样就可以很好的消除训练过程中的噪声点。

　　5. 在测试阶段，对新来像素点的值与B个高斯中的每一个均值进行比较，如果其差值在2倍的方差之间的话，则认为是背景，否则认为是前景。并且只要其中有一个高斯分量满足该条件就认为是前景。前景赋值为255，背景赋值为0。这样就形成了一副前景二值图。

　　6. 由于前景二值图中含有很多噪声，所以采用了形态学的开操作将噪声缩减到0，紧接着用闭操作重建由于开操作丢失的边缘部分的信息。消除了不连通的小噪声点。

　　上面是该算法实现的大概流程，但是当我们在具体编程时，还是有很多细节的地方需要注意，比如有些参数值的选择。在这里，经过试验将一些常用的参数值声明如下：

1. 3个参数的值的更新方差中，取其中的学习率为0.005。也就是说T等于200。
2. 定义每个像素的最大混合高斯个数取7。
3. 取视频的前200帧进行训练。
4. 取Cf为0.3。即满足权值大于0.7的个数为自适应高斯的个数B。
5. 训练过程中，需要新建立一个高斯时，其权值取值设为与学习率大小值相等，即0.005。
6. 训练过程中，需要新建立一个高斯时，取该高斯的均值为该输入像素值大小本身。方差为15。
7. 训练过程中，需要新建立一个高斯时，取该高斯的方差15。

 

   程序流程框图

 

　　该工程的流程框图如下图所示：

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　

　　　　　　　　　　　　　　　　　　　　　　　　　　

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　   　　　　　　　　　　

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　试验结果

 

　　本次试验的数据为摇摆的树枝作为背景，Waving Trees，其来源网址为：http://research.microsoft.com/en-us/um/people/jckrumm/WallFlower/TestImages.htm 由于该数据是286张bmp格式的图片，所以用的前200张图片来作为GMM参数训练，后186张作为测试。训练的过程中树枝被很大幅度的摆动，测试过程中有行人走动，该行人是需要迁就检测的部分。

　　下图为训练过程中动态背景截图

　　

 

　　由上图可以看出，树枝在不断摇摆，可见其背景是动态的。

　　没有形态学处理的结果如下：

　　 

     下图是有简单形态学的试验结果：

　　 

      (上面2幅图的左边为测试原始图片，右图为检测效果图)

 

总结

 

　　通过本次试验，不仅学习到了GMM的相关理论知识，以及背景减图法在前景检测中的应用。更重要的时，对opencv在计算机视觉中的使用更进一步的熟悉了。另外对于目标检测的难点有了更深一层的理解。

 

参考文献

 

1. Stauffer, C. and W. E. L. Grimson (1999). Adaptive background mixture models for real-time tracking, IEEE.
2. KaewTraKulPong, P. and R. Bowden (2001). An improved adaptive background mixture model for real-time tracking with shadow detection.
3. Zivkovic, Z. and F. van der Heijden (2004). "Recursive unsupervised learning of finite mixture models." Pattern Analysis and Machine Intelligence, IEEE Transactions on 26(5): 651-656.
4. Zivkovic, Z. and F. van der Heijden (2006). "Efficient adaptive density estimation per image pixel for the task of background subtraction." Pattern recognition letters 27(7): 773-780.
5. Bouzerdoum, A., A. Beghdadi, et al. (2010). "On the analysis of background subtraction techniques using Gaussian mixture models."
6. Lin, H. H., J. H. Chuang, et al. (2011). "Regularized background adaptation: a novel learning rate control scheme for Gaussian mixture modeling." Image Processing, IEEE Transactions on 20(3): 822-836.

 

                                                                                                         后续改进

 

　　需要改进的地方：

　　1. 程序运行的速度太慢，很多参数都是浮点数，每个像素都要建立一个gmm，gmm个数本身比较多，所以训练过程中速度比较慢，代码需要优化。

　　2. 最后生成的前景图需要用连通域处理算法进行修整，即需要形态学操作，然后找出连通域大小满足要求的轮廓，用多边形拟合来进行处理。这种算法在《learning opencv》一书中有提到。后续有时间加入该算法，效果会好很多的。

 

MOG_BGS.h

[cpp]  view plain  copy

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 plain  copy

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 plain  copy

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!\n"<<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. }