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 opened
if (!capture.isOpened())
return 0;

// current video frame
cv::Mat frame, frameGray;
// foreground binary image
cv::Mat foreground;

cv::namedWindow("Extracted Foreground");

// The Mixture of Gaussian object
// used with all default parameters
cv::BackgroundSubtractorMOG mog;

bool stop(false);
// for all frames in video
while (!stop) {

// read next frame if any
if (!capture.read(frame))
break;

// update the background
// and return the foreground
cvtColor(frame, frameGray, CV_BGR2GRAY);
mog(frameGray,foreground,0.01);

// show foreground
cv::imshow("Extracted Foreground",foreground);
cv::imshow("image",frame);

// introduce a delay
// or press key to stop
if (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
.w = 0;
ppgmm[i*cols + j].gd
.u = 0;
ppgmm[i*cols + j].gd
.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;
else
ppgmm[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;
else
ppgmm[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;
b3d6
n >= 0; n--)
{
if(ppgmm[i*cols + j].gd
.w/ppgmm[i*cols + j].gd
.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
.sigma;
ppgmm[i*cols + j].gd
.sigma = ppgmm[i*cols + j].gd[n+1].sigma;
ppgmm[i*cols + j].gd[n+1].sigma = temp;
temp = ppgmm[i*cols + j].gd
.u;
ppgmm[i*cols + j].gd
.u = ppgmm[i*cols + j].gd[n+1].u;
ppgmm[i*cols + j].gd[n+1].u = temp;
temp = ppgmm[i*cols + j].gd
.w;
ppgmm[i*cols + j].gd
.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
.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
.w *= weightScale;
weightSum += ppgmm[i*cols + j].gd
.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年。。。不得不佩服这群人的专业和强大，大神就是大神。。。