混合高斯背景建模
在运动目标检测提取中,背景目标对于目标的识别和跟踪至关重要。而建模正是背景目标提取的一个重要环节。
前景是指在假设背景为静止的情况下,任何有意义的运动物体即为前景。
运动物体检测的问题主要分为两类,摄像机固定和摄像机运动。对于摄像机运动的运动物体检测问题,比较著名的解决方案是光流法,通过求解偏微分方程求的图像序列的光流场,从而预测摄像机的运动状态。对于摄像机固定的情形,当然也可以用光流法,但是由于光流法的复杂性,往往难以实时的计算。而混合高斯背景建模适合于在摄像机固定的情况下从图像序列中分离出背景和前景。在摄像机固定的情况下,背景的变化是缓慢的,而且大都是光照,风等等的影响,通过对背景建模,对一幅给定图像分离前景和背景,一般来说,前景就是运动物体,从而达到运动物体检测的目的。
混合高斯模型已经广泛应用于鲁棒的复杂场景背景建模,特别是在有微小重复运动的场合,如摇动的树叶、灌木丛、旋转的风扇、海面波涛、雨雪天气、光线反射等。基于像素的混合高斯模型对多峰分布背景进行建模很有效,能适应背景的变化如光线渐变,并能基本满足实际应用中对算法的实时性要求。
算法基本思想
用GMM对背景建模的基本思想是把每一个像素点所呈现的颜色用K个高斯分布的叠加来表示,通常K取3-5之间。将像素点所呈现的颜色X认为是随机变量,则在每时刻t=1,…,T所得到视频帧图像的像素值只是随机变量X的采样值(观值)。
在进行前景检测前,先对背景进行训练,对每一帧图像中每个背景采用一个混合高斯模型进行模拟,背景一旦提取出来,前景的检测就简单了,检查像素是否与背景的高斯模型匹配,匹配是背景,不匹配就是前景。
所以关键就是混合高斯背景模型的建立。
GMM之所以能够将前景和背景分开是基于如下两点事实的:
(1)在长期观测的场景中,背景占大多数时间,更多的数据是支持背景分布的
(2)即使是相对颜色一致的运动物体也会比背景产生更多变化,况且一般情况下物体都是带有不同颜色的
算法大致这样:
首先初始化预先定义的几个高斯模型,对高斯模型中的参数进行初始化,并求出之后将要用到的参数。其次,对于每一帧中的每一个像素进行处理,看其是否匹配某个模型,若匹配,则将其归入该模型中,并对该模型根据新的像素值进行新,若不匹配,则以该像素建立一个高斯模型,初始化参数,代替原有模型中最不可能的模型。最后选择前面几个最有可能的模型作为背景模型,为背景目标提取做铺垫。
算法原理
背景图像的每一个像素分别用由K个高斯分布构成的混合高斯模型来建模, 即:
Xj表示像素j在t时刻的取值,如果像素为RGB三通道,则xj为向量,xj=[xjR xjG xjB],表示时刻t混合高斯模型中第i个高斯分布的权系数的估计值,和分别表示时刻t混合高斯模型中第i个高斯分布的均值向量和协方差矩阵(此处假定像素的红、绿、蓝分量相互独立);η表示高斯分布概率密度函数。
参数初始化
在第一帧图像时每个像素对应的第一个高斯分布进行初始化,均值赋为当前像素的值,权值赋为1,除第一以外的高斯分布函数的均值、权值和都初始化零。
参数更新
在时刻t对图像帧的每个像素Xt与它对应的高斯模型进行匹配,
匹配规则为:
如果像素值Xt 与混合高斯模型中第i个高斯分布Gi均值的距离小于其标准差的2.5倍, 则定义该高斯分布Gi 与像素值Xt 匹配。
如果检验出该像素混合高斯模型中至少有一个高斯分布与像素值Xt 匹配, 那么混合高斯模型的参数更新规则为:
1)对于不匹配的高斯分布, 它们的均值μ和协方差矩阵保持不变;
2)匹配的高斯分布Gi 的均值μ和协方差矩阵按下式更新:
α为参数估计的学习速率。
如果该像素对应的混合高斯模型中没有高斯分布与像素值Xt 匹配, 那么将最不可能代表背景过程的高斯分布Gj重新赋值, 即
式中, W0 和V0是预先给定的正值;I为一个3 ×3单位矩阵。
然后按下式更新所有K 个高斯分布在时刻t的权系数ωi,t:
式中, 如果高斯分布Gi与t时刻像素值Xt匹配, 则Mi,t取值1,否则取值为0。后一式完成权重的归一化,保证权重和为1。
背景模型估计及运动分割
首先按ωk / |k|的比值由大到小将组成每个像素混合高斯模型的K 个高斯分布排序, ωk代表了第k个分布产生的数据所占的比例,k代表第k个分布的方差,ωk / |k|的值越大,第k个分布是背景模型的可能性就越大。一个分布的权值增大或方差减小,都可以使ωk / |k|的值增大。
那么最有可能描述稳定背景过程的高斯分布将位于序列的顶部, 而由背景暂态扰动产生的分布将滑向序列的底部, 最终被新赋值的高斯分布所取代。
这样, 我们选择上述序列中前B 个高斯分布作为背景像素模型:
式中, T为预定的阈值(0. 5≤ T ≤ 1);B的含义是, 排序后的K个高斯分布中的前B 个高斯分布是该背景像素的最佳描述。
前景(运动)分割
现在, 从新检验t时刻每一个像素值Xt与其得到的前B 个高斯分布的匹配关系, 如果该像素值Xt与前B个高斯分布之一匹配, 则该像素为背景点, 否则该像素被分类为前景, 即运动物体。这样, 就完成了将运动物体从背景中分割出来的任务。
////////////////////////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[n].g_values = bg_model->g_point[0].g_values + n*params.n_gauss; bg_model->g_point[n].g_values[0].weight = 1; //权值赋为1 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; } } } } //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[n]; 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* , 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[n].g_values[0].mean[m]+0.5); for( b = 0; b < bg_model->params.n_gauss; b++) { weight_sum += g_point[n].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; }
下面是我的检测结果,分别为原视频,前景和背景
版权声明: