混合高斯背景建模——opencv

#include<opencv2/highgui.hpp>
#include<opencv2/core.hpp>
#include<opencv2/opencv.hpp>
using namespace std;
using namespace cv;
int main()
{
CvCapture*capture = cvCreateFileCapture("C:\\Users\\aoe\\Desktop\\avi\\walk.avi");//读取视频
IplImage*mframe = cvQueryFrame(capture);//读取视频中的一帧

int height = mframe->height;
int width = mframe->width;
int C = 4;//number of gaussian components
int M = 4;//number of background components
int std_init = 6;//initial standard deviation
double D = 2.5;
double T = 0.7;
double alpha = 0.01;
double p = alpha / (1 / C);
double thresh = 0.25;
int min_index = 0;
int*rank_ind=0;
int i, j,k,m;
int rand_temp=0;
int rank_ind_temp = 0;
CvRNG state;
IplImage*current = cvCreateImage(cvSize(mframe->width, mframe->height), IPL_DEPTH_8U, 1);
IplImage*test = cvCreateImage(cvSize(mframe->width, mframe->height), IPL_DEPTH_8U, 1);
IplImage*frg = cvCreateImage(cvSize(mframe->width, mframe->height), IPL_DEPTH_8U, 1);

double*mean = (double*)malloc(sizeof(double)*width*height*C);//pixelmeans
double*std = (double*)malloc(sizeof(double)*width*height*C);//pixel standard deviations
double*w = (double*)malloc(sizeof(double)*width*height*C);//权值
double*u_diff = (double*)malloc(sizeof(double)*width*height*C);//存放像素值与每一个单高斯模式的均值的差值
int*bg_bw = (int*)malloc(sizeof(int)*width*height);
double*rank = (double*)malloc(sizeof(double) * 1 * C);

//初始化
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-255之间的随机数
w[i*width*C + j*C + k] = (double)1 / C;//每一个单高斯模型的权值系数
std[i*width*C + j*C + k] = std_init;
}
}
}

while (1)
{
rank_ind = (int*)malloc(sizeof(int)*C);

cvCvtColor(mframe, current, CV_RGB2GRAY);//灰度化

//对于每一个像素，分别计算它和每一个单高斯模型的均值的差值
for (i = 0; i < height; i++)//对于每一个像素
{
for (j = 0; j < width; j++)
{
for (k = 0; k < C; k++)
{
u_diff[i*width*C + j*C + k] = abs((uchar)current->imageData[i*width + j] - mean[i*width*C + j*C + k]);

}
}
}

for (i = 0; i < height; i++)
{
for (j = 0; j < width; j++)
{
int match = 0;
double temp = 0;
double single_temp = 0;
//遍历所有的单高斯模式，如果此像素满足任一单高斯模式，则匹配；如果此像素不满足任何的单高斯模式，则不匹配
for (k = 0; k < C; k++)
{
if (abs(u_diff[i*width*C + j*C + k]) < D*std[i*width*C + j*C + k])//如果像素匹配某单个高斯模式，则对其权值、均值和标准差进行更新
{
match = 1;
w[i*width*C + j*C + k] += alpha*(1 - w[i*width*C + j*C + k]);//更新权值
p = alpha / 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];//更新均值
std[i*width*C + j*C + k] = sqrt((1 - p)*(std[i*width*C + j*C + k] * std[i*width*C + j*C + k]) + p*(pow((uchar)current->imageData[i*width + j] - mean[i*width*C + j*C + k], 2)));//更新标准差

}
else
{
w[i*width*C + j*C + k] = (1 - alpha)*w[i*width*C + j*C + k];//如果像素不符合某单个高斯模型，则将此单高斯模型的权值降低
}
}

if (match == 1)//如果和任一单高斯模式匹配，则将权值归一化
{
for (k = 0; k < C; k++)
{
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;//权值归一化，使得所有权值和为1
}
}
else//如果和所有单高斯模式都不匹配，则寻找权值最小的高斯模式并删除，然后增加一个新的高斯模式
{
single_temp = w[i*width*C + j*C];
for (k = 0; k < C; k++)
{
if (w[i*width*C + j*C + k] < single_temp)
{
min_index = k;//寻找权值最小的高斯模式
single_temp = w[i*width*C + j*C + k];
}

}
mean[i*width*C + j*C + min_index] = (uchar)current->imageData[i*width + j];//建立一个新的高斯模式，均值为当前像素值
std[i*width*C + j*C + min_index] = std_init;//标准差为初始值

for (k = 0; k < C; k++)
{
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;//权值归一化，使得所有权值和为1
}

}

for (k = 0; k < C; k++)//计算每个单高斯模式的重要性
{
rank[k] = w[i*width*C + j*C + k] / std[i*width*C + j*C + k];
rank_ind[k] = k;
}

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

bg_bw[i*width + j] = 0;
for (k = 0; k < C; k++)//如果前几个单高斯模式的重要性之和大于T，则将这前几个单高斯模式认为为背景模型
{
temp += w[i*width*C + j*C + rank_ind[k]];
bg_bw[i*width + j] += mean[i*width*C + j*C + rank_ind[k]] * w[i*width*C + j*C + rank_ind[k]];
if (temp >= T)
{
M = k;
break;
}
}

test->imageData[i*width + j] = (uchar)bg_bw[i*width + j];//背景图像

match = 0; k = 0;
while ((match == 0) && (k <= M))//如果某像素不符合背景模型中任一单高斯模型，则此像素为前景像素
{
if (abs(u_diff[i*width*C + j*C + rank_ind[k]]) <= D*std[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 += 1;

}

}
}
mframe = cvQueryFrame(capture);
if (mframe == NULL)
return -1;
cvNamedWindow("frg");
cvShowImage("frg", frg);
cvNamedWindow("back");
cvShowImage("back", test);
char s = cvWaitKey(33);
//if (s == 27)
//  break;
free(rank_ind);

}
cvWaitKey();
return 0;
}