OPENCV中混合高斯背景模型的实现

看高斯混合模型的论文的时候感觉一头雾水，主要是那些概率公式比较难懂，看了几遍论文和代码才有了一点感觉，以下是个人的对混合高斯背景模型代码的理解，后面会通过进一步理解不断对此贴修正.....

注：本人对该代码的理解也有很多不到位的地方，希望以后会慢慢矫正。程序中代码的实现与论文中公式有些出入，不过不影响对全局的把握和理解。

混合高斯背景模型理论的论文：点击下载

////////////////////////cvCreateGaussianBGModel///////////////////////////////////////////

CV_IMPL CvBGStatModel *cvCreateGaussianBGModel( IplImage*first_frame,CvGaussBGStatModelParams* parameters )

{

CvGaussBGModel* bg_model = 0;

CV_FUNCNAME( "cvCreateGaussianBGModel" );

__BEGIN__;

double var_init;

CvGaussBGStatModelParams params;

int i, j, k, m, n;

// init parameters

if( parameters == NULL )

{

params.win_size = CV_BGFG_MOG_WINDOW_SIZE; // 初始化阶段的帧数；用户自定义模型学 习率a=1/win_size;

params.bg_threshold = CV_BGFG_MOG_BACKGROUND_THRESHOLD;

params.std_threshold = CV_BGFG_MOG_STD_THRESHOLD;

params.weight_init = CV_BGFG_MOG_WEIGHT_INIT;

params.variance_init = CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT;
//方差

params.minArea = CV_BGFG_MOG_MINAREA;

params.n_gauss = CV_BGFG_MOG_NGAUSSIANS;
//高斯分布函数的个数

}

else

{

params = *parameters; //用户自定义参数

}

if( !CV_IS_IMAGE(first_frame) )

CV_ERROR( CV_StsBadArg, "Invalid or NULL first_frame parameter" );

CV_CALL( bg_model = (CvGaussBGModel*)cvAlloc( sizeof(*bg_model) ));

memset( bg_model, 0, sizeof(*bg_model) );

bg_model->type = CV_BG_MODEL_MOG; //CV_BG_MODEL_MOG为高斯背景模型

bg_model->release = (CvReleaseBGStatModel)icvReleaseGaussianBGModel;

bg_model->update = (CvUpdateBGStatModel)icvUpdateGaussianBGModel;

bg_model->params = params;

//prepare storages

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());

//initializing

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)));

for( i = 0, n = 0; i < first_frame->height; i++ )
//行

{

for( j = 0; j < first_frame->width; j++, n++ )
//列

{

const int p = i*first_frame->widthStep+j*first_frame->nChannels;

//以下几步是对第一个高斯函数做初始化

bg_model->g_point
.g_values = bg_model->g_point[0].g_values + n*params.n_gauss;

bg_model->g_point
.g_values[0].weight = 1; //权值赋为1

bg_model->g_point
.g_values[0].match_sum = 1;
//高斯函数被匹配的次数

for( m = 0; m < first_frame->nChannels; m++)

{

bg_model->g_point
.g_values[0].variance[m] = var_init;

//均值赋为当前像素的值

bg_model->g_point
.g_values[0].mean[m] = (unsigned char)first_frame->imageData[p + m];

}

//除第一以外的高斯分布函数的初始化(均值、权值和匹配次数都置零)

for( k = 1; k < params.n_gauss; k++)

{

bg_model->g_point
.g_values[k].weight = 0;

bg_model->g_point
.g_values[k].match_sum = 0;

for( m = 0; m < first_frame->nChannels; m++){

bg_model->g_point
.g_values[k].variance[m] = var_init;

bg_model->g_point
.g_values[k].mean[m] = 0;

}

}

}

} //g_point[]:像素，g_values[]:高斯分布函数，mean[]:通道

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;

}

cvUpdateBGStatModel(videoFrame,bgModel);

typedef int (CV_CDECL * CvUpdateBGStatModel)( IplImage* curr_frame, struct CvBGStatModel* bg_model );

/////////////////////////cvUpdateBGStatModel//////////////////////////////////

//函数功能：背景模型的更新，不仅要更新高斯分布函数的参数，还要更新各高斯函数的权重

static int CV_CDECL icvUpdateGaussianBGModel( IplImage* curr_frame, CvGaussBGModel* bg_model )

{

int i, j, k, n;

int region_count = 0;

CvSeq *first_seq = NULL, *prev_seq = NULL, *seq = NULL;

bg_model->countFrames++;

for( i = 0, n = 0; i < curr_frame->height; i++ )

{

for( j = 0; j < curr_frame->width; j++, n++ )

{

int match[CV_BGFG_MOG_MAX_NGAUSSIANS];
//对高斯函数做标记，match[m]=1表示函数m为匹配的高斯分布函数

double sort_key[CV_BGFG_MOG_MAX_NGAUSSIANS];
//此数组存贮每个高斯函数的均值与方差比值

const int nChannels = curr_frame->nChannels;

const int p = curr_frame->widthStep*i+j*nChannels;

CvGaussBGPoint* g_point = &bg_model->g_point
;

const CvGaussBGStatModelParams bg_model_params = bg_model->params;

double pixel[4]; // pixel[]存贮当前像素的各通道RGB值

int no_match;

for( k = 0; k < nChannels; k++ )

pixel[k] = (uchar)curr_frame->imageData[p+k];

no_match = icvMatchTest( pixel, nChannels, match, g_point, &bg_model_params );//检查是否有与当前像素匹配的高斯函数

if( bg_model->countFrames >= bg_model->params.win_size ) ？？？？？？？？？？？？？

{

icvUpdateFullWindow( pixel, nChannels, match, g_point, &bg_model->params );

if( no_match == -1)

icvUpdateFullNoMatch( curr_frame, p, match, g_point, &bg_model_params );

}

else

{

icvUpdatePartialWindow( pixel, nChannels, match, g_point, &bg_model_params );

if( no_match == -1)

icvUpdatePartialNoMatch( pixel, nChannels, match, g_point, &bg_model_params );

}

icvGetSortKey( nChannels, sort_key, g_point, &bg_model_params );

icvInsertionSortGaussians( g_point, sort_key, (CvGaussBGStatModelParams *)&bg_model_params );

icvBackgroundTest( nChannels, n, i, j, match, bg_model );

}

}

//foreground filtering

//filter small regions

cvClearMemStorage(bg_model->storage);

//cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_OPEN, 1 );

//cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_CLOSE, 1 );

cvFindContours( bg_model->foreground, bg_model->storage, &first_seq, sizeof(CvContour), CV_RETR_LIST );

for( seq = first_seq; seq; seq = seq->h_next )

{

CvContour* cnt = (CvContour*)seq;

if( cnt->rect.width * cnt->rect.height < bg_model->params.minArea )

{

//delete small contour

prev_seq = seq->h_prev;

if( prev_seq )

{

prev_seq->h_next = seq->h_next;

if( seq->h_next ) seq->h_next->h_prev = prev_seq;

}

else

{

first_seq = seq->h_next;

if( seq->h_next ) seq->h_next->h_prev = NULL;

}

}

else

{

region_count++;

}

}

bg_model->foreground_regions = first_seq;

cvZero(bg_model->foreground);

cvDrawContours(bg_model->foreground, first_seq, CV_RGB(0, 0, 255), CV_RGB(0, 0, 255), 10, -1);

return region_count;

}

/////////////////////////////////////icvMatchTest//////////////////////////////////////////

//函数功能：拿当前像素的值与已存在的高斯分布函数比较，查找是否存在匹配的的高斯分布函数，如果有则返回 k值(高斯分布函数的序号)

static int icvMatchTest( double* src_pixel, int nChannels, int* match,

const CvGaussBGPoint* g_point,

const CvGaussBGStatModelParams *bg_model_params )

{ //参数的传递：src_pixel为piexl[]:即当前像素的各通道值

int k;

int matchPosition=-1;

for ( k = 0; k < bg_model_params->n_gauss; k++)

match[k]=0;

for ( k = 0; k < bg_model_params->n_gauss; k++)

if (g_point->g_values[k].match_sum > 0)

{

double sum_d2 = 0.0;

double var_threshold = 0.0;

for(int m = 0; m < nChannels; m++)

{

double d = g_point->g_values[k].mean[m]- src_pixel[m];//通道m的原始模型值与当前像素的值之差

sum_d2 += (d*d);

var_threshold += g_point->g_values[k].variance[m];

}

//当前sum_d2为d0,d1,d2的平方和，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; //如果和第k个高斯函数匹配，则终止与后续函数的匹配

}

}

return matchPosition;

}

///////////////////////icvUpdateFullWindow////////////////////////////////////////

//函数功能：更新每个高斯分布的权值(对匹配的高斯函数k加大权值，其余的则减小权值)，如果前面的结果中存在匹配的高斯分布函数k,则需要再对第k个高斯分布函数的均值mean和方差variance做修正

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

CvGaussBGPoint* g_point,

const CvGaussBGStatModelParams *bg_model_params )

{ //参数的传递：src_pixel为piexl[]:即当前帧中该像素的RGB值

const double learning_rate_weight = (1.0/(double)bg_model_params->win_size);
//用户自定义模型学习率a

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->g_values[k].weight +

(learning_rate_weight*((double)match[k] -g_point->g_values[k].weight));

if(match[k]) //如果存在匹配的高斯分布函数k（当前像素为背景像素）,则需要再对第k个高斯分布函数的均值mean和方差variance更新

{

double learning_rate_gaussian = (double)match[k]/(g_point->g_values[k].weight*

(double)bg_model_params->win_size);//参数学习率p(p=a/w)

for(int m = 0; m < nChannels; m++)

{ //参数更新公式：u=(1-p)*u0+p*x; o*o=(1-p)*o*o+p*tmpDiff*tmpDiff

const double tmpDiff = src_pixel[m] - g_point->g_values[k].mean[m];
//当前像素的通道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]));

}

}

}

}

/////////////////////////icvUpdatePartialWindow/////////////////////////////

//函数功能：对所有的高斯分布函数做更新.至少每个高斯分布的权值必须修正，如果前面的结果中存在匹配的高斯分布函数k,则需要再对第k个高斯分布函数的match_sum修改，最终对那些匹配的高斯分布函数k的参数match_sum>0的做均值mean和方差variance修正

static void icvUpdatePartialWindow( double* src_pixel, int nChannels, int* match, CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params )

{

int k, m;

int window_current = 0;

for( k = 0; k < bg_model_params->n_gauss; k++ )

window_current += g_point->g_values[k].match_sum; //window_current为k个高斯分布函数的match_sum值之和

for( k = 0; k < bg_model_params->n_gauss; k++ )

{

g_point->g_values[k].match_sum += match[k]; //修正匹配的高斯分布函数k的match_sum值

double learning_rate_weight = (1.0/((double)window_current + 1.0));
//increased by one since sum

//修正每个高斯分布的权值

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( g_point->g_values[k].match_sum > 0 && match[k] )

{

double learning_rate_gaussian = (double)match[k]/((double)g_point->g_values[k].match_sum);

for( m = 0; m < nChannels; m++ )

{

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]));

}

}

}

}

//////////////////////////icvUpdateFullNoMatch//////////////////////////

//函数功能：当所有的高斯函数均不匹配时，说明有新的分布出现，需要将原高斯函数中sort_key最小的替换为新的高斯函数(权值小，方差大),其余的高斯函数对应的只需更新权值

static void icvUpdateFullNoMatch( IplImage* gm_image, int p, int* match,

CvGaussBGPoint* g_point,

const CvGaussBGStatModelParams *bg_model_params)

{ //参数的传递：gm_image为当前帧curr_frame

int k, m;

double alpha;

int match_sum_total = 0;

g_point->g_values[bg_model_params->n_gauss - 1].match_sum = 1;
//将新的高斯分布函数的match_sum置为1

for( k = 0; k < bg_model_params->n_gauss ; k++ )

match_sum_total += g_point->g_values[k].match_sum;

g_point->g_values[bg_model_params->n_gauss - 1].weight = 1./(double)match_sum_total;
//要给新的高斯分布函数赋一个较小的权值

//将新的高斯分布函数的variance[m]全部置为variance_init；mean[m]的值置为当前像素各通道的值

for( m = 0; m < gm_image->nChannels ; m++ )

{

g_point->g_values[bg_model_params->n_gauss - 1].variance[m] = bg_model_params->variance_init;

g_point->g_values[bg_model_params->n_gauss - 1].mean[m] = (unsigned char)gm_image->imageData[p + m];

}

//对其他的高斯分布函数做权值更新：w=(1-a)*w+a*m (a:模型学习率，m是匹配，匹配就是1，不匹配就是0)

alpha = 1.0 - (1.0/bg_model_params->win_size);
//alpha=1-a;

for( k = 0; k < bg_model_params->n_gauss - 1; k++ )

{

g_point->g_values[k].weight *= alpha;

if( match[k] )

g_point->g_values[k].weight += alpha;

}

}

////////////////////////////icvUpdatePartialNoMatch////////////////////////////////

static void

icvUpdatePartialNoMatch(double *pixel,

int nChannels,

int* /*match*/,

CvGaussBGPoint* g_point,

const CvGaussBGStatModelParams *bg_model_params)

{

int k, m;

//new value of last one

g_point->g_values[bg_model_params->n_gauss - 1].match_sum = 1;

//get sum of all but last value of match_sum

int match_sum_total = 0;

for(k = 0; k < bg_model_params->n_gauss ; k++)

match_sum_total += g_point->g_values[k].match_sum;

for(m = 0; m < nChannels; m++)

{

//first pass mean is image value

g_point->g_values[bg_model_params->n_gauss - 1].variance[m] = bg_model_params->variance_init;

g_point->g_values[bg_model_params->n_gauss - 1].mean[m] = pixel[m];

}

for(k = 0; k < bg_model_params->n_gauss; k++)

{ //更新所有高斯分布的权值

g_point->g_values[k].weight = (double)g_point->g_values[k].match_sum /

(double)match_sum_total;

}

}

/////////////////////////////////icvGetSortKey///////////////////////////////////

//函数功能：计算各个高斯分布weight/sqrt(variance_sum)的值，后面将对该值进行排序（该值越大则表示背景的可能性就越大）

static void icvGetSortKey( const int nChannels, double* sort_key, const CvGaussBGPoint* g_point,

const CvGaussBGStatModelParams *bg_model_params )

{

int k, m;

for( k = 0; k < bg_model_params->n_gauss; k++ )

{

// Avoid division by zero

if( g_point->g_values[k].match_sum > 0 )

{

// Independence assumption between components

double variance_sum = 0.0;

for( m = 0; m < nChannels; m++ )

variance_sum += g_point->g_values[k].variance[m];

sort_key[k] = g_point->g_values[k].weight/sqrt(variance_sum); //sort_key=w/(o*o)

}

else

sort_key[k]= 0.0;

}

}

//////////////////////////////icvInsertionSortGaussians////////////////////////////

static void icvInsertionSortGaussians( CvGaussBGPoint* g_point, double* sort_key, CvGaussBGStatModelParams *bg_model_params )

{

int i, j;

for( i = 1; i < bg_model_params->n_gauss; i++ )

{

double index = sort_key[i];

for( j = i; j > 0 && sort_key[j-1] < index; j-- )
//对sort_key[]按降序排序

{

double temp_sort_key = sort_key[j];

sort_key[j] = sort_key[j-1];

sort_key[j-1] = temp_sort_key;

CvGaussBGValues temp_gauss_values = g_point->g_values[j];

g_point->g_values[j] = g_point->g_values[j-1];

g_point->g_values[j-1] = temp_gauss_values;

}

// sort_key[j] = index;

}

}

///////////////////////////////////icvBackgroundTest/////////////////////////

static void icvBackgroundTest( const int nChannels, int n, int i, int j, int *match, CvGaussBGModel* bg_model )

{

int m, b;

uchar pixelValue = (uchar)255; // 像素默认都为前景

double weight_sum = 0.0;

CvGaussBGPoint* g_point = bg_model->g_point;

for( m = 0; m < nChannels; m++)？？？？？？？？？？？？？

bg_model->background->imageData[ bg_model->background->widthStep*i + j*nChannels + m] = (unsigned char)(g_point
.g_values[0].mean[m]+0.5);

for( b = 0; b < bg_model->params.n_gauss; b++)

{

weight_sum += g_point
.g_values[b].weight;

if( match[b] )

pixelValue = 0; //if为真，说明该像素已与某高斯函数匹配，该像素为背景

if( weight_sum > bg_model->params.bg_threshold )

break; //如果if语句为真，则前b个高斯分布被选为描述背景的函数

}

bg_model->foreground->imageData[ bg_model->foreground->widthStep*i + j] = pixelValue;

}

更新于2011年6月2日

上个星期又一次和师兄讨论了高斯背景的这块知识，此时距这篇博文的发布已经是一年之后，除了觉得之前对高斯的理解太过肤浅之外，真的是不得不感慨岁月蹉跎啊！在这一年当中，也让我对图像处理有了一些不同的认识：图像研究到最后，都划归了数学的范畴，所以大学里面学到的高等数学、概率、线性代数这三门数学基础学科都是必须掌握的。

听了师兄对高斯背景模型的讲述之后，让我意识到自己搞图像处理的一个弊病：图像处理中相关的理念并没有真正的领会，比如本文涉及到的核心概念高斯分布，在研究之前完全没有真正理解高斯和图像处理的背景建模为什么能对应起来，所以在看程序时就是一味的去按照论文的顺序去理解代码设计的初衷，我们已经完全局限在了这篇已有的论文里面，这样的思路怎么可能有理论上的突破和创新！

建议所有将要研究高斯背景模型算法的，一定要搞清楚高斯分布的特点，为什么自然中的背景可以用高斯分布来模拟？我们建立并更新高斯背景模型的理由又是什么？这些我觉得都应该是研究算法和代码之前必须搞明白的问题！

也欢迎对高斯背景模型有更深认识和了解的留言，一起交流！

本文转载自：http://hi.baidu.com/mailrabbit/blog/item/5225b0d27a7d6b38970a16cf.html/cmtid/46c05b17e3b4e74720a4e985