混合高斯背景建模原理及实现

前些日子一直在忙答辩的事情,毕业后去了华为,图像处理什么的都派不上用场了。打算分3-4篇文章,把我研究生阶段学过的常用算法为大家和4107的师弟师妹们分享下。本次介绍混合高斯背景建模算法,还是老样子,首先介绍理论部分,然后给出代码,最后实验贴图。

 

一、理论

混合高斯背景建模是基于像素样本统计信息的背景表示方法,利用像素在较长时间内大量样本值的概率密度等统计信息(如模式数量、每个模式的均值和标准差)表示背景,然后使用统计差分(如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,tt时刻第i个高斯分布的权重。

 

详细算法流程:

 

二、代码实现

 

 
  1. // my_mixgaussians.cpp : 定义控制台应用程序的入口点。

  2. //

  3.  
  4. #include "stdafx.h"

  5. #include "cv.h"

  6. #include "highgui.h"

  7.  
  8. int _tmain(int argc, _TCHAR* argv[])

  9. {

  10. CvCapture *capture=cvCreateFileCapture("test.avi");

  11. IplImage *mframe,*current,*frg,*test;

  12. int *fg,*bg_bw,*rank_ind;

  13. double *w,*mean,*sd,*u_diff,*rank;

  14. 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;

  15. double D,alph,thresh,p,temp;

  16. CvRNG state;

  17. int match,height,width;

  18. mframe=cvQueryFrame(capture);

  19.  
  20. frg = cvCreateImage(cvSize(mframe->width,mframe->height),IPL_DEPTH_8U,1);

  21. current = cvCreateImage(cvSize(mframe->width,mframe->height),IPL_DEPTH_8U,1);

  22. test = cvCreateImage(cvSize(mframe->width,mframe->height),IPL_DEPTH_8U,1);

  23.  
  24. C = 4; //number of gaussian components (typically 3-5)

  25. M = 4; //number of background components

  26. sd_init = 6; //initial standard deviation (for new components) var = 36 in paper

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

  28. D = 2.5; //positive deviation threshold

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

  30. p = alph/(1/C); //initial p variable (used to update mean and sd)

  31.  
  32. height=current->height;width=current->widthStep;

  33.  
  34. fg = (int *)malloc(sizeof(int)*width*height); //foreground array

  35. bg_bw = (int *)malloc(sizeof(int)*width*height); //background array

  36. rank = (double *)malloc(sizeof(double)*1*C); //rank of components (w/sd)

  37. w = (double *)malloc(sizeof(double)*width*height*C); //weights array

  38. mean = (double *)malloc(sizeof(double)*width*height*C); //pixel means

  39. sd = (double *)malloc(sizeof(double)*width*height*C); //pixel standard deviations

  40. u_diff = (double *)malloc(sizeof(double)*width*height*C); //difference of each pixel from mean

  41.  
  42. for (i=0;i<height;i++)

  43. {

  44. for (j=0;j<width;j++)

  45. {

  46. for(k=0;k<C;k++)

  47. {

  48. mean[i*width*C+j*C+k] = cvRandReal(&state)*255;

  49. w[i*width*C+j*C+k] = (double)1/C;

  50. sd[i*width*C+j*C+k] = sd_init;

  51. }

  52. }

  53. }

  54.  
  55. while(1){

  56. rank_ind = (int *)malloc(sizeof(int)*C);

  57. cvCvtColor(mframe,current,CV_BGR2GRAY);

  58. // calculate difference of pixel values from mean

  59. for (i=0;i<height;i++)

  60. {

  61. for (j=0;j<width;j++)

  62. {

  63. for (m=0;m<C;m++)

  64. {

  65. u_diff[i*width*C+j*C+m] = abs((uchar)current->imageData[i*width+j]-mean[i*width*C+j*C+m]);

  66. }

  67. }

  68. }

  69. //update gaussian components for each pixel

  70. for (i=0;i<height;i++)

  71. {

  72. for (j=0;j<width;j++)

  73. {

  74. match = 0;

  75. temp = 0;

  76. for(k=0;k<C;k++)

  77. {

  78. if (abs(u_diff[i*width*C+j*C+k]) <= D*sd[i*width*C+j*C+k]) //pixel matches component

  79. {

  80. match = 1; // variable to signal component match

  81.  
  82. //update weights, mean, sd, p

  83. w[i*width*C+j*C+k] = (1-alph)*w[i*width*C+j*C+k] + alph;

  84. p = alph/w[i*width*C+j*C+k];

  85. mean[i*width*C+j*C+k] = (1-p)*mean[i*width*C+j*C+k] + p*(uchar)current->imageData[i*width+j];

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

  87. }else{

  88. w[i*width*C+j*C+k] = (1-alph)*w[i*width*C+j*C+k]; // weight slighly decreases

  89. }

  90. temp += w[i*width*C+j*C+k];

  91. }

  92.  
  93. for(k=0;k<C;k++)

  94. {

  95. w[i*width*C+j*C+k] = w[i*width*C+j*C+k]/temp;

  96. }

  97.  
  98. temp = w[i*width*C+j*C];

  99. bg_bw[i*width+j] = 0;

  100. for (k=0;k<C;k++)

  101. {

  102. bg_bw[i*width+j] = bg_bw[i*width+j] + mean[i*width*C+j*C+k]*w[i*width*C+j*C+k];

  103. if (w[i*width*C+j*C+k]<=temp)

  104. {

  105. min_index = k;

  106. temp = w[i*width*C+j*C+k];

  107. }

  108. rank_ind[k] = k;

  109. }

  110.  
  111. test->imageData[i*width+j] = (uchar)bg_bw[i*width+j];

  112.  
  113. //if no components match, create new component

  114. if (match == 0)

  115. {

  116. mean[i*width*C+j*C+min_index] = (uchar)current->imageData[i*width+j];

  117. //printf("%d ",(uchar)bg->imageData[i*width+j]);

  118. sd[i*width*C+j*C+min_index] = sd_init;

  119. }

  120. for (k=0;k<C;k++)

  121. {

  122. rank[k] = w[i*width*C+j*C+k]/sd[i*width*C+j*C+k];

  123. //printf("%f ",w[i*width*C+j*C+k]);

  124. }

  125.  
  126. //sort rank values

  127. for (k=1;k<C;k++)

  128. {

  129. for (m=0;m<k;m++)

  130. {

  131. if (rank[k] > rank[m])

  132. {

  133. //swap max values

  134. rand_temp = rank[m];

  135. rank[m] = rank[k];

  136. rank[k] = rand_temp;

  137.  
  138. //swap max index values

  139. rank_ind_temp = rank_ind[m];

  140. rank_ind[m] = rank_ind[k];

  141. rank_ind[k] = rank_ind_temp;

  142. }

  143. }

  144. }

  145.  
  146. //calculate foreground

  147. match = 0;k = 0;

  148. //frg->imageData[i*width+j]=0;

  149. while ((match == 0)&&(k<M)){

  150. if (w[i*width*C+j*C+rank_ind[k]] >= thresh)

  151. if (abs(u_diff[i*width*C+j*C+rank_ind[k]]) <= D*sd[i*width*C+j*C+rank_ind[k]]){

  152. frg->imageData[i*width+j] = 0;

  153. match = 1;

  154. }

  155. else

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

  157. k = k+1;

  158. }

  159. }

  160. }

  161.  
  162. mframe = cvQueryFrame(capture);

  163. cvShowImage("fore",frg);

  164. cvShowImage("back",test);

  165. char s=cvWaitKey(33);

  166. if(s==27) break;

  167. free(rank_ind);

  168. }

  169.  
  170. free(fg);free(w);free(mean);free(sd);free(u_diff);free(rank);

  171. cvNamedWindow("back",0);

  172. cvNamedWindow("fore",0);

  173. cvReleaseCapture(&capture);

  174. cvDestroyWindow("fore");

  175. cvDestroyWindow("back");

  176. return 0;

  177. }


实验结果:

 

前景:

背景:

  • 1
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值