opencv中的GMM（混合高斯分布）算法原理及C++实现（BackgroundSubtractorMOG）

1.opencv中的GMM算法

GMM（Gaussian Mixture Model）是一种经典的背景提取算法，opencv中也把它引入并封装为算法类。使用opencv2413版本时，通过BackgroundSubtractorMOG类即可调用。这里先给出一个调用例子和算法效果，代码如下所示。

#include <iostream>#include <opencv2/core/core.hpp>#include <opencv2/imgproc/imgproc.hpp>#include <opencv2/highgui/highgui.hpp>#include <opencv2/video/background_segm.hpp>int main(){// Open the video file    cv::VideoCapture capture("768X576.avi");// check if video successfully openedif (!capture.isOpened())return 0;// current video framecv::Mat frame, frameGray; // foreground binary imagecv::Mat foreground;cv::namedWindow("Extracted Foreground");// The Mixture of Gaussian object// used with all default parameterscv::BackgroundSubtractorMOG mog;bool stop(false);// for all frames in videowhile (!stop) {// read next frame if anyif (!capture.read(frame))break;// update the background// and return the foregroundcvtColor(frame, frameGray, CV_BGR2GRAY);mog(frameGray,foreground,0.01);// show foregroundcv::imshow("Extracted Foreground",foreground);cv::imshow("image",frame);// introduce a delay// or press key to stopif (cv::waitKey(10)>=0)stop= true;}return 0;}

其中“768X576.avi”是opencv源码文件中的一个包含多个行人的视频文件，算法效果如下图所示，可以看到行人被很好的用白色标注出来。

2.算法原理

BackgroundSubtractorMOG的实现在2413版本的opencv_2.4.13\opencv\sources\modules\video\src\bgfg_gaussmix.cpp中，cpp中讲到算法原理参考《An Improved Adaptive Background Mixture Model for Real-time Tracking with Shadow Detection》文章，但是同时说明了算法只实现了检测部分，并没有实现文章中的阴影检测部分。后来查阅一些其他资料之后，个人认为该算法的基本原理应该来自于，或者说更加符合文章《Adaptive-background-mixture-models-for-real-time-tracking》（这篇文章的下载和 翻译参考点击打开链接）。其原理这里简单描述如下。

1）背景

混合高斯模型的必要性：
（1）如果场景中光照条件稳定且背景固定，那么某一点处的像素值会比较稳定，考虑到噪声的存在，使用一个高斯分布即可比较好的描述该点像素值
（2）实际场景中，光照条件往往不稳定，且存在一些局部运动的背景物体（譬如摇摆的树枝），此时使用多个高斯分布表达比较合适
（3）实际场景中背景区域也可能实时改变，所以需要实时更新混合高斯分布参数的策略

2）算法步骤

文章《Adaptive-background-mixture-models-for-real-time-tracking》中描述的十分清楚了，精简版本的推荐参考https://wenku.baidu.com/view/273e56855f0e7cd1842536c4.html?from=search。不重复造轮子，这里把主要步骤截图如下：

这里把个人写代码实践过程中的一些经验分享下：

（1）上图中的“排序及删减”有些不确切，实际过程中不会“删减”N之后索引的高斯分布，而是作为前景点判断的一种条件

（2）计算过程中前景点的判断标准包括与当前高斯模型集合都不匹配的点，和虽然匹配上了但是匹配的模型索引超出了由T计算得到的N；后者一般是场景中新出现没有多久的像素点

3.C++实现

算法步骤已经明确，后面就可以着手实现了，自己写的代码如下：

#include "core/core.hpp"    #include "highgui/highgui.hpp"    #include "imgproc/imgproc.hpp"#include "video/tracking.hpp"#include <opencv2/video/background_segm.hpp>#include <iostream>    #include <numeric>#include <vector>using namespace cv;    using namespace std; struct GaussianDistribution{float w;// 权重float u;// 期望float sigma;// 标准差};struct PixelGMM{int num;// 高斯模型最大数量int index;// 当前使用的高斯模型数量GaussianDistribution* gd;// 高斯模型数组指针};bool modelInit = false;const int GAUSSIAN_MODULE_NUMS = 5;// 模型数量const float ALPHA = 0.005;// 学习速率const float SIGMA = 30;// 标准差初始值const float WEIGHT = 0.05;// 高斯模型权重初始值const float T = 0.7;// 有效高斯分布阈值int rows, cols;PixelGMM* ppgmm;int main(){Mat image;Mat imageGray, imageFG, imageMog;cv::BackgroundSubtractorMOG mog;// 打开视频文件VideoCapture cap("768x576.avi");if(!cap.isOpened()){cout<<"cannot open avi file"<<endl;return -1;}double fps = cap.get(CV_CAP_PROP_FPS);// 获取图像帧率int pauseTime = (int)(1000.f/fps);namedWindow("video");while(1){if(!cap.read(image))break;cvtColor(image, imageGray, CV_BGR2GRAY);// 转化为灰度图处理mog(imageGray,imageMog,0.005);// 初始化各个像素的高斯分布if(!modelInit){/* 高斯分布参数分配空间*/rows = image.rows;cols = image.cols;ppgmm = (PixelGMM*)malloc(rows*cols*sizeof(PixelGMM));for(int i = 0; i < rows*cols; i++){ppgmm[i].num = GAUSSIAN_MODULE_NUMS;ppgmm[i].index = 0;ppgmm[i].gd = (GaussianDistribution*)malloc(GAUSSIAN_MODULE_NUMS*sizeof(GaussianDistribution));}/* 初始化高斯分布参数 */for(int i = 0; i < rows; i++){for(int j = 0; j < cols; j++){for(int n = 0; n < GAUSSIAN_MODULE_NUMS; n++){ppgmm[i*cols + j].gd[n].w = 0;ppgmm[i*cols + j].gd[n].u = 0;ppgmm[i*cols + j].gd[n].sigma = SIGMA;}}}imageFG.create(rows, cols, CV_8UC1);modelInit = true;}{// 结果图像初始化为全黑imageFG.setTo(Scalar(0));for(int i = 0; i < rows; i++){for(int j = 0; j < cols; j++){int kHit = -1;// 匹配高斯模型索引int gray = imageGray.at<unsigned char>(i,j);//判断是否属于当前高斯模型for(int m = 0; m < ppgmm[i*cols + j].index; m++){if(fabs(gray - ppgmm[i*cols + j].gd[m].u) < 2.5*ppgmm[i*cols + j].gd[m].sigma)// 满足分布{// 更新该高斯分布的标准差、期望、权值，标准差变小，权值增加if(ppgmm[i*cols + j].gd[m].w > 1)ppgmm[i*cols + j].gd[m].w = 1;elseppgmm[i*cols + j].gd[m].w = (1-ALPHA)*ppgmm[i*cols + j].gd[m].w + ALPHA*1;ppgmm[i*cols + j].gd[m].u = (1-ALPHA)*ppgmm[i*cols + j].gd[m].u + ALPHA*gray;if(ppgmm[i*cols + j].gd[m].sigma < SIGMA/2)ppgmm[i*cols + j].gd[m].sigma = SIGMA;elseppgmm[i*cols + j].gd[m].sigma = sqrt((1-ALPHA)*ppgmm[i*cols + j].gd[m].sigma*ppgmm[i*cols + j].gd[m].sigma + ALPHA*(gray - ppgmm[i*cols + j].gd[m].u)*(gray - ppgmm[i*cols + j].gd[m].u));// 若同一高斯分布被重复匹配多次，其标准差会一直下降，这里需要设置最低阈值// 根据w/sigma降序排序int n;for(n = m-1; n >= 0; n--){if(ppgmm[i*cols + j].gd[n].w/ppgmm[i*cols + j].gd[n].sigma <= ppgmm[i*cols + j].gd[n+1].w/ppgmm[i*cols + j].gd[n+1].sigma){float temp;temp = ppgmm[i*cols + j].gd[n].sigma;ppgmm[i*cols + j].gd[n].sigma = ppgmm[i*cols + j].gd[n+1].sigma;ppgmm[i*cols + j].gd[n+1].sigma = temp;temp = ppgmm[i*cols + j].gd[n].u;ppgmm[i*cols + j].gd[n].u = ppgmm[i*cols + j].gd[n+1].u;ppgmm[i*cols + j].gd[n+1].u = temp;temp = ppgmm[i*cols + j].gd[n].w;ppgmm[i*cols + j].gd[n].w = ppgmm[i*cols + j].gd[n+1].w;ppgmm[i*cols + j].gd[n+1].w = temp;}else{break;}}kHit = n + 1;// 匹配高斯分布的最终索引break;}else{// 没有被匹配到的高斯分布权重降低ppgmm[i*cols + j].gd[m].w *= (1-ALPHA);}}// 增加新的高斯分布，属于前景if(kHit == -1){// 需要去除影响最小的高斯分布if(ppgmm[i*cols + j].index == GAUSSIAN_MODULE_NUMS){ppgmm[i*cols + j].gd[ppgmm[i*cols + j].index - 1].sigma = SIGMA;ppgmm[i*cols + j].gd[ppgmm[i*cols + j].index - 1].u = gray;ppgmm[i*cols + j].gd[ppgmm[i*cols + j].index - 1].w = WEIGHT;}else{// 增加新的高斯分布ppgmm[i*cols + j].gd[ppgmm[i*cols + j].index].sigma = SIGMA;ppgmm[i*cols + j].gd[ppgmm[i*cols + j].index].u = gray;ppgmm[i*cols + j].gd[ppgmm[i*cols + j].index].w = WEIGHT;ppgmm[i*cols + j].index++;}}// 高斯分布权值归一化float weightSum = 0;for(int n = 0; n < ppgmm[i*cols + j].index; n++){weightSum += ppgmm[i*cols + j].gd[n].w;}float weightScale = 1.f/weightSum;// 根据T得到有效高斯分布的截止索引weightSum = 0;int kForeground = -1;for(int n = 0; n < ppgmm[i*cols + j].index; n++){ppgmm[i*cols + j].gd[n].w *= weightScale;weightSum += ppgmm[i*cols + j].gd[n].w;if(weightSum > T && kForeground < 0){kForeground = n + 1;}}// 属于前景点的判断条件if((kHit > 0&&kHit >= kForeground) || kHit == -1){imageFG.at<unsigned char>(i,j) = 255;}}}}imshow("imageMog", imageMog);imshow("video", imageGray);imshow("imageFG", imageFG);waitKey(10);}return 0;}

算法效果和opencv中自带算法效果类似，不过执行速度相差很多。。。

4.其他

（1）opencv中的BackgroundSubtractorMOG类算法流程和上述相同，不过opencv考虑了单通道图像和三通道图像两种情况，第二种情况中高斯集合中的每个高斯单元由三个高斯分布组成，分别对应图像的三个通道数据。

（2）《Adaptive-background-mixture-models-for-real-time-tracking》文章中作者以此算法为背景提取的基本步骤，之后结合连通域计算等操作，实现了包括车辆和行人的跟踪和分类等功能，这还是在1999年。。。不得不佩服这群人的专业和强大，大神就是大神。。。