混合高斯背景建模

一、理论

混合高斯背景建模是基于像素样本统计信息的背景表示方法，利用像素在较长时间内大量样本值的概率密度等统计信息(如模式数量、每个模式的均值和标准差)表示背景，然后使用统计差分(如3σ原则)进行目标像素判断，可以对复杂动态背景进行建模，计算量较大。

在混合高斯背景模型中，认为像素之间的颜色信息互不相关，对各像素点的处理都是相互独立的。对于视频图像中的每一个像素点，其值在序列图像中的变化可看作是不断产生像素值的随机过程，即用高斯分布来描述每个像素点的颜色呈现规律【单模态(单峰)，多模态(多峰)】。
对于多峰高斯分布模型，图像的每一个像素点按不同权值的多个高斯分布的叠加来建模，每种高斯分布对应一个可能产生像素点所呈现颜色的状态，各个高斯分布的权值和分布参数随时间更新。当处理彩色图像时，假定图像像素点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个高斯分布的权重。

详细算法流程：





由于能力有限，代码来自实验室的师兄写的代码，我从中学到了很多内容。

//混合高斯背景建模

/*****************************************************************************************************

* File name: MixtureGaussian.cpp

* Author: source code by RamiYim

* Version: cs-hhu-cn0.0

* Date: 20150208

* Description: background modelling by mixture gaussian

* Others:

* Function List: mixtureGaussian();

* History: <rami>,<20150208>,<version>,<desc>

*

******************************************************************************************************/

/*

#include "cv.h"

#include "highgui.h"

#include <iostream>

using namespace std;

/*************************************************************************************************

Function:void mixtureGaussian()

Description:

Calls:

Called By:main()

Input:

Output:

Return:

Others:建立彩色的背景

Histroy:<author><time><version><desc>

<RamiYim><20150505><0.0><第一次创建函数>

***************************************************************************************************/

/*

void mixtureGaussianColor();

void mixtureGaussianGrey();

int main()

{

int i =0;

mixtureGaussianColor();

//mixtureGaussianGrey();

std::cout << i << endl;

return 0;

}

void mixtureGaussianColor()

{

//读入视频

CvCapture *capture = cvCreateFileCapture("E:\\新建文件夹\\cs.AVI");

IplImage *sourceFrame = cvQueryFrame(capture);

//创建背景帧和前景帧，并分配内存

IplImage *backgroundFrame = cvCreateImage(cvSize(sourceFrame->width, sourceFrame->height), IPL_DEPTH_8U, 3);

IplImage *foregroundFrame = cvCreateImage(cvSize(sourceFrame->width, sourceFrame->height), IPL_DEPTH_8U, 3);

/*----------------初始化各种变量--------------*/

/*

int liv_componets = 4;
//高斯分布个数number of gaussian components (typically 3-5)

int liv_backgroundComponents = 4;
//背景分布个数number of background components

int liv_stdDeviation = 6;
//初始标准差initial standard deviation (for new components) var = 36 in paper

double lfv_alph = 0.01;
//学习速率a learning rate (between 0 and 1) (from paper 0.01)

double lfv_lamda = 2.5;
//偏差阈值λ positive deviation threshold

double lfv_thresh = 0.25;
//前景阈值 foreground threshold (0.25 or 0.75 in paper)

//初始rho，用于更新均值和标准差

double lfv_rho = lfv_alph / (1 / liv_componets);
//initial rho variable (used to update lpv_mean and lpv_std)

//初始化高度宽度，步长

int liv_height = sourceFrame->height;

int liv_width = sourceFrame->width;

int liv_widthStep = sourceFrame->widthStep;

/*-------------------申请空间，分配内存-----------*/

/*

//foreground array 前景数组

uchar *lpv_foregroundImageData = (uchar *) malloc(sizeof(char) * liv_height * liv_widthStep);

//background array 背景数组

uchar *lpv_backgroundImageData = (uchar *) malloc(sizeof(char) * liv_height * liv_widthStep);

//lpv_weights array 权值数组

double *lpv_weights = (double *) malloc(sizeof(double) * liv_height * liv_widthStep * liv_componets);

//pixel means 像素均值

double *lpv_mean = (double *) malloc(sizeof(double) * liv_height * liv_widthStep * liv_componets);

//pixel standard deviations 像素标准差

double *lpv_std = (double *) malloc(sizeof(double) * liv_height * liv_widthStep * liv_componets);

//difference of each pixel from lpv_mean

double *lpv_uDiff = (double *) malloc(sizeof(double) * liv_height * liv_widthStep * liv_componets);

//lpv_rank of components (weights / lpv_std)

//double *lpv_rank = (double *) malloc(sizeof(double) * liv_componets);

int *lpv_rankIndex = (int *) malloc(sizeof(int) * liv_componets);

// 为均值、标准差、权重赋初值

//opencv中c的rng类可以压缩一个64位的i整数并可以得到scalar和array的随机数。

//目前的版本支持均匀分布随机数和gaussian分布随机数。

CvRNG state; //初始化随机数生成器状态

for (int i = 0; i < liv_height; i++)

{

for (int j = 0; j < liv_widthStep; j++)

{

for(int k = 0; k < liv_componets; k++)

{

// 每一个高斯分量的下标

int liv_eachCompoent = i * liv_widthStep * liv_componets + j * liv_componets + k;

// 每一个通道的下标

int liv_eachChannels = i * liv_widthStep + j;

// 初始化均值、标准差、权重

lpv_mean[liv_eachCompoent] = cvRandReal(&state) * 255;

lpv_std[liv_eachCompoent] = liv_stdDeviation;

lpv_weights[liv_eachCompoent] = (double)1 / liv_componets;

}

}

}

while(sourceFrame = cvQueryFrame(capture))

{

for (int i = 0; i < liv_height; i++)

{

for (int j = 0; j < liv_widthStep; j++)

{

// 用于判断是否有高斯分量相匹配

int liv_match = 0;

float lfv_weightTotal = 0;

// 每一个通道的下标

int liv_eachChannel = i * liv_widthStep + j;

for (int k = 0; k < liv_componets; k++)

{

// 每一个高斯分量的下标

int liv_eachCompoent = i * liv_widthStep * liv_componets + j * liv_componets + k;

// 每一个高斯分量与其相对应通道的差值数组 |I-u|

lpv_uDiff[liv_eachCompoent] = abs((uchar)sourceFrame->imageData[liv_eachChannel]- lpv_mean[liv_eachCompoent]);

//匹配的

// |I-u| < lamda * std 时认为该像素是背景，需更新

if (lpv_uDiff[liv_eachCompoent] < lfv_lamda * lpv_std[liv_eachCompoent])

{

liv_match = 1; //liv_match变成了1

// 更新高斯分量权重, w = (1 - alph) * weights + alph * M (M = 1)

lpv_weights[liv_eachCompoent] = (1 - lfv_alph) * lpv_weights[liv_eachCompoent] + lfv_alph;

lfv_rho = lfv_alph / lpv_weights[liv_componets];

//更新高斯分量期望 u = (1 - rho) * u + rho * I

lpv_mean[liv_eachCompoent] = (1 - lfv_rho) * lpv_mean[liv_eachCompoent] + lfv_rho * (uchar)sourceFrame->imageData[liv_eachChannel];

//更新高斯分量标准差std = sqrt（(1 - rho) * std^2 + rho * (I - u)^2）

lpv_std[liv_eachCompoent] = sqrt((1 - lfv_rho) * pow(lpv_std[liv_eachCompoent], 2)

+ lfv_rho * pow((uchar)sourceFrame->imageData[liv_eachChannel] - lpv_mean[liv_eachCompoent], 2));

}

//未匹配的 M=0,这个时候均值和方差不变

else

{

// 更新高斯分量权重, w = (1 - alph) * weights + alph * M(M = 0)

// weight slighly decreases

lpv_weights[liv_eachCompoent] = (1 - lfv_alph) * lpv_weights[liv_eachCompoent];

}

//权重总和

lfv_weightTotal += lpv_weights[liv_eachCompoent];

} // end for (int k = 0; k < liv_componets; k++)

// 归一化权值

for (int k = 0; k < liv_componets; k++)

{

// 每一个高斯分量的下标

int liv_eachCompoent = i * liv_widthStep * liv_componets + j * liv_componets + k;

//归一化权值，权值变为权值除以总的权值

lpv_weights[liv_eachCompoent] = lpv_weights[liv_eachCompoent] / lfv_weightTotal;

}

// 背景提取

lpv_backgroundImageData[liv_eachChannel] = 0;

for (int k = 0; k < liv_componets; k++)

{

// 每一个高斯分量的下标

int liv_eachCompoent = i * liv_widthStep * liv_componets + j * liv_componets + k;

//背景数组=归一化的权值*均值

lpv_backgroundImageData[liv_eachChannel] += lpv_weights[liv_eachCompoent] * lpv_mean[liv_eachCompoent];

}

//处理过后的背景数组里面的数据指回背景帧内存空间

backgroundFrame->imageData[liv_eachChannel] = lpv_backgroundImageData[liv_eachChannel];

//没有一个高斯分量匹配的

// 如果没有高斯分量匹配，则将权值最小的替换，即该高斯分量的均值变为当前像素值，标准差为初始较大值，权重为较小值

if (liv_match == 0)

{

// 最小权值下标

int liv_minWeightsIndex = 0;

// 权值最小的值,初始值为第一个高斯分量值

float lfv_minWeights = lpv_weights[i * liv_widthStep * liv_componets + j * liv_componets + liv_minWeightsIndex];

//找出最小的权值lfv_minWeights，以及最小权值的下标liv_minWeightsIndex

for (int k = 0; k < liv_componets; k++)

{

int liv_eachCompoent = i * liv_widthStep * liv_componets + j * liv_componets + k;

if (lpv_weights[liv_eachCompoent]< lfv_minWeights)

{

lfv_minWeights = lpv_weights[liv_eachCompoent];

liv_minWeightsIndex = k;

}

}

//当前最小权值的高斯分量的均值变为当前像素值

lpv_mean[i * liv_widthStep * liv_componets + j * liv_componets + liv_minWeightsIndex] = (uchar)sourceFrame->imageData[liv_eachChannel];

//最小权值的高斯分量的标准差是初始的标准差值

lpv_std[i * liv_widthStep * liv_componets + j * liv_componets + liv_minWeightsIndex] = liv_stdDeviation;

}

} // end (int j = 0; j < liv_widthStep; j++)

} // end for (int i = 0; i < liv_height; i++)

cvShowImage("source", sourceFrame);

cvShowImage("background", backgroundFrame);

char s = cvWaitKey(1);

if(s == 27) break;

} // end while(sourceFrame = cvQueryFrame(capture))

}

/*************************************************************************************************

Function:void mixtureGaussianGrey()

Description:

Calls:

Called By:main()

Input:

Output:

Return:

Others:从金圣韬博客copy过来，建立灰度背景

Histroy:<author><time><version><desc>

<RamiYim><20150208><0.0><第一次创建函数>

***************************************************************************************************/

/*

void mixtureGaussianGrey()

{

CvCapture *capture = cvCreateFileCapture("D:\\VideoData\\SDC11945.AVI");

IplImage *lpv_sourceFrame = cvQueryFrame(capture);

//当前，单通道

IplImage *current = cvCreateImage(cvSize(lpv_sourceFrame->width, lpv_sourceFrame->height), IPL_DEPTH_8U, 1);

//前景帧，单通道

IplImage *lpv_foregroundFrame = cvCreateImage(cvSize(lpv_sourceFrame->width, lpv_sourceFrame->height), IPL_DEPTH_8U, 1);

//背景帧，单通道

IplImage *lpv_backgroundFrame = cvCreateImage(cvSize(lpv_sourceFrame->width, lpv_sourceFrame->height), IPL_DEPTH_8U, 1);

//设置参数

int liv_componets = 4;
//number of gaussian components (typically 3-5)

int liv_backgroundComponents = 4;
//number of background components

int liv_stdDeviation = 6;
//initial standard deviation (for new components) var = 36 in paper

double alph = 0.01;
//learning rate (between 0 and 1) (from paper 0.01)

double lamda = 2.5;
//positive deviation threshold

double thresh = 0.25;
//foreground threshold (0.25 or 0.75 in paper)

double rho = alph / (1 / liv_componets);
//initial rho variable (used to update lpv_mean and lpv_std)

int width = current->widthStep;

int height = current->height;

//申请空间，分配内存

int *lpv_foregroundImageData = (int *) malloc(sizeof(int) * width * height);
//foreground array

int *lpv_backgroundImageData = (int *) malloc(sizeof(int) * width * height);
//background array

double *lpv_weights = (double *) malloc(sizeof(double) * width * height * liv_componets);//lpv_weights array

double *lpv_mean = (double *) malloc(sizeof(double) * width * height * liv_componets);
//pixel means

double *lpv_std = (double *) malloc(sizeof(double) * width * height * liv_componets);
//pixel standard deviations

double *u_diff = (double *) malloc(sizeof(double) * width * height * liv_componets);
//difference of each pixel from lpv_mean

////lpv_rank of components (weights / lpv_std)

//分布的排序，按照weights / lpv_std排序

double *lpv_rank = (double *) malloc(sizeof(double) * 1 * liv_componets);

//排序的序列号

int *rank_index = (int *) malloc(sizeof(int) * liv_componets);

CvRNG state;

for (int i = 0; i < height; i++)

{

for (int j = 0; j < width; j++)

{

for(int k = 0; k < liv_componets; k++)

{

lpv_mean[i * width * liv_componets + j * liv_componets + k] = cvRandReal(&state) * 255;

lpv_weights[i * width * liv_componets + j * liv_componets + k] = (double)1 / liv_componets;

lpv_std[i * width * liv_componets + j * liv_componets + k] = liv_stdDeviation;

}

}

}

while(lpv_sourceFrame = cvQueryFrame(capture))

{

//转灰度图像，当前帧图像lpv_sourceFrame，变成灰度图像current

cvCvtColor(lpv_sourceFrame, current, CV_BGR2GRAY);

//update gaussian components for each pixel

for (int i = 0; i < height; i++)

{

for (int j = 0; j < width; j++)

{

int match = 0;

double temp = 0; //权重总和初始值0

for(int k = 0; k < liv_componets; k++)

{

//calculate difference of pixel values from lpv_mean

u_diff[i * width * liv_componets + j * liv_componets + k] = abs((uchar)current->imageData[i * width + j]

- lpv_mean[i * width * liv_componets + j * liv_componets + k]);

if (abs(u_diff[i * width * liv_componets + j * liv_componets + k]) <= lamda * lpv_std[i * width * liv_componets + j * liv_componets + k])//pixel matches component

{

match = 1;
// variable to single component match

//update lpv_weights, lpv_mean, lpv_std, rho

lpv_weights[i * width * liv_componets + j * liv_componets + k] = (1 - alph) * lpv_weights[i * width * liv_componets + j * liv_componets + k] + alph;

rho = alph / lpv_weights[i * width * liv_componets + j * liv_componets + k];

lpv_mean[i * width * liv_componets + j * liv_componets + k] = (1 - rho) * lpv_mean[i * width * liv_componets + j * liv_componets + k]

+ rho * (uchar)current->imageData[i * width + j];

lpv_std[i * width * liv_componets + j * liv_componets + k] = sqrt((1 - rho) * (lpv_std[i * width * liv_componets + j * liv_componets + k]

* lpv_std[i * width * liv_componets + j * liv_componets + k]) + rho * (pow((uchar)current->imageData[i * width + j]

- lpv_mean[i * width * liv_componets + j * liv_componets + k], 2)));

}

else

{

lpv_weights[i * width * liv_componets + j * liv_componets + k] = (1 - alph)

* lpv_weights[i * width * liv_componets+ j * liv_componets + k];// weight slighly decreases

}

//计算出权重的总和

temp += lpv_weights[i * width * liv_componets + j * liv_componets + k];

} //end for(int k = 0; k < liv_componets; k++)

//归一化权值

for(int k = 0; k < liv_componets; k++)

{

lpv_weights[i * width * liv_componets + j * liv_componets + k] = lpv_weights[i * width * liv_componets + j * liv_componets + k] / temp;

}

//temp值为第一个高斯分布的归一化之后的权值总和

temp = lpv_weights[i * width * liv_componets + j * liv_componets];

lpv_backgroundImageData[i * width + j] = 0;

int min_index = 0;

for (int k = 0; k < liv_componets; k++)

{

//背景数据=归一化后的权值*均值

lpv_backgroundImageData[i * width + j] = lpv_backgroundImageData[i * width + j]

+ lpv_weights[i * width * liv_componets + j * liv_componets + k] * lpv_mean[i * width * liv_componets + j * liv_componets + k];

//找出权值最小的高斯分布还有它的下标值

if (lpv_weights[i * width * liv_componets + j * liv_componets + k] <= temp)

{

min_index = k;

temp = lpv_weights[i * width * liv_componets + j * liv_componets + k];

}

}

//背景数据指回它对应的数据块

lpv_backgroundFrame->imageData[i * width + j] = (uchar)lpv_backgroundImageData[i * width + j];

//如果没有任何模式匹配的情况

//if no components match, create new component

if (match == 0)

{

lpv_mean[i * width * liv_componets + j * liv_componets + min_index] = (uchar)current->imageData[i * width + j];

lpv_std[i * width * liv_componets + j * liv_componets + min_index] = liv_stdDeviation;

}

for (int k = 0; k < liv_componets; k++)

{

rank_index[k] = k;

lpv_rank[k] = lpv_weights[i * width * liv_componets + j * liv_componets + k] / lpv_std[i * width * liv_componets + j * liv_componets + k];

}

//sort lpv_rank values

int rand_temp = 0, rank_index_temp = 0;

//冒泡排序

//各模式按照降序排列

for (int k = 1; k < liv_componets; k++)

{

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

{

if (lpv_rank[k] > lpv_rank[m])

{

// swap max values

rand_temp = lpv_rank[m];

lpv_rank[m] = lpv_rank[k];

lpv_rank[k] = rand_temp;

// swap max index values

rank_index_temp = rank_index[m];

rank_index[m] = rank_index[k];

rank_index[k] = rank_index_temp;

}

}

}

// calculate foreground

match = 0; int k = 0;

// lpv_foregroundFrame->imageData[i * width + j] = 0;

while ((match == 0) && (k < liv_backgroundComponents))

{

//选择前B个作为背景，B=arg(min(权值和>阈值））

if (lpv_weights[i * width * liv_componets + j * liv_componets + rank_index[k]] >= thresh)

{

if (abs(u_diff[i * width * liv_componets + j * liv_componets + rank_index[k]]) <=

lamda * lpv_std[i * width * liv_componets + j * liv_componets + rank_index[k]])

{

lpv_foregroundFrame->imageData[i * width + j] = 0; //背景，即前景中黑色部分，

match = 1;

}

else

{

lpv_foregroundFrame->imageData[i * width + j] = (uchar)current->imageData[i * width + j];

}

}

k = k + 1;

}

} // end for (int j = 0; j < width; j++)

} // end for (int i = 0; i < height; i++)

cvShowImage("source", lpv_sourceFrame);

cvShowImage("fore",lpv_foregroundFrame);

cvShowImage("back",lpv_backgroundFrame);

char s = cvWaitKey(1);

if(s == 27) break;

} // end while

free(rank_index);

free(lpv_foregroundImageData);

free(lpv_weights);

free(lpv_mean);

free(lpv_std);

free(u_diff);

free(lpv_rank);

cvReleaseCapture(&capture);

cvDestroyWindow("fore");

cvDestroyWindow("back");

}

*/