背景建模gmm + 追踪sort原理及基于opencv的c++代码实现

目录

GMM(高斯混合模型)

GMM实现过程

EM算法

sort原理

gmm+sort 基于opencv的 c++工程代码

        最近做高空抛物项目，因采用深度学习的方法不适用，看了帧差法，光流法与前背景建模的方式，遂写下前背景建模的方法，也对比了比较通用的gmm和vibe的效果，于是将gmm算法原理+sort追踪及基于opencv的c++代码写上，附上简单的工程demo，工程中提供了关联的少量几个opencv动态库，如有需要可以拿来即用。

 测试效果如下图所示：

图1  gmm+sort 测试示意图

 注：绿框为检测到的前景框（原本是检测到前景物体分割的像素，后采用最小外接矩形将其轮廓刻画出来），中心红点为检测框的中心点，绿点为sort追踪时kalman的预测点，53*87为矩形框的像素宽高，绿线为sort追踪时kalman预测的中心点（取前50帧的预测中心点做的连线）

GMM(高斯混合模型)

         GMM 模型由多个高斯模型线性叠加混合而成，常用于聚类，每个高斯模型可以称为一个component，component的个数可以认为是类别的数量，即将原样本数据聚类为指定的类别数，其描述的为数据本身存在的一种分布。如下图示，图2为单高斯模型示意图，图3为混合高斯模型示意图，从图1，2对比可以发现，单高斯模型的中心处数据较少，一般远离椭圆中心的概率密度越小，表征模型的高斯分布越不好，同样预测效果也较差。而混合高斯模型更能够将样本数据分割开，使更多数据分布在椭圆中心，数据呈现高斯分布状态。

   

                     图2 单个高斯模型示意图                              图3 混合高斯模型示意图

 图片源自：林立民爱洗澡博文

对于单高斯混合模型其公式原理如下：

这是一个概率函数， x是维度为d的列向量，u是模型期望，即样本数据的均值，σ是模型方差。通过公式可以判断出一个随机样本x是否属于某个类别A。对于每个类别都有自己的均值u和方差σ，将任意一个样本x代入上式，当概率大于我们设定的阈值时就认为x属于其中某个类别A。

对于多个高斯模型混合的概率公式如下：

其中各参数表示含义为：

k：高斯模型的个数；

πk：为样本数据x属于第k个高斯模型的概率，所有的πk之和为1；

p(x，μk ， Σk )：为第k个高斯模型的概率密度，其中μk 为均值，Σk为方差。

混合高斯数据分布的示意图，如下图所示：

图4  高斯混合模型数据分布示意图 

GMM实现过程

       1. 首先设置初始化高斯模型的模型个数，随机初始化每个高斯模型的均值,方差，权值，即初始化模型矩阵参数。

        2.计算每个样本数据属于每个高斯模型的概率（后验概率，越靠近高斯分布中心，概率越大），选择基于视频中的前N帧用来训练背景。对每一个像素而言，建立其模型个数最大GMM_MAX_COMPONT个高斯的GMM模型。当第一个像素出现时，随机初始化均值，方差等参数值。

       3.计算 α,μ,Σ 参数使得数据有最大概率，使用数据点概率的加权来计算，更新这些参数，权重就是数据点属于该高斯模型的概率。当后面来的像素值，与前面已有的高斯的均值，方差比较，基于后验概率判断是否属于该高斯模型，进而对高斯模型参数进行再次更新；

对于参数的评估通过EM算法（最大期望算法），即在概率模型中寻找参数，使相应的后验概率最大化。

EM算法

① 随机初始化模型参数θ的初始值θ0；

②E步骤：估计隐藏变量的概率分布期望函数；

②M步骤：极大化期望函数，重新估计分布参数。

重复②、③步骤，直至 已经收敛，则算法结束，输出最终模型参数θ。

其原理如下：

基于opencv的c++的程序代码实现：

#include "opencv2/opencv.hpp"  
#include "opencv2/video/background_segm.hpp"
#include <opencv2/core/core.hpp>   
#include <opencv2/highgui/highgui.hpp>  

Ptr<BackgroundSubtractorMOG2> bg = cv::createBackgroundSubtractorMOG2(history, varThreshold, detectShadows)

其中参数history为用于训练背景的帧数，默认帧数为500帧，varThreshold为方差阈值，用于判断当前像素是前景还是背景，默认为16。detectShadows为bool值，表示是否检测影子，true为检测，false为不检测，检测影子会增加程序时间复杂度，默认为false；

对于检测到的前景物体采用kalman进行预测

sort原理

        SORT算法采用卡尔曼滤波和匈牙利算法来处理运动预测和数据关联这两个分量。检测、当前帧到未来帧的转换、当前帧与对象的关联、管理追踪目标的生命周期。 

kalman预测模型，首先随机初始化将目标的标识等信息传播到下一帧。即采用一个匀速模型来近似物体每一帧的位移，这个模型独立于其他物体。每个目标的状态被建模为：

其中，u，v 代表预测目标中心的横纵x、y坐标；

s、r表示bounding box的面积和长宽比；

u^,v^,s^表示对下一帧的预测中心坐标和检测框面积。bounding box用于更新目标状态，其中的速度分量使用卡尔曼滤波进行求解。如果没有和目标关联的检测框，就使用线性的预测模型而不需要修正。

匈牙利滤波的级联匹配在这儿不做讲解，下一章会单独再做分析。

gmm+sort 基于opencv的 c++工程代码

#include <iostream> 
#include "opencv2/opencv.hpp"  
#include "opencv2/video/background_segm.hpp"
#include <opencv2/core/core.hpp>   
#include <opencv2/highgui/highgui.hpp>  

using namespace std;
using namespace cv;

vector<Point> mousev, kalmanv;
KalmanFilter KF;
Mat_<float> measurement(2, 1);   //模板矩阵
Mat_<float> state(4, 1); // (x, y, Vx, Vy)  
int incr = 0;
string num2str(int i)
{
	stringstream ss;
	ss << i;
	return ss.str();
}

void initKalman(float x, float y)
{ 
	KF.init(4, 2,   0);   //卡尔曼的4个动量参数和2个测量参数

	measurement = Mat_<float>::zeros(2, 1);
	measurement.at<float>(0, 0) = x;
	measurement.at<float>(0, 0) = y;

	KF.statePre.setTo(0);
	KF.statePre.at<float>(0, 0) = x;
	KF.statePre.at<float>(1, 0) = y;
	KF.statePost.setTo(0);   // 设置部分值为指定值
	KF.statePost.at<float>(0, 0) = x;
	KF.statePost.at<float>(1, 0) = y;

	setIdentity(KF.transitionMatrix);  //初始化一个缩放矩阵 s：分配给对角线的元素值
	setIdentity(KF.measurementMatrix);
	setIdentity(KF.processNoiseCov, Scalar::all(.005)); //adjust this for faster convergence - but higher noise  快速收敛，噪音较大
	setIdentity(KF.measurementNoiseCov, Scalar::all(1e-1));
	setIdentity(KF.errorCovPost, Scalar::all(.1));
}

Point kalmanPredict()
{
	Mat prediction = KF.predict();  //计算预测状态
	Point predictPt(prediction.at<float>(0), prediction.at<float>(1));

	KF.statePre.copyTo(KF.statePost);
	KF.errorCovPre.copyTo(KF.errorCovPost);

	return predictPt;
}

Point kalmanCorrect(float x, float y)
{
	measurement(0) = x;
	measurement(1) = y;
	Mat estimated = KF.correct(measurement);  //更新预测状态
	Point statePt(estimated.at<float>(0), estimated.at<float>(1));  //更新后的值再传入
	return statePt;
}

int main()
{
	Mat frame, thresh_frame;
	int op = 0;
	vector<Mat> channels;
	
	// 设置视频及保存视频路径
	string save_path = "../video_save/test_result.mp4";
	string video_path = "../test.mp4";
	VideoCapture capture(video_path);

	// 添加，保存视频
	VideoWriter writer;
	Size size = Size(capture.get(CAP_PROP_FRAME_WIDTH), capture.get(CAP_PROP_FRAME_HEIGHT));
	int codec = VideoWriter::fourcc('m', 'p', '4', 'v');
	int rate = capture.get(cv::CAP_PROP_FPS);
	std::cout << "rate" << rate;	
	writer.open(save_path, codec, rate, size, true);

	int lo = 0;
	map<string,int> dict;
	// vector<int> vec1;
 
	if (!capture.isOpened())  // 成功打开视频 capture.isOpened() 返回true
		cerr << "Problem opening video source" << endl;  
	vector<Vec4i> hierarchy;  
	vector<vector<Point> > contours;   //vector 轮廓点储存再contours

	// 设置初始的背景，前景
	Mat back,img;
	Mat fore;
	int history=300;
	double varThreshold=16;
	bool detectShadows=false;

	Ptr<BackgroundSubtractorMOG2> bg = createBackgroundSubtractorMOG2(history, varThreshold, detectShadows);  //创建MOG背景减法器； history：用于训练背景的帧数，默认帧数为500帧； varThreshold:方差阈值，默认16.判断当前像素是前景还是背景

	int track = 0;
	bool update_bg_model = true;

	mousev.clear();
	kalmanv.clear();

	initKalman(0, 0);

	while (capture.grab())   // capture.grab() 获取视频帧
	{
		Point s, p;   // opencv 2维点模板s，p

		capture.retrieve(frame);
		bg->apply(frame, fore, -1);  //高斯混合，frame 原图（当前帧），fore：前景；update_bg_model为true(-1)（自动更新），为false：（0）不更新
		bg->getBackgroundImage(back);   // 背景图片结果
		erode(fore, fore, Mat());    //腐蚀源图像，取最小值的像素邻域的形状 
		dilate(fore, fore, Mat());  // 膨胀；

		normalize(fore, fore, 0, 1., NORM_MINMAX);
		threshold(fore, fore, .9, 1., THRESH_BINARY);  // 阈值 (阈值类型THRESH_BINARY，当大于设定值0.9，设为最大值，否则设为0)

		split(frame, channels);  // 分离通道
		add(channels[0], channels[1], channels[1]); 
		subtract(channels[2], channels[1], channels[2]);  //c2 = c2 - c1 -c0
		// threshold(channels[2], thresh_frame, 50, 255, CV_THRESH_BINARY);
		threshold(channels[2], thresh_frame, 50, 255, THRESH_BINARY);  // 将减完的c2通道根据阈值投影到（50，255）的像素间
		medianBlur(thresh_frame, thresh_frame, 7);  // 中值滤波；孔径线性尺寸，可以修改数值，来调整

		findContours(fore, contours, hierarchy, RETR_EXTERNAL, CHAIN_APPROX_SIMPLE, Point(0, 0)); // RETR_EXTERNAL：只检测最外围轮廓 ；RETR_LIST：检测所有轮廓（内，外）
		                                                                                    	 // CHAIN_APPROX_SIMPLE：保存物体边界上所有连续的轮廓点到contours向量内
		
																								 // Point(0, 0)：Point偏移量，所有的轮廓信息相对于原始图像对应点的偏移量，相当于在每一个检测出的轮廓点上加上该偏移量
																								 // contours 保存向量每组point的点集轮廓
		vector<vector<Point>> contours_poly(contours.size());
		vector<Rect> boundRect(contours.size());

		//  Get the mass centers:  
		for (size_t i = 0; i < contours.size(); i++)
		{
			approxPolyDP(Mat(contours[i]), contours_poly[i], 3, true); // 将连续光滑的线条曲折化，即采用最小多边形包围物体（类似外接圆）； 第三个参数 epsilon越小，折线的形状越“接近”曲线
			boundRect[i] = boundingRect(Mat(contours_poly[i]));  // 用最小矩形将外围曲线包裹起来，返回最小矩形的(x,y,w,h)
		}

		p = kalmanPredict();
		// cout << "kalman prediction: " << p << "p.x" << p.x << "p.y " << p.y << endl;  
		mousev.push_back(p);
		string text;

		for (size_t i = 0; i < contours.size(); i++)
		{
			int point_value = 0;
			
			if (contourArea(contours[i]) > 500)       // 前景轮廓的面积大于1000就绘出相应的图形  // contours[i] 代表被检测物体的轮廓点
			{
				cout << "contourArea(contours[i]" << contourArea(contours[i]) << std::endl;
				rectangle(frame, boundRect[i].tl(), boundRect[i].br(), Scalar(0, 255, 0), 2, 8, 0);   // 绘制物体的矩形框

				Point pt = Point(boundRect[i].x, boundRect[i].y+boundRect[i].height+15);
				Point pt1 = Point(boundRect[i].x, boundRect[i].y - 10);
				Point center = Point(boundRect[i].x + (boundRect[i].width / 2), boundRect[i].y + (boundRect[i].height / 2));
				circle(frame, center, 8, Scalar(0, 0, 255), -1, 1, 0);                               //绘制物体中心点
				circle(frame, Point(p.x, p.y), 8, Scalar(0, 255, 0), -1, 1, 0);      //绘制kalman 预测点
				Scalar color = CV_RGB(255, 0, 0);
				text = num2str(boundRect[i].width) + "*" + num2str(boundRect[i].height);
				putText(frame,"object", pt, cv::FONT_HERSHEY_DUPLEX, 1.0f, color);               // 绘制框下面 物体名称
				putText(frame,text, pt1, cv::FONT_HERSHEY_DUPLEX, 1.0f, color);                 // 绘制框上面 物体尺寸

				s = kalmanCorrect(center.x, center.y);   //s 卡尔曼更新后的中心坐标
				// drawCross(frame, s, Scalar(255, 255, 255), 5);

				int center_y1 = boundRect[i].y + (boundRect[i].height / 2);
				for (int i = mousev.size() - 50; i < mousev.size() - 1; i++)  //  mousev 经过卡尔曼预测的中心点变化的整个过程曲线,选取前50此kalman预测中心进行轨迹连线
				{
					line(frame, mousev[i], mousev[i + 1], Scalar(0, 255, 0), 1);
				}
			}
		}
		
		// 添加，保存视频
		writer << frame;

		// 视频展示
		// imshow("Video", frame);
	}
	return 0;

	
}

对代码函数等细节在上述进行了注释，方便读者加深对原理和过程的理解，此外可以在个人github工程上拉取：GitHub - ywysf120806/gmm-sort: Background Modeling and Multiple Object Tracking

在工程目录下执行以下代码

1. cd gmm-sort

2.mkdir build

3. cmake ..

4. make

5. ./gmm_sort