tornadomeet
前景檢測算法_4(opencv自帶GMM)
http://www.cnblogs.com/tornadomeet/archive/2012/06/02/2531705.html
前面已經有3篇博文介紹了背景減圖方面相關知識(見下面的鏈接)，在第3篇博文中自己也實現了gmm簡單算法，但效果不是很好，下面來體驗下opencv自帶2個gmm算法。
opencv實現背景減圖法1(codebook和平均背景法)
http://www.cnblogs.com/tornadomeet/archive/2012/04/08/2438158.html
opencv實現背景減圖法2(幀差法)
http://www.cnblogs.com/tornadomeet/archive/2012/05/01/2477629.html
opencv實現背景減圖法3(GMM)
http://www.cnblogs.com/tornadomeet/archive/2012/06/02/2531565.html
工程環境opencv2.3.1+vs2010
實現功能：與上面第三篇博文一樣，完成動態背景的訓練，來檢測前景。
數據來源和前面的一樣： http://research.microsoft.com/en-us/um/people/jckrumm/WallFlower/TestImages.htm 由於該數據是286張bmp格式的圖片，所以用的前200張圖片來作為GMM參數訓練，后186張作為測試。訓練的過程中樹枝被很大幅度的擺動，測試過程中有行人走動，該行人是需要遷就檢測的部分。
Opencv自帶的gmm算法1的實驗結果如下：
其工程代碼如下：
1 // gmm_wavetrees.cpp : 定義控制台應用程序的入口點。
2 //
3
4 #include "stdafx.h"
5
6 #include "opencv2/core/core.hpp"
7 #include "opencv2/video/background_segm.hpp"
8 #include "opencv2/highgui/highgui.hpp"
9 #include "opencv2/imgproc/imgproc.hpp"
10 #include <stdio.h>
11
12 using namespace std;
13 using namespace cv;
14
15 //this is a sample for foreground detection functions
16 string src_img_name="WavingTrees/b00";
17 const char *src_img_name1;
18 Mat img, fgmask, fgimg;
19 int i=-1;
20 char chari[500];
21 bool update_bg_model = true;
22 bool pause=false;
23
24 //第一種gmm,用的是KaewTraKulPong, P. and R. Bowden (2001).
25 //An improved adaptive background mixture model for real-time tracking with shadow detection.
26 BackgroundSubtractorMOG bg_model;
27
28 void refineSegments(const Mat& img, Mat& mask, Mat& dst)
29 {
30 int niters = 3;
31
32 vector<vector<Point> > contours;
33 vector<Vec4i> hierarchy;
34
35 Mat temp;
36
37 dilate(mask, temp, Mat(), Point(-1,-1), niters);//膨脹，3*3的element，迭代次數為niters
38 erode(temp, temp, Mat(), Point(-1,-1), niters*2);//腐蝕
39 dilate(temp, temp, Mat(), Point(-1,-1), niters);
40
41 findContours( temp, contours, hierarchy, CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE );//找輪廓
42
43 dst = Mat::zeros(img.size(), CV_8UC3);
44
45 if( contours.size() == 0 )
46 return;
47
48 // iterate through all the top-level contours,
49 // draw each connected component with its own random color
50 int idx = 0, largestComp = 0;
51 double maxArea = 0;
52
53 for( ; idx >= 0; idx = hierarchy[idx][0] )//這句沒怎么看懂
54 {
55 const vector<Point>& c = contours[idx];
56 double area = fabs(contourArea(Mat(c)));
57 if( area > maxArea )
58 {
59 maxArea = area;
60 largestComp = idx;//找出包含面積最大的輪廓
61 }
62 }
63 Scalar color( 0, 255, 0 );
64 drawContours( dst, contours, largestComp, color, CV_FILLED, 8, hierarchy );
65 }
66
67 int main(int argc, const char** argv)
68 {
69 bg_model.noiseSigma = 10;
70 img=imread("WavingTrees/b00000.bmp");
71 if(img.empty())
72 {
73 namedWindow("image",1);//不能更改窗口
74 namedWindow("foreground image",1);
75 namedWindow("mean background image", 1);
76 }
77 for(;;)
78 {
79 if(!pause)
80 {
81 i++;
82 itoa(i,chari,10);
83 if(i<10)
84 {
85 src_img_name+="00";
86 }
87 else if(i<100)
88 {
89 src_img_name+="0";
90 }
91 else if(i>285)
92 {
93 i=-1;
94 }
95 if(i>=230)
96 update_bg_model=false;
97 else update_bg_model=true;
98
99 src_img_name+=chari;
100 src_img_name+=".bmp";
101
102 img=imread(src_img_name);
103 if( img.empty() )
104 break;
105
106 //update the model
107 bg_model(img, fgmask, update_bg_model ? 0.005 : 0);//計算前景mask圖像，其中輸出fgmask為8-bit二進制圖像，第3個參數為學習速率
108 refineSegments(img, fgmask, fgimg);
109
110 imshow("image", img);
111 imshow("foreground image", fgimg);
112
113 src_img_name="WavingTrees/b00";
114
115 }
116 char k = (char)waitKey(80);
117 if( k == 27 ) break;
118
119 if( k == ' ' )
120 {
121 pause=!pause;
122 }
123 }
124
125 return 0;
126 }
Opencv自帶的gmm算法2的實驗結果如下：
其工程代碼如下：
1 // gmm2_wavetrees.cpp : 定義控制台應用程序的入口點。
2 //
3
4 #include "stdafx.h"
5
6 #include "opencv2/core/core.hpp"
7 #include "opencv2/video/background_segm.hpp"
8 #include "opencv2/highgui/highgui.hpp"
9 #include "opencv2/imgproc/imgproc.hpp"
10 #include <stdio.h>
11
12 using namespace std;
13 using namespace cv;
14
15 //this is a sample for foreground detection functions
16 string src_img_name="WavingTrees/b00";
17 const char *src_img_name1;
18 Mat img, fgmask, fgimg;
19 int i=-1;
20 char chari[500];
21 bool update_bg_model = true;
22 bool pause=false;
23
24 //第一種gmm,用的是KaewTraKulPong, P. and R. Bowden (2001).
25 //An improved adaptive background mixture model for real-time tracking with shadow detection.
26 BackgroundSubtractorMOG2 bg_model;
27
28 void refineSegments(const Mat& img, Mat& mask, Mat& dst)
29 {
30 int niters = 3;
31
32 vector<vector<Point> > contours;
33 vector<Vec4i> hierarchy;
34
35 Mat temp;
36
37 dilate(mask, temp, Mat(), Point(-1,-1), niters);
38 erode(temp, temp, Mat(), Point(-1,-1), niters*2);
39 dilate(temp, temp, Mat(), Point(-1,-1), niters);
40
41 findContours( temp, contours, hierarchy, CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE );
42
43 dst = Mat::zeros(img.size(), CV_8UC3);
44
45 if( contours.size() == 0 )
46 return;
47
48 // iterate through all the top-level contours,
49 // draw each connected component with its own random color
50 int idx = 0, largestComp = 0;
51 double maxArea = 0;
52
53 for( ; idx >= 0; idx = hierarchy[idx][0] )
54 {
55 const vector<Point>& c = contours[idx];
56 double area = fabs(contourArea(Mat(c)));
57 if( area > maxArea )
58 {
59 maxArea = area;
60 largestComp = idx;
61 }
62 }
63 Scalar color( 255, 0, 0 );
64 drawContours( dst, contours, largestComp, color, CV_FILLED, 8, hierarchy );
65 }
66
67 int main(int argc, const char** argv)
68 {
69 img=imread("WvingTrees/b00000.bmp");
70 if(img.empty())
71 {
72 namedWindow("image",1);//不能更改窗口
73 //cvNamedWindow("image",0);
74 namedWindow("foreground image",1);
75 // namedWindow("mean background image", 1);
76 }
77 for(;;)
78 {
79 if(!pause)
80 {
81 i++;
82 itoa(i,chari,10);
83 if(i<10)
84 {
85 src_img_name+="00";
86 }
87 else if(i<100)
88 {
89 src_img_name+="0";
90 }
91 else if(i>285)
92 {
93 i=-1;
94 }
95 // if(i>=230)
96 // update_bg_model=false;
97 // else update_bg_model=true;
98
99 src_img_name+=chari;
100 src_img_name+=".bmp";
101
102 img=imread(src_img_name);
103 if( img.empty() )
104 break;
105
106 //update the model
107 bg_model(img, fgmask, update_bg_model ? 0.005 : 0);//計算前景mask圖像，其中輸出fgmask為8-bit二進制圖像，第3個參數為學習速率
108 refineSegments(img, fgmask, fgimg);
109
110 imshow("foreground image", fgimg);
111 imshow("image", img);
112
113 src_img_name="WavingTrees/b00";
114
115 }
116 char k = (char)waitKey(100);
117 if( k == 27 ) break;
118
119 if( k == ' ' )
120 {
121 pause=!pause;
122 }
123 }
124
125 return 0;
126 }
可以看出gmm1效果比gmm2的好，但是研究發現，gmm2是在gmm1上改進的，不會越該越差吧，除非2個函數的使用方法不同(雖然函數形式一樣),也就是說相同的參數值對函數功能的影響不同。以后有時間在研究了。
前景檢測算法_3(GMM)
http://www.cnblogs.com/tornadomeet/archive/2012/06/02/2531565.html
摘要
本文通過opencv來實現一種前景檢測算法——GMM，算法采用的思想來自論文[1][2][4]。在進行前景檢測前，先對背景進行訓練，對圖像中每個背景采用一個混合高斯模型進行模擬，每個背景的混合高斯的個數可以自適應。然后在測試階段，對新來的像素進行GMM匹配，如果該像素值能夠匹配其中一個高斯，則認為是背景，否則認為是前景。由於整個過程GMM模型在不斷更新學習中，所以對動態背景有一定的魯棒性。最后通過對一個有樹枝搖擺的動態背景進行前景檢測，取得了較好的效果。
關鍵字：GMM,opencv，前景檢測
前言
前景檢測主要分為幀差法，平均背景法，光流法，前景建模法，背景非參數估計，背景建模法等。而本文要實現的方法屬於背景建模法中的一種——GMM，也稱混合高斯模型。
混合高斯模型最早在計算機視覺中的應用是Stauffer et al.[1],作者將其用來做前景檢測，主要是用於視頻監控領域，這個系統和穩定且有自學能力，能在戶外環境跑16個多月。KaewTraKulPong et al.[2]將GMM的訓練過程做了改進，將訓練過程分為2步進行，前L幀采用EM算法進行權值，均值，方差更新，后面的過程就采用[1]中的方法進行更新，取得了更好的檢測效果。Zivkovic et al.在[3]中對GMM理論做了全面的論述，使得GMM理論的使用不僅金限於計算機視覺領域。並且該作者在[4]將該理論進一步具體到背景減圖的前景檢測中來，即加入了參數估計的先驗知識，取得了很好的效果和穩定性。
最近幾年陸續有學者對GMM的背景見圖中的應用做了更深一步的研究，其代表性貢獻見論文[5][6]。
實現過程
本文中主要是根據[2][4]中提出的算法，采用其中的更新方差來根性模型中的3個參數，最后結合用opencv基本底層函數實現該算法。其主要過程如下：
- 首先將每個高斯的均值,方差，權值都設置為0，即初始化個模型矩陣參數。
- 采用視頻中的T幀用來訓練GMM模型。對每一個像素而言，建立其模型個數最大GMM_MAX_COMPONT個高斯的GMM模型。當第一個像素來，單獨為其在程序中設置好其固定的初始均值，方差，並且權值設置為1。
- 非第一幀訓練過程中，當后面來的像素值時，與前面已有的高斯的均值比較，如果該像素點的值與其模型均值差在3倍的方差內，則任務屬於該高斯。此時用如下方程進行更新：
4. 當到達訓練的幀數T后，進行不同像素點GMM個數自適應的選擇。首先用權值除以方差對各個高斯進行從大到小排序，然后選取最前面B個高斯，使
這樣就可以很好的消除訓練過程中的噪聲點。
5. 在測試階段，對新來像素點的值與B個高斯中的每一個均值進行比較，如果其差值在2倍的方差之間的話，則認為是背景，否則認為是前景。並且只要其中有一個高斯分量滿足該條件就認為是前景。前景賦值為255，背景賦值為0。這樣就形成了一副前景二值圖。
6. 由於前景二值圖中含有很多噪聲，所以采用了形態學的開操作將噪聲縮減到0，緊接着用閉操作重建由於開操作丟失的邊緣部分的信息。消除了不連通的小噪聲點。
上面是該算法實現的大概流程，但是當我們在具體編程時，還是有很多細節的地方需要注意，比如有些參數值的選擇。在這里，經過試驗將一些常用的參數值聲明如下：
- 3個參數的值的更新方差中，取其中的學習率為0.005。也就是說T等於200。
- 定義每個像素的最大混合高斯個數取7。
- 取視頻的前200幀進行訓練。
- 取Cf為0.3。即滿足權值大於0.7的個數為自適應高斯的個數B。
- 訓練過程中，需要新建立一個高斯時，其權值取值設為與學習率大小值相等，即0.005。
- 訓練過程中，需要新建立一個高斯時，取該高斯的均值為該輸入像素值大小本身。方差為15。
- 訓練過程中，需要新建立一個高斯時，取該高斯的方差15。
程序流程框圖
該工程的流程框圖如下圖所示：
試驗結果
本次試驗的數據為搖擺的樹枝作為背景，Waving Trees，其來源網址為：http://research.microsoft.com/en-us/um/people/jckrumm/WallFlower/TestImages.htm 由於該數據是286張bmp格式的圖片，所以用的前200張圖片來作為GMM參數訓練，后186張作為測試。訓練的過程中樹枝被很大幅度的擺動，測試過程中有行人走動，該行人是需要遷就檢測的部分。
下圖為訓練過程中動態背景截圖
由上圖可以看出，樹枝在不斷搖擺，可見其背景是動態的。
沒有形態學處理的結果如下：
下圖是有簡單形態學的試驗結果：
(上面2幅圖的左邊為測試原始圖片，右圖為檢測效果圖)
總結
通過本次試驗，不僅學習到了GMM的相關理論知識，以及背景減圖法在前景檢測中的應用。更重要的時，對opencv在計算機視覺中的使用更進一步的熟悉了。另外對於目標檢測的難點有了更深一層的理解。
參考文獻
- Stauffer, C. and W. E. L. Grimson (1999). Adaptive background mixture models for real-time tracking, IEEE.
- KaewTraKulPong, P. and R. Bowden (2001). An improved adaptive background mixture model for real-time tracking with shadow detection.
- Zivkovic, Z. and F. van der Heijden (2004). "Recursive unsupervised learning of finite mixture models." Pattern Analysis and Machine Intelligence, IEEE Transactions on 26(5): 651-656.
- Zivkovic, Z. and F. van der Heijden (2006). "Efficient adaptive density estimation per image pixel for the task of background subtraction." Pattern recognition letters 27(7): 773-780.
- Bouzerdoum, A., A. Beghdadi, et al. (2010). "On the analysis of background subtraction techniques using Gaussian mixture models."
- Lin, H. H., J. H. Chuang, et al. (2011). "Regularized background adaptation: a novel learning rate control scheme for Gaussian mixture modeling." Image Processing, IEEE Transactions on 20(3): 822-836.
后續改進
需要改進的地方：
1. 程序運行的速度太慢，很多參數都是浮點數，每個像素都要建立一個gmm，gmm個數本身比較多，所以訓練過程中速度比較慢，代碼需要優化。
2. 最后生成的前景圖需要用連通域處理算法進行修整，即需要形態學操作，然后找出連通域大小滿足要求的輪廓，用多邊形擬合來進行處理。這種算法在《learning opencv》一書中有提到。后續有時間加入該算法，效果會好很多的。
附錄：工程代碼
1 //gmm.cpp : 定義控制台應用程序的入口點。
2 #include "stdafx.h"
3 #include "cv.h"
4 #include "highgui.h"
5 #include <opencv2/imgproc/imgproc.hpp>
6 #include <opencv2/core/core.hpp>
7 #include <opencv2/highgui/highgui.hpp>
8 #include <stdio.h>
9 #include <iostream>
10
11 using namespace cv;
12 using namespace std;
13
14 //定義gmm模型用到的變量
15 #define GMM_MAX_COMPONT 6
16 #define GMM_LEARN_ALPHA 0.005 //該學習率越大的話，學習速度太快，效果不好
17 #define GMM_THRESHOD_SUMW 0.7 //如果取值太大了的話，則更多樹的部分都被檢測出來了
18 #define END_FRAME 200
19
20 bool pause=false;
21
22 Mat w[GMM_MAX_COMPONT];
23 Mat u[GMM_MAX_COMPONT];
24 Mat sigma[GMM_MAX_COMPONT];
25 Mat fit_num,gmask,foreground;
26 vector<Mat> output_m;
27 Mat output_img;
28
29 float temp_w,temp_sigma;
30 unsigned char temp_u;
31 int i=-1;
32
33
34 //For connected components:
35 int CVCONTOUR_APPROX_LEVEL = 2; // Approx.threshold - the bigger it is, the simpler is the boundary
36 int CVCLOSE_ITR = 1;
37
38 //Just some convienience macros
39 #define CV_CVX_WHITE CV_RGB(0xff,0xff,0xff)
40 #define CV_CVX_BLACK CV_RGB(0x00,0x00,0x00)
41
42 //gmm整體初始化函數聲明
43 void gmm_init(Mat img);
44
45 //gmm第一幀初始化函數聲明
46 void gmm_first_frame(Mat img);
47
48 //gmm訓練過程函數聲明
49 void gmm_train(Mat img);
50
51 //對輸入圖像每個像素gmm選擇合適的個數函數聲明
52 void gmm_fit_num(Mat img);
53
54 //gmm測試函數的實現
55 void gmm_test(Mat img);
56
57 //連通域去噪函數聲明
58 void find_connected_components(Mat img);
59 //void cvconnectedComponents(IplImage *mask, int poly1_hull0, float perimScale, int *num, CvRect *bbs, CvPoint *centers);
60
61 int main(int argc, const char* argv[])
62 {
63 Mat img,img_gray;
64 char str_num[5];
65
66
67 // char *str_num;//why does this definition not work?
68 String str="WavingTrees/b00";//string,the 's' can be a captial or lower-caseletters
69
70 /****read the first image,and reset the array w,u,sigma****/
71 img=imread("WavingTrees/b00000.bmp");
72 if(img.empty())
73 return -1;
74
75 output_img=Mat::zeros(img.size(),img.type());
76 cvtColor(img,img_gray,CV_BGR2GRAY);//covert the colorful image to the corresponding gray-level image
77
78 /****initialization the three parameters ****/
79 gmm_init(img_gray);
80 fit_num=Mat(img.size(),CV_8UC1,-1);//初始化為1
81 gmask=Mat(img.size(),CV_8UC1,-1);
82 foreground=img.clone();
83 split(img,output_m);
84
85 output_m[0]=Mat::zeros(img.size(),output_m[0].type());
86 output_m[1]=Mat::zeros(img.size(),output_m[0].type());
87 output_m[2]=Mat::zeros(img.size(),output_m[0].type());
88
89 namedWindow("src",WINDOW_AUTOSIZE);
90 namedWindow("gmask",WINDOW_AUTOSIZE);
91
92 //在定義視頻輸出對象時，文件名一定后面要加后綴，比如這里的.avi,否則是輸出不了視頻的！並且這里只能是avi格式的，當參數為('P','I','M','1')時
93 VideoWriter output_src("src.avi",CV_FOURCC('P','I','M','1'),20,Size(160,120),1);//c++版本的opencv用Size函數即可，c版本的用cvSize()函數
94 //VideoWriter output_src("src.avi",CV_FOURCC('M','J','P','G'),5,Size(160,120),1);//c++版本的opencv用Size函數即可，c版本的用cvSize()函數
95 VideoWriter output_dst("dst.avi",CV_FOURCC('P','I','M','1'),20,Size(160,120),1);//這樣輸出的是3個通道的數據
96 while(1)
97
98 {
99 if(!pause)
100 {
101 /****read image from WavingTrees****/
102 i++;
103 _itoa_s(i,str_num,10);//the latest name is _itoa_s or _itoa,not the itoa,although iota can be used,deprecated
104 if(i<10)
105 str+="00";
106 else if(i<100)
107 str+="0";
108 else if(i>285)//we used the first 285 frames to learn the gmm model
109 i=-1;
110 str+=str_num;
111 str+=".bmp";
112
113 img=imread(str);
114 if(img.empty())
115 break;
116 str="WavingTrees/b00";//after read,str must be reseted ;
117
118 cvtColor(img,img_gray,CV_BGR2GRAY);//covert the colorful image to the corresponding gray-level image
119
120 /****when it is the first frame,set the default parameter****/
121 if(1==i)
122 {
123 gmm_first_frame(img_gray);
124 }
125
126 //the train of gmm phase
127 //if(1<i&&i<5&&i!=3)//由此可知當i大於等於3以后，就會一直出現錯誤，且錯誤在內部排序的部分
128 if(1<i<END_FRAME)
129 {
130 gmm_train(img_gray);
131 }//end the train phase
132
133 cout<<i<<endl;
134 /****chose the fitting number of component in gmm****/
135 if(END_FRAME==i)
136 {
137 gmm_fit_num(img_gray);
138 // cout<<fit_num<<endl;//其輸出值有4個高斯的，但也有0個高斯的，why？照理說不可能的啊！
139 }
140
141 /****start the test phase****/
142 if(i>=END_FRAME)
143 {
144 output_src<<img;
145 gmm_test(img_gray);
146 find_connected_components(img_gray);
147
148 output_m[0]=gmask.clone();
149 output_m[1]=gmask.clone();
150 output_m[2]=gmask.clone();
151
152 merge(output_m,output_img);
153 output_dst<<output_img;
154 }
155 if(285==i)
156 {
157 return 0;
158 }
159 imshow("src",img);
160 imshow("gmask",gmask);
161 }
162
163 char c=(char)waitKey(1);
164 if(c==27)//if press the ESC key,the exit the proggram
165 break;
166 if(c==' ')
167 // pause=~pause;//if use '~',then the pause key cannot work,why?
168 pause=!pause;
169 }
170 return 0;
171 }
172
173
174 //gmm初始化函數實現
175 void gmm_init(Mat img)
176 {
177 /****initialization the three parameters ****/
178 for(int j=0;j<GMM_MAX_COMPONT;j++)
179 {
180 w[j]=Mat(img.size(),CV_32FC1,0.0);//CV_32FC1本身體現了正負符號
181 u[j]=Mat(img.size(),CV_8UC1,-1);//為什么這里賦值為0時，后面的就一直出錯？暫時還不知道原因，先賦值-1，其實內部存儲的也是0
182 sigma[j]=Mat(img.size(),CV_32FC1,0.0);//float類型
183 }
184
185 //為什么一下語句不能放在這個函數里面呢
186 // output_m[0]=Mat(img.size(),CV_8UC1,0);
187 // output_m[1]=Mat(img.size(),CV_8UC1,0);
188 // output_m[2]=Mat(img.size(),CV_8UC1,0);
189 }
190
191
192 //gmm第一幀初始化函數實現
193 void gmm_first_frame(Mat img)
194 {
195 for(int m=0;m<img.rows;m++)
196 for(int n=0;n<img.cols;n++)
197 {
198 w[0].at<float>(m,n)=1.0;
199
200 //if the pixvel is gray-clever,then we should use unsigned char,not the unsigned int
201 u[0].at<unsigned char>(m,n)=img.at<unsigned char>(m,n);// 一定要注意其類型轉換，否則會得不得預期的結果
202 sigma[0].at<float>(m,n)=15.0;//opencv 自帶的gmm代碼中用的是15.0
203
204 for(int k=1;k<GMM_MAX_COMPONT;k++)
205 {
206 /****when assigment this,we must be very carefully****/
207 w[k].at<float>(m,n)=0.0;
208 u[k].at<unsigned char>(m,n)=-1;
209 sigma[k].at<float>(m,n)=15.0;//防止后面排序時有分母為0的情況
210 }
211 }
212 }
213
214
215 //gmm訓練過程函數實現
216 void gmm_train(Mat img)
217 {
218 for(int m=0;m<img.rows;m++)
219 for(int n=0;n<img.cols;n++)
220 {
221 int k=0;
222 int nfit=0;
223 for(;k<GMM_MAX_COMPONT;k++)
224 {
225 // if(w[k].at<float>(m,n)!=0)//只有在權值不為0的情況下才進行比較
226 // {
227 int delam=abs(img.at<unsigned char>(m,n)-u[k].at<unsigned char>(m,n));//防止溢出
228 long dis=delam*delam;
229 if(dis<3.0*sigma[k].at<float>(m,n))//the present pixpel is fit the component
230 {
231 /****update the weight****/
232 w[k].at<float>(m,n)=w[k].at<float>(m,n)+GMM_LEARN_ALPHA*(1-w[k].at<float>(m,n));
233
234 /****update the average****/
235 u[k].at<unsigned char>(m,n)=u[k].at<unsigned char>(m,n)+(GMM_LEARN_ALPHA/w[k].at<float>(m,n))*delam;
236
237 /****update the variance****/
238 sigma[k].at<float>(m,n)=sigma[k].at<float>(m,n)+(GMM_LEARN_ALPHA/w[k].at<float>(m,n))*(dis-sigma[k].at<float>(m,n));
239
240 // break;
241 }
242 else{
243 w[k].at<float>(m,n)=w[k].at<float>(m,n)+GMM_LEARN_ALPHA*(0-w[k].at<float>(m,n));
244 nfit++;
245 }
246 // }
247 }
248
249 ////訓練過程加速算法
250 //for(int bb=k+1;bb<GMM_MAX_COMPONT;bb++)
251 //{
252 // w[bb].at<float>(m,n)=w[bb].at<float>(m,n)+GMM_LEARN_ALPHA*(0-w[bb].at<float>(m,n));
253 // nfit++;
254 //}
255
256 //對gmm各個高斯進行排序,從大到小排序,排序依據為w/sigma
257 for(int kk=0;kk<GMM_MAX_COMPONT;kk++)
258 {
259 for(int rr=kk;rr<GMM_MAX_COMPONT;rr++)
260 {
261 //怎樣才能做到gmm結構體整體排序呢？
262 if(w[rr].at<float>(m,n)/sigma[rr].at<float>(m,n)>w[kk].at<float>(m,n)/sigma[kk].at<float>(m,n))
263 {
264 //權值交換
265 temp_w=w[rr].at<float>(m,n);
266 w[rr].at<float>(m,n)=w[kk].at<float>(m,n);
267 w[kk].at<float>(m,n)=temp_w;
268
269 //均值交換
270 temp_u=u[rr].at<unsigned char>(m,n);
271 u[rr].at<unsigned char>(m,n)=u[kk].at<unsigned char>(m,n);
272 u[kk].at<unsigned char>(m,n)=temp_u;
273
274 //方差交換
275 temp_sigma=sigma[rr].at<float>(m,n);
276 sigma[rr].at<float>(m,n)=sigma[kk].at<float>(m,n);
277 sigma[kk].at<float>(m,n)=temp_sigma;
278 }
279 }
280 }
281
282 //****如果沒有滿足條件的高斯，則重新開始算一個高斯分布****/
283 if(nfit==GMM_MAX_COMPONT&&0==w[GMM_MAX_COMPONT-1].at<float>(m,n))//if there is no exit component fit,then start a new componen
284 {
285 //不能寫為for(int h=0;h<MAX_GMM_COMPONT&&((0==w[h].at<float>(m,n)));h++)，因為這樣明顯h不會每次都加1
286 for(int h=0;h<GMM_MAX_COMPONT;h++)
287 {
288 if((0==w[h].at<float>(m,n)))
289 {
290 w[h].at<float>(m,n)=GMM_LEARN_ALPHA;//按照論文的參數來的
291 u[h].at<unsigned char>(m,n)=(unsigned char)img.at<unsigned char>(m,n);
292 sigma[h].at<float>(m,n)=15.0;//the opencv library code is 15.0
293 for(int q=0;q<GMM_MAX_COMPONT&&q!=h;q++)
294 {
295 /****update the other unfit's weight,u and sigma remain unchanged****/
296 w[q].at<float>(m,n)*=1-GMM_LEARN_ALPHA;//normalization the weight,let they sum to 1
297 }
298 break;//找到第一個權值不為0的即可
299 }
300 }
301 }
302 //如果GMM_MAX_COMPONT都曾經賦值過，則用新來的高斯代替權值最弱的高斯，權值不變，只更新均值和方差
303 else if(nfit==GMM_MAX_COMPONT&&w[GMM_MAX_COMPONT-1].at<float>(m,n)!=0)
304 {
305 u[GMM_MAX_COMPONT-1].at<unsigned char>(m,n)=(unsigned char)img.at<unsigned char>(m,n);
306 sigma[GMM_MAX_COMPONT-1].at<float>(m,n)=15.0;//the opencv library code is 15.0
307 }
308
309
310 }
311 }//end the train phase
312
313
314 //對輸入圖像每個像素gmm選擇合適的個數
315 void gmm_fit_num(Mat img)
316 {
317 //float sum_w=0.0;//重新賦值為0，給下一個像素做累積
318 for(int m=0;m<img.rows;m++)
319 for(int n=0;n<img.cols;n++)
320 {
321 float sum_w=0.0;//重新賦值為0，給下一個像素做累積
322 //chose the fittest number fit_num
323 for(unsigned char a=0;a<GMM_MAX_COMPONT;a++)
324 {
325 //cout<<w[a].at<float>(m,n)<<endl;
326 sum_w+=w[a].at<float>(m,n);
327 if(sum_w>=GMM_THRESHOD_SUMW)//如果這里THRESHOD_SUMW=0.6的話，那么得到的高斯數目都為1，因為每個像素都有一個權值接近1
328 {
329 fit_num.at<unsigned char>(m,n)=a+1;
330 break;
331 }
332 }
333 }
334 }
335
336
337 //gmm測試函數的實現
338 void gmm_test(Mat img)
339 {
340 for(int m=0;m<img.rows;m++)
341 for(int n=0;n<img.cols;n++)
342 {
343 unsigned char a=0;
344 for(;a<fit_num.at<unsigned char>(m,n);a++)
345 {
346 //如果對sigma取根號的話，樹枝當做前景的概率會更大，不過人被檢測出來的效果也更好些；用2相乘，不用開根號效果還不錯
347 // if(abs(img.at<unsigned char>(m,n)-u[a].at<unsigned char>(m,n))<(unsigned char)(2*(sigma[a].at<float>(m,n))))
348 if(abs(img.at<unsigned char>(m,n)-u[a].at<unsigned char>(m,n))<(unsigned char)(2.5*(sigma[a].at<float>(m,n))))
349 {
350 gmask.at<unsigned char>(m,n)=1;//背景
351 break;
352 }
353 }
354 if(a==fit_num.at<unsigned char>(m,n))
355 gmask.at<unsigned char>(m,n)=255;//前景
356 }
357 }
358
359 //連通域去噪函數實現
360 void find_connected_components(Mat img)
361 {
362 morphologyEx(gmask,gmask,MORPH_OPEN,Mat());
363 // morphologyEx(gmask,gmask,MORPH_CLOSE,Mat());
364 }
365
366 ////連通域去噪函數實現
367 //void find_connected_components(Mat img)
368 //{
369 // morphologyEx(gmask,gmask,MORPH_OPEN,Mat());
370 // morphologyEx(gmask,gmask,MORPH_CLOSE,Mat());
371 //// erode(gmask,gmask,Mat());//只腐蝕是不行的，人來了也被腐蝕掉了
372 //
373 // vector<vector<Point>> contours;//由點向量組成的向量，所以有2個層次結構
374 // vector<Vec4i> hierarchy;//typedef Vec<int,4>Vec4i;即由4個整數組成的向量
375 //
376 // //找到gmask的輪廓，存儲在contours中，其拓撲結構存儲在hierarchy中，且僅僅找出最外面的輪廓，用壓縮算法只存儲水平，垂直，斜對角線的端點
377 // //其中暗含了hierarchy[i][2]=hierarchy[3]=-1,即其父輪廓和嵌套輪廓不用考慮
378 // findContours(gmask,contours,hierarchy,CV_RETR_EXTERNAL,CV_CHAIN_APPROX_SIMPLE);
379 // if(contours.size()==0)
380 // return;
381 //
382 // int idex=0;
383 // for(idex=0;idex<contours.size();idex++)
384 // {
385 // const vector<Point>& c=contours[idex];
386 //// const vector<Point>& cnull::zeros();
387 // double len=arcLength(c,false);//求出輪廓的周長，並不一定要求其是封閉的
388 // double q=(img.rows+img.cols)/4;
389 // if(q>=len)
390 // {
391 // const vector<Point> &cnew=contours[idex];
392 // // Mat mcnew=Mat(cnew);
393 // // Mat mcnew;
394 // // approxPolyDP(Mat(c),mcnew,CVCONTOUR_APPROX_LEVEL,false);//多邊形曲線擬合，並不一定要求其輪廓閉合
395 // // approxPolyDP(Mat(c),Mat(cnew),CVCONTOUR_APPROX_LEVEL,false);//多邊形曲線擬合，並不一定要求其輪廓閉合
396 // approxPolyDP(Mat(c),cnew,CVCONTOUR_APPROX_LEVEL,false);//多邊形曲線擬合，並不一定要求其輪廓閉合
397 // // cnew=vector<Point>(mcnew);
398 // // contours[idex]=cnew;
399 // }
400 //// else contours[idex]=vector<Point(0,0,0)>;
401 // }
402 //
403 //}
404
405 ///////////////////////////////////////////////////////////////////////////////////////////
406 //void cvconnectedComponents(IplImage *mask, int poly1_hull0, float perimScale, int *num, CvRect *bbs, CvPoint *centers)
407 // This cleans up the forground segmentation mask derived from calls to cvbackgroundDiff
408 //
409 // mask Is a grayscale (8 bit depth) "raw" mask image which will be cleaned up
410 //
411 // OPTIONAL PARAMETERS:
412 // poly1_hull0 If set, approximate connected component by (DEFAULT) polygon, or else convex hull (0)
413 // perimScale Len = image (width+height)/perimScale. If contour len < this, delete that contour (DEFAULT: 4)
414 // num Maximum number of rectangles and/or centers to return, on return, will contain number filled (DEFAULT: NULL)
415 // bbs Pointer to bounding box rectangle vector of length num. (DEFAULT SETTING: NULL)
416 // centers Pointer to contour centers vectore of length num (DEFULT: NULL)
417 //
418 //void cvconnectedComponents(IplImage *mask, int poly1_hull0, float perimScale, int *num, CvRect *bbs, CvPoint *centers)
419 //{
420 // static CvMemStorage* mem_storage = NULL;
421 // static CvSeq* contours = NULL;
422 //// static CvSeq** firstContour;
423 //
424 // //CLEAN UP RAW MASK
425 // //開運算作用：平滑輪廓，去掉細節,斷開缺口
426 // cvMorphologyEx( mask, mask, NULL, NULL, CV_MOP_OPEN, CVCLOSE_ITR );//對輸入mask進行開操作，CVCLOSE_ITR為開操作的次數，輸出為mask圖像
427 // //閉運算作用：平滑輪廓，連接缺口
428 // cvMorphologyEx( mask, mask, NULL, NULL, CV_MOP_CLOSE, CVCLOSE_ITR );//對輸入mask進行閉操作，CVCLOSE_ITR為閉操作的次數，輸出為mask圖像
429 //
430 // //FIND CONTOURS AROUND ONLY BIGGER REGIONS
431 // if( mem_storage==NULL ) mem_storage = cvCreateMemStorage(0);
432 // else cvClearMemStorage(mem_storage);
433 //
434 // //CV_RETR_EXTERNAL=0是在types_c.h中定義的，CV_CHAIN_APPROX_SIMPLE=2也是在該文件中定義的
435 // CvContourScanner scanner = cvStartFindContours(mask,mem_storage,sizeof(CvContour),CV_RETR_EXTERNAL,CV_CHAIN_APPROX_SIMPLE);
436 //// CvContourScanner scanner = cvFindContours(mask,mem_storage,firstContour,sizeof(CvContour),CV_RETR_EXTERNAL,CV_CHAIN_APPROX_SIMPLE);
437 // CvSeq* c;
438 // int numCont = 0;
439 // while( (c = cvFindNextContour( scanner )) != NULL )
440 // {
441 // double len = cvContourPerimeter( c );
442 // double q = (mask->height + mask->width) /perimScale; //calculate perimeter len threshold
443 // if( len < q ) //Get rid of blob if it's perimeter is too small
444 // {
445 // cvSubstituteContour( scanner, NULL );
446 // }
447 // else //Smooth it's edges if it's large enough
448 // {
449 // CvSeq* c_new;
450 // if(poly1_hull0) //Polygonal approximation of the segmentation
451 // c_new = cvApproxPoly(c,sizeof(CvContour),mem_storage,CV_POLY_APPROX_DP, CVCONTOUR_APPROX_LEVEL,0);
452 // else //Convex Hull of the segmentation
453 // c_new = cvConvexHull2(c,mem_storage,CV_CLOCKWISE,1);
454 // cvSubstituteContour( scanner, c_new );
455 // numCont++;
456 // }
457 // }
458 // contours = cvEndFindContours( &scanner );
459 //
460 // // PAINT THE FOUND REGIONS BACK INTO THE IMAGE
461 // cvZero( mask );
462 // IplImage *maskTemp;
463 // //CALC CENTER OF MASS AND OR BOUNDING RECTANGLES
464 // if(num != NULL)
465 // {
466 // int N = *num, numFilled = 0, i=0;
467 // CvMoments moments;
468 // double M00, M01, M10;
469 // maskTemp = cvCloneImage(mask);
470 // for(i=0, c=contours; c != NULL; c = c->h_next,i++ )
471 // {
472 // if(i < N) //Only process up to *num of them
473 // {
474 // cvDrawContours(maskTemp,c,CV_CVX_WHITE, CV_CVX_WHITE,-1,CV_FILLED,8);
475 // //Find the center of each contour
476 // if(centers != NULL)
477 // {
478 // cvMoments(maskTemp,&moments,1);
479 // M00 = cvGetSpatialMoment(&moments,0,0);
480 // M10 = cvGetSpatialMoment(&moments,1,0);
481 // M01 = cvGetSpatialMoment(&moments,0,1);
482 // centers[i].x = (int)(M10/M00);
483 // centers[i].y = (int)(M01/M00);
484 // }
485 // //Bounding rectangles around blobs
486 // if(bbs != NULL)
487 // {
488 // bbs[i] = cvBoundingRect(c);
489 // }
490 // cvZero(maskTemp);
491 // numFilled++;
492 // }
493 // //Draw filled contours into mask
494 // cvDrawContours(mask,c,CV_CVX_WHITE,CV_CVX_WHITE,-1,CV_FILLED,8); //draw to central mask
495 // } //end looping over contours
496 // *num = numFilled;
497 // cvReleaseImage( &maskTemp);
498 // }
499 // //ELSE JUST DRAW PROCESSED CONTOURS INTO THE MASK
500 // else
501 // {
502 // for( c=contours; c != NULL; c = c->h_next )
503 // {
504 // cvDrawContours(mask,c,CV_CVX_WHITE, CV_CVX_BLACK,-1,CV_FILLED,8);
505 // }
506 // }
507 //}
計算機視覺小菜鳥的專欄
http://www.cnblogs.com/tornadomeet/archive/2012/06/02/2531565.html
運動目標檢測可以分為攝像機固定和攝像機運動兩類；對於攝像機運動情況下的運動目標檢測，光流法是比較常用的解決方法，通過求解偏微分方程求得圖像序列的光流場，從而預測攝像機的運動狀態。對於攝像機固定的情形，可以采用光流法也可以采用高斯背景模型，考慮到光流法計算量巨大，故而，高斯背景模型相對更常用一些。需要提醒的是，這里所謂的“背景”是指用戶不需要的目標，而“前景”自然指代用戶需要的特定目標了。背景模型有很多種，其中很多方法對光照的的突變和其它因素的適應能力不夠，而高斯混合模型是最為成功的一種背景建模方法。
高斯背景模型是由Stauffer等人提出的經典的自適應混合高斯背景提取方法，是一種基於背景建模的方法，它是根據視頻中的每個像素在時域上的分布情況來構建各個像素的顏色分布模型，依次來達到背景建模的目的。混合高斯背景模型是有限個高斯函數的加權和，它能描述像素的多峰狀態，適用於對光照漸變、樹木搖擺等復雜背景進行准確建模。此后經過很多研究人員的不斷改進，該方法目前已經成為比較常用的背景提取方法。
