混合高斯背景建模程序分析

（转）混合高斯背景建模程序分析

 

标签： 杂谈

分类： 图像处理

将高斯建模改成了用一个亮度分量信息建立，但是发现，修改那个权值，还有那个多少个高斯分布，好像对检查没有什么影响一样！

#include <stdio.h>
#include <cv.h>
#include <cxcore.h>
#include <highgui.h>
#include <cvaux.h>//必须引此头文件

void main( )
{
//参数初始化定义
IplImage* pFrame = NULL;
IplImage* pFrImg = NULL;
IplImage* pBkImg = NULL;
CvCapture* pCapture = NULL;
IplImage* origin_rgb = NULL ;//定义rgb空间的存储
IplImage* origin_ycc = NULL ;//定义转换成YCrCb空间的存储
IplImage* lumi = NULL ;//定义亮度分量的存储空间

//定义窗口
cvNamedWindow("lumi",1);//定义显示窗口的名字，显示原始的视频
cvMoveWindow("lumi",30,0);//定义显示窗口的位置
cvNamedWindow("background",1);//显示经过转换格式的视频
cvMoveWindow("background",360,0);
cvNamedWindow("foreground",1);//显示经过亮度提取的视频
cvMoveWindow("foreground",690,0);

//读取一帧视频文件作为初始化
pCapture = cvCaptureFromFile("video.long.xvid.avi") ;
pFrame = cvQueryFrame(pCapture);
int i ;
for (i=0;i<2;i++)
{
pFrame = cvQueryFrame(pCapture);
}
//RGB转换成亮度
origin_rgb = cvCreateImage(cvSize(pFrame->width, pFrame->height), IPL_DEPTH_8U, 3 );
origin_ycc = cvCreateImage(cvSize(pFrame->width, pFrame->height), IPL_DEPTH_8U, 3 );
lumi = cvCreateImage(cvSize(pFrame->width, pFrame->height), IPL_DEPTH_8U, 1 );
pFrImg = cvCreateImage(cvSize(pFrame->width, pFrame->height), IPL_DEPTH_8U, 1 );
pBkImg = cvCreateImage(cvSize(pFrame->width, pFrame->height), IPL_DEPTH_8U, 1 );

//origin_rgb = cvCloneImage(pFrame) ;//或者 cvCopy(pFrame,origin_rgb,NULL) ;
cvCopy(pFrame,origin_rgb,NULL) ;//复制视频
cvCvtColor(origin_rgb,origin_ycc,CV_BGR2YCrCb) ; //实现视频格式转换
cvSplit(origin_ycc,lumi,NULL,NULL,NULL); //获取亮度分量

//为高斯模型设置初时参数
CvGaussBGStatModelParams* params = new CvGaussBGStatModelParams;
params->win_size = 50;
params->n_gauss = 3;
params->bg_threshold = 0.7;
params->std_threshold = 3.5;
params->minArea = 15;
params->weight_init = 0.333;
params->variance_init = 30;

CvBGStatModel* bgModel = cvCreateGaussianBGModel(lumi,params);

int key=-1;
while(key != 'q')
{
//获取下一帧视频
pFrame = cvQueryFrame(pCapture);
cvCopy(pFrame,origin_rgb,NULL) ;//复制视频
cvCvtColor(origin_rgb,origin_ycc,CV_BGR2YCrCb) ; //实现视频格式转换
cvSplit(origin_ycc,lumi,NULL,NULL,NULL); //获取亮度分量
if( !pFrame )
break;

//更新高斯模型
cvUpdateBGStatModel(lumi,bgModel);
pFrImg = bgModel->foreground ;//前景图象
pBkImg = bgModel->background ; //背景图象

//将图象倒转过来
pBkImg->origin = 1 ;
pFrImg->origin = 1 ;
lumi->origin = 1 ;

//显示结果
cvShowImage("lumi",lumi);
cvShowImage("background",pBkImg);
cvShowImage("foreground",pFrImg);
key = cvWaitKey(10);
}
// cvWaitKey(0) ;//窗口的回调函数，必须要的，不然窗口的显示会不正常
//释放窗口内存资源
cvDestroyWindow("lumi");
cvDestroyWindow("background");
cvDestroyWindow("foreground");

//释放图象占用的内存资源
cvReleaseImage(&lumi);
cvReleaseImage(&pBkImg);
cvReleaseImage(&pFrImg);
cvReleaseCapture(&pCapture);
cvReleaseBGStatModel( &bgModel );
}

创建高斯背景模型

cvCreateGaussianBGModel( IplImage* first_frame, CvGaussBGStatModelParams* parameters )
{
//CvGaussBGModel在cvaux.h中有定义，是一个结构体
CvGaussBGModel* bg_model = 0;

CV_FUNCNAME( "cvCreateGaussianBGModel" );//在cxerror.h定义，定义cvFuncName宏变量
//cvFuncName定义为和函数名称相同cvCreateGaussianBGModel

__BEGIN__;//开始处理（是必须接在这个CV_FUNCNAME之后的）

double var_init;
CvGaussBGStatModelParams params;//定义初始化变量，在cvaux.h中定义的结构体CvGaussBGStatModelParams
int i, j, k, n, m, p;

//init parameters
if( parameters == NULL )
{
params.win_size = CV_BGFG_MOG_WINDOW_SIZE;//CV_BGFG_MOG_WINDOW_SIZE＝200，和学习率的关系1/200＝0.005
params.bg_threshold = CV_BGFG_MOG_BACKGROUND_THRESHOLD;//CV_BGFG_MOG_BACKGROUND_THRESHOLD＝0.7（判断是否为背景点的阈值）
params.std_threshold = CV_BGFG_MOG_STD_THRESHOLD;//CV_BGFG_MOG_STD_THRESHOLD＝2.5（标准阈值）
params.weight_init = CV_BGFG_MOG_WEIGHT_INIT;//CV_BGFG_MOG_WEIGHT_INIT＝0.05（权值）
params.variance_init = CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT;//CV_BGFG_MOG_SIGMA_INIT＝30（方差）
params.minArea = CV_BGFG_MOG_MINAREA;//CV_BGFG_MOG_MINAREA＝15.f（这个不知道？）
params.n_gauss = CV_BGFG_MOG_NGAUSSIANS;//CV_BGFG_MOG_NGAUSSIANS＝5（高斯模型数量）
}
else
{
params = *parameters;
}
//CV_IS_IMAGE在cxtypes.h中定义，在这里估计是判断是否有读入图象帧
//CV_StsBadArg＝-5，代表函数有问题，或者输入的参数有问题
if( !CV_IS_IMAGE(first_frame) )
CV_ERROR( CV_StsBadArg, "Invalid or NULL first_frame parameter" );

//CV_CALL在cxerror.h中有定义，这里用来确认一下调用是否正确
CV_CALL( bg_model = (CvGaussBGModel*)cvAlloc( sizeof(*bg_model) ));
memset( bg_model, 0, sizeof(*bg_model) );
bg_model->type = CV_BG_MODEL_MOG;//这个在cvGaussBGModel中定义的CV_BG_STAT_MODEL_FIELDS()函数中都有，存在type，release，update，foreground，background等
bg_model->release = (CvReleaseBGStatModel)icvReleaseGaussianBGModel;
bg_model->update = (CvUpdateBGStatModel)icvUpdateGaussianBGModel;

bg_model->params = params;

//分配存储空间
CV_CALL( bg_model->g_point = (CvGaussBGPoint*)cvAlloc(sizeof(CvGaussBGPoint)*
((first_frame->width*first_frame->height) + 256)));//这个是参与的点数，以及存放这些点需要的空间
CV_CALL( bg_model->background = cvCreateImage(cvSize(first_frame->width,
first_frame->height), IPL_DEPTH_8U, first_frame->nChannels));//给背景分配存储空间
CV_CALL( bg_model->foreground = cvCreateImage(cvSize(first_frame->width,
first_frame->height), IPL_DEPTH_8U, 1));//给前景分配存储空间

CV_CALL( bg_model->storage = cvCreateMemStorage());//分配存储空间

//初始化

var_init = 2 * params.std_threshold * params.std_threshold;
CV_CALL( bg_model->g_point[0].g_values =
(CvGaussBGValues*)cvAlloc( sizeof(CvGaussBGValues)*params.n_gauss*
(first_frame->width*first_frame->height + 128)));//这个是给g_value分配足够的存储空间

//程序说明
//g_values[0],g_values[1],g_values[2]存放3个高斯混合模型的变量
//g_values[].weight(权重) g_values[].mean[](均值) g_values[].variance[](方差)
//具体安排是每一个象素点都有3个模型，
//然后每一个象素点的三个模型
//模型0的权重为1，方差为2倍的标准差的平方，均值为当前象素点的值
//模型1的权重为0，方差为2倍的标准差的平方，均值为0
//模型2的权重为0，方差为2倍的标准差的平方，均值为0
//g_point指的是参加高斯背景建模的象素点的个数
for( i = 0, p = 0, n = 0; i < first_frame->height; i++ )
{
for( j = 0; j < first_frame->width; j++, n++ )
{
bg_model->g_point[n].g_values =
bg_model->g_point[0].g_values + n*params.n_gauss;
bg_model->g_point[n].g_values[0].weight = 1; //the first value seen has weight one
bg_model->g_point[n].g_values[0].match_sum = 1;
for( m = 0; m < first_frame->nChannels; m++)
{
bg_model->g_point[n].g_values[0].variance[m] = var_init;
bg_model->g_point[n].g_values[0].mean[m] = (unsigned char)first_frame->imageData[p + m];
}
for( k = 1; k < params.n_gauss; k++)
{
bg_model->g_point[n].g_values[k].weight = 0;
bg_model->g_point[n].g_values[k].match_sum = 0;
for( m = 0; m < first_frame->nChannels; m++){
bg_model->g_point[n].g_values[k].variance[m] = var_init;
bg_model->g_point[n].g_values[k].mean[m] = 0;
}
}
p += first_frame->nChannels;
}
}

bg_model->countFrames = 0;

__END__;

if( cvGetErrStatus() < 0 )
{
CvBGStatModel* base_ptr = (CvBGStatModel*)bg_model;

if( bg_model && bg_model->release )
bg_model->release( &base_ptr );
else
cvFree( &bg_model );
bg_model = 0;
}

return (CvBGStatModel*)bg_model;
}

这整个函数就是对结构体cvGaussBGModel里面的参数：CV_BG_STAT_MODEL_FIELDS(),params，g_point，countFrames赋值，实际上也是可以把其复制过来，自己修改初时的参数。//对每个象素点进行匹配
static int icvMatchTest( double* src_pixel, int nChannels, int* match,

const CvGaussBGPoint* g_point,
const CvGaussBGStatModelParams *bg_model_params )
{
int k;
int matchPosition=-1; //这个匹配的位置只能是0，1，2（在假设3个高斯的前提下）
for ( k = 0; k < bg_model_params->n_gauss; k++) match[k]=0;

for ( k = 0; k < bg_model_params->n_gauss; k++)
{
double sum_d2 = 0.0;
double var_threshold = 0.0;
for(int m = 0; m < nChannels; m++)
{ //src_pixel的值是由pixel传送
//pixel = (uchar)curr_frame->imageData[p+k]
//即是如果是单通道或者是多通道，都是把一个象素点的所有分量全部读出来进行比较
//g_point->g_values[].mean[]就是在cvCreateGaussianBGModel对bg_model进行赋值的结果
double d = g_point->g_values[k].mean[m]- src_pixel[m];
sum_d2 += (d*d);
var_threshold += g_point->g_values[k].variance[m];
} //difference < STD_LIMIT*STD_LIMIT or difference**2 < STD_LIMIT*STD_LIMIT*VAR
//我觉得这里有问题，我觉得应该是：
//var_threshold = bg_model_params->std_threshold*var_threshold*var_threshold;
var_threshold = bg_model_params->std_threshold*bg_model_params->std_threshold*var_threshold;
if(sum_d2 < var_threshold)
{
match[k] = 1;
matchPosition = k;
break;
}
}

return matchPosition;
}
觉得opencv本身携带的模型并不好，感觉可以根据他的模块改一下程序，比如cvCreateGaussianModel和icvMatchTest都可以在进行改进。

记录下来现在学习的想法，不然以后忘记了，呵呵。

将来看看自己的博客，就是自己学习的历程，然后在整理，比较容易。

static void icvUpdateFullWindow( double* src_pixel, int nChannels, int* match,

CvGaussBGPoint* g_point,
const CvGaussBGStatModelParams *bg_model_params )
{
const double learning_rate_weight = (1.0/(double)bg_model_params->win_size);
for(int k = 0; k < bg_model_params->n_gauss; k++)
{
//这个就是更新权重，w=(1-a)w+a*m (a:更新率，m是匹配，匹配就是1，不匹配就是0)
//可以改写：
// g_point->g_values[k].weight = g_point->(1-learning_rate_weight)*g_values[k].weight
// +learning_rate_weight*(double)match[k]
g_point->g_values[k].weight = g_point->g_values[k].weight +
(learning_rate_weight*((double)match[k] -
g_point->g_values[k].weight));
//更新方差，均值等参数
if(match[k])
{
//learning_rate_gaussian实际就是p（在更新方差和均值要用的参数）
//实际上可以这样设置：
//double learning_rate_gaussian = ((double)match[k]*learning_rate_weight)/(g_point->g_values[k].weight
double learning_rate_gaussian = (double)match[k]/(g_point->g_values[k].weight*
(double)bg_model_params->win_size);
for(int m = 0; m < nChannels; m++)
{
//均值和方差的更新公式
//u = (1-p)*u+p*x = u+p*(x-u)
//o = (1-p)*o*o+p*(x-u)*(x-u)
const double tmpDiff = src_pixel[m] - g_point->g_values[k].mean[m];
g_point->g_values[k].mean[m] = g_point->g_values[k].mean[m] +
(learning_rate_gaussian * tmpDiff);
g_point->g_values[k].variance[m] = g_point->g_values[k].variance[m]+
(learning_rate_gaussian*((tmpDiff*tmpDiff) - g_point->g_values[k].variance[m]));
}
}
}
}