高斯混合背景建模

1、高斯混合背景建模原理

混合高斯背景建模是基于像素样本统计信息的背景表示方法，利用像素在较长时间内大量样本值的概率密度等统计信息(如模式数量、每个模式的均值和标准差)表示背景，然后使用统计差分(如3σ原则)进行目标像素判断，可以对复杂动态背景进行建模，计算量较大。说白了就是在一个有限的时间系列，对一个位置的像素值进行统计，得到高斯模型，而后来一个新的像素值，判断是否符合3σ原则,如果是，表示背景，不是，符合前景。而后更新背景模型。在混合高斯背景模型及opencv实现中，有详细的介绍，博主这里给出整理。

在这里，我们给出假设：在混合高斯背景模型中，认为像素之间的颜色信息互不相关，对各像素点的处理都是相互独立的。对于视频图像中的每一个像素点，其值在序列图像中的变化可看作是不断产生像素值的随机过程，即用高斯分布来描述每个像素点的颜色呈现规律【单模态(单峰)，多模态(多峰)】。

对于多峰高斯分布模型，图像的每一个像素点按不同权值的多个高斯分布的叠加来建模，每种高斯分布对应一个可能产生像素点所呈现颜色的状态，各个高斯分布的权值和分布参数随时间更新。当处理彩色图像时，假定图像像素点R、G、B三色通道相互独立并具有相同的方差。对于随机变量X的观测数据集{x1,x2,…,xN}，xt=(rt,gt,bt)为t时刻像素的样本，则单个采样点xt其服从的混合高斯分布概率密度函数：

其中k为分布模式总数，η(xt,μi,t, τi,t)为t时刻第i个高斯分布，μi,t为其均值，τi,t为其协方差矩阵，δi,t为方差，I为三维单位矩阵，ωi,t为t时刻第i个高斯分布的权重。

详细算法流程为：

opencv实现：

#include "stdafx.h"  
#include "cv.h"  
#include "highgui.h"  

int _tmain(int argc, _TCHAR* argv[])
{
    CvCapture *capture = cvCreateFileCapture("F:\\毕业相关的程序\\测试数据集\\events_14-33-View_001.mov");
    IplImage *mframe, *current, *frg, *test;
    int *fg, *bg_bw, *rank_ind;
    double *w, *mean, *sd, *u_diff, *rank;
    int C, M, sd_init, i, j, k, m, rand_temp = 0, rank_ind_temp = 0, min_index = 0, x = 0, y = 0, counter_frame = 0;
    double D, alph, thresh, p, temp;
    CvRNG state;//基本随机数
    int match, height, width;
    mframe = cvQueryFrame(capture);

    frg = cvCreateImage(cvSize(mframe->width, mframe->height), IPL_DEPTH_8U, 1);
    current = cvCreateImage(cvSize(mframe->width, mframe->height), IPL_DEPTH_8U, 1);
    test = cvCreateImage(cvSize(mframe->width, mframe->height), IPL_DEPTH_8U, 1);

    C = 4;                      //number of gaussian components (typically 3-5)  
    M = 4;                      //number of background components  
    sd_init = 6;                //initial standard deviation (for new components) var = 36 in paper  初始化偏差
    alph = 0.01;                //learning rate (between 0 and 1) (from paper 0.01)  学习率
    D = 2.5;                    //positive deviation threshold                       积极偏差阈值
    thresh = 0.25;              //foreground threshold (0.25 or 0.75 in paper)       前景阈值  
    p = alph / (1 / C);         //initial p variable (used to update mean and sd)  

    height = current->height; width = current->widthStep;

    fg = (int *)malloc(sizeof(int)*width*height);                   //foreground array  前景数组
    bg_bw = (int *)malloc(sizeof(int)*width*height);                //background array  背景数组
    rank = (double *)malloc(sizeof(double) * 1 * C);                    //rank of components (w/sd) 级别的组件 
    w = (double *)malloc(sizeof(double)*width*height*C);            //weights array  权重数组
    mean = (double *)malloc(sizeof(double)*width*height*C);         //pixel means  均值
    sd = (double *)malloc(sizeof(double)*width*height*C);           //pixel standard deviations 像素标准偏差 
    u_diff = (double *)malloc(sizeof(double)*width*height*C);       //difference of each pixel from mean  每个像素与平均值的差异

    for (i = 0; i<height; i++)
    {
        for (j = 0; j<width; j++)
        {
            for (k = 0; k<C; k++)
            {
                mean[i*width*C + j*C + k] = cvRandReal(&state) * 255;//返回均匀分布，0-1之间的小数
                w[i*width*C + j*C + k] = (double)1 / C;//
                sd[i*width*C + j*C + k] = sd_init;//偏差=6
            }
        }
    }

    while (1){
        rank_ind = (int *)malloc(sizeof(int)*C);
        cvCvtColor(mframe, current, CV_BGR2GRAY);//rgb空间转换为GRAY空间
        // calculate difference of pixel values from mean  从均值中计算像素的差
        for (i = 0; i<height; i++)
        {
            for (j = 0; j<width; j++)
            {
                for (m = 0; m<C; m++)
                {
                    u_diff[i*width*C + j*C + m] = abs((uchar)current->imageData[i*width + j] - mean[i*width*C + j*C + m]);//计算当前像素值和均值的绝对值
                }
            }
        }
        //update gaussian components for each pixel  问每一个像素更新高斯参数  
        for (i = 0; i<height; i++)
        {
            for (j = 0; j<width; j++)
            {
                match = 0;//匹配值
                temp = 0;
                for (k = 0; k<C; k++)
                {
                    if (abs(u_diff[i*width*C + j*C + k]) <= D*sd[i*width*C + j*C + k])      //pixel matches component像素匹配部分 如果像素的偏差在2.5个方差之内，那么认为是背景 
                    {
                        match = 1;                                                 // variable to signal component match  

                        //update weights, mean, sd, p  更新权值，均值和像素标准偏差，
                        w[i*width*C + j*C + k] = (1 - alph)*w[i*width*C + j*C + k] + alph;//更新权重
                        p = alph / w[i*width*C + j*C + k];//更新
                        mean[i*width*C + j*C + k] = (1 - p)*mean[i*width*C + j*C + k] + p*(uchar)current->imageData[i*width + j];//更新均值
                        sd[i*width*C + j*C + k] = sqrt((1 - p)*(sd[i*width*C + j*C + k] * sd[i*width*C + j*C + k]) + p*(pow((uchar)current->imageData[i*width + j] - mean[i*width*C + j*C + k], 2)));//sqrt'非负实数的平方根函数
                    }
                    else{
                        w[i*width*C + j*C + k] = (1 - alph)*w[i*width*C + j*C + k];           // weight slighly decreases如果不属于，权重减少  
                    }
                    temp += w[i*width*C + j*C + k];//临时存储权重
                }

                for (k = 0; k<C; k++)
                {
                    w[i*width*C + j*C + k] = w[i*width*C + j*C + k] / temp;
                }

                temp = w[i*width*C + j*C];
                bg_bw[i*width + j] = 0;//背景数组
                for (k = 0; k<C; k++)
                {
                    bg_bw[i*width + j] = bg_bw[i*width + j] + mean[i*width*C + j*C + k] * w[i*width*C + j*C + k];//每一个像素的背景数组为当前的值+权重*均值
                    if (w[i*width*C + j*C + k] <= temp)//如果权重小于此时的背景数组
                    {
                        min_index = k;
                        temp = w[i*width*C + j*C + k];
                    }
                    rank_ind[k] = k;
                }

                test->imageData[i*width + j] = (uchar)bg_bw[i*width + j];//测数据中存入数值

                //if no components match, create new component  
                if (match == 0)
                {
                    mean[i*width*C + j*C + min_index] = (uchar)current->imageData[i*width + j];//均值为当前像素值
                    //printf("%d ",(uchar)bg->imageData[i*width+j]);  
                    sd[i*width*C + j*C + min_index] = sd_init;//标准差为为初始最大值
                    //权重为最小值没体现
                }
                for (k = 0; k<C; k++)
                {
                    rank[k] = w[i*width*C + j*C + k] / sd[i*width*C + j*C + k];//把权重除以方差从大到小排序
                    //printf("%f ",w[i*width*C+j*C+k]);  
                }

                //sort rank values  
                for (k = 1; k<C; k++)
                {
                    for (m = 0; m<k; m++)
                    {
                        if (rank[k] > rank[m])
                        {
                            //swap max values  
                            rand_temp = rank[m];
                            rank[m] = rank[k];
                            rank[k] = rand_temp;

                            //swap max index values  
                            rank_ind_temp = rank_ind[m];
                            rank_ind[m] = rank_ind[k];
                            rank_ind[k] = rank_ind_temp;
                        }
                    }
                }

                //calculate foreground  计算背景
                match = 0; k = 0;
                //frg->imageData[i*width+j]=0;  
                while ((match == 0) && (k<M)){
                    if (w[i*width*C + j*C + rank_ind[k]] >= thresh)
                        if (abs(u_diff[i*width*C + j*C + rank_ind[k]]) <= D*sd[i*width*C + j*C + rank_ind[k]]){
                            frg->imageData[i*width + j] = 0;
                            match = 1;
                        }
                        else
                            frg->imageData[i*width + j] = (uchar)current->imageData[i*width + j];
                    k = k + 1;
                }
            }
        }

        mframe = cvQueryFrame(capture);
        cvShowImage("fore", frg);
        cvShowImage("back", test);
        char s = cvWaitKey(33);
        if (s == 27) break;
        free(rank_ind);
    }

    free(fg); free(w); free(mean); free(sd); free(u_diff); free(rank);
    cvNamedWindow("back", 0);
    cvNamedWindow("fore", 0);
    cvReleaseCapture(&capture);
    cvDestroyWindow("fore");
    cvDestroyWindow("back");
    return 0;
}

代码说明：从单点像素复现了高斯混合背景建模，算法基础，选择的是固定的高斯函数个数。从运行结果看，对分辨率低的视频检测效果不好，对大面积移动区域有ghost区域，有噪声。在前景检测部分，可以做出改进。对检测的结果进行形态学处理，以达到滤波和消去空洞的目的。