c++版本的高斯混合模型的源代码完全注释

之前看到过C版本的，感觉写的很长，没有仔细看，但是C++版本的写的还是很不错的。我仔细看了一下，并对内容进行了仔细的注释，如果有人没有看懂，欢迎留言讨论。

先看一眼头文件，在background_segm.hpp中

[cpp] view plain copy

print ?

1. class CV_EXPORTS_W BackgroundSubtractorMOG : public BackgroundSubtractor  
2. {  
3. public:  
4.     //! the default constructor  
5.     CV_WRAP BackgroundSubtractorMOG();  
6.     //! the full constructor that takes the length of the history, the number of gaussian mixtures, the background ratio parameter and the noise strength  
7.     CV_WRAP BackgroundSubtractorMOG(int history, int nmixtures, double backgroundRatio, double noiseSigma=0);  
8.     //! the destructor  
9.     virtual ~BackgroundSubtractorMOG();  
10.     //! the update operator  
11.     virtual void operator()(InputArray image, OutputArray fgmask, double learningRate=0);  
12.       
13.     //! re-initiaization method  
14.     virtual void initialize(Size frameSize, int frameType);  
15.       
16.     virtual AlgorithmInfo* info() const;  
17.   
18. protected:      
19.     Size frameSize;  
20.     int frameType;  
21.     Mat bgmodel;  
22.     int nframes;  
23.     int history;//利用历史帧数计算学习速率，不是主要参数  
24.     int nmixtures;//高斯模型的个数  
25.     double varThreshold;//方差门限  
26.     double backgroundRatio;//背景门限  
27.     double noiseSigma;//噪声方差  
28. };    

class CV_EXPORTS_W BackgroundSubtractorMOG : public BackgroundSubtractor
{
public:
    //! the default constructor
    CV_WRAP BackgroundSubtractorMOG();
    //! the full constructor that takes the length of the history, the number of gaussian mixtures, the background ratio parameter and the noise strength
    CV_WRAP BackgroundSubtractorMOG(int history, int nmixtures, double backgroundRatio, double noiseSigma=0);
    //! the destructor
    virtual ~BackgroundSubtractorMOG();
    //! the update operator
    virtual void operator()(InputArray image, OutputArray fgmask, double learningRate=0);
    
    //! re-initiaization method
    virtual void initialize(Size frameSize, int frameType);
    
    virtual AlgorithmInfo* info() const;

protected:    
    Size frameSize;
    int frameType;
    Mat bgmodel;
    int nframes;
    int history;//利用历史帧数计算学习速率，不是主要参数
    int nmixtures;//高斯模型的个数
    double varThreshold;//方差门限
    double backgroundRatio;//背景门限
    double noiseSigma;//噪声方差
};	

再看一眼源文件，在bgfg_gaussmix.cpp中：

[cpp] view plain copy

print ?

1. static const int defaultNMixtures = 5;//默认混合模型个数  
2. static const int defaultHistory = 200;//默认历史帧数  
3. static const double defaultBackgroundRatio = 0.7;//默认背景门限  
4. static const double defaultVarThreshold = 2.5*2.5;//默认方差门限  
5. static const double defaultNoiseSigma = 30*0.5;//默认噪声方差  
6. static const double defaultInitialWeight = 0.05;//默认初始权值  
7.  //不带参数的构造函数，使用默认值    
8. BackgroundSubtractorMOG::BackgroundSubtractorMOG()  
9. {  
10.     frameSize = Size(0,0);  
11.     frameType = 0;  
12.       
13.     nframes = 0;  
14.     nmixtures = defaultNMixtures;  
15.     history = defaultHistory;  
16.     varThreshold = defaultVarThreshold;  
17.     backgroundRatio = defaultBackgroundRatio;  
18.     noiseSigma = defaultNoiseSigma;  
19. }  
20. //带参数的构造函数，使用参数传进来的值      
21. BackgroundSubtractorMOG::BackgroundSubtractorMOG(int _history, int _nmixtures,  
22.                                                  double _backgroundRatio,  
23.                                                  double _noiseSigma)  
24. {  
25.     frameSize = Size(0,0);  
26.     frameType = 0;  
27.       
28.     nframes = 0;  
29.     nmixtures = min(_nmixtures > 0 ? _nmixtures : defaultNMixtures, 8);//不能超过8个，否则就用默认的  
30.     history = _history > 0 ? _history : defaultHistory;//不能小于0，否则就用默认的  
31.     varThreshold = defaultVarThreshold;//门限使用默认的  
32.     backgroundRatio = min(_backgroundRatio > 0 ? _backgroundRatio : 0.95, 1.);//背景门限必须大于0，小于1，否则使用0.95  
33.     noiseSigma = _noiseSigma <= 0 ? defaultNoiseSigma : _noiseSigma;//噪声门限大于0  
34. }  
35.       
36. BackgroundSubtractorMOG::~BackgroundSubtractorMOG()  
37. {  
38. }  
39.   
40.   
41. void BackgroundSubtractorMOG::initialize(Size _frameSize, int _frameType)  
42. {  
43.     frameSize = _frameSize;  
44.     frameType = _frameType;  
45.     nframes = 0;  
46.       
47.     int nchannels = CV_MAT_CN(frameType);  
48.     CV_Assert( CV_MAT_DEPTH(frameType) == CV_8U );  
49.       
50.     // for each gaussian mixture of each pixel bg model we store ...  
51.     // the mixture sort key (w/sum_of_variances), the mixture weight (w),  
52.     // the mean (nchannels values) and  
53.     // the diagonal covariance matrix (another nchannels values)  
54.     bgmodel.create( 1, frameSize.height*frameSize.width*nmixtures*(2 + 2*nchannels), CV_32F );//初始化一个1行*XX列的矩阵  
55.     //XX是这样计算的：图像的行*列*混合模型的个数*（1（优先级） + 1（权值） + 2（均值 + 方差） * 通道数）  
56.     bgmodel = Scalar::all(0);//清零  
57. }  
58.   
59. //设为模版，就是考虑到了彩色图像与灰度图像两种情况      
60. template<typename VT> struct MixData  
61. {  
62.     float sortKey;  
63.     float weight;  
64.     VT mean;  
65.     VT var;  
66. };  
67.   
68.       
69. static void process8uC1( const Mat& image, Mat& fgmask, double learningRate,  
70.                          Mat& bgmodel, int nmixtures, double backgroundRatio,  
71.                          double varThreshold, double noiseSigma )  
72. {  
73.     int x, y, k, k1, rows = image.rows, cols = image.cols;  
74.     float alpha = (float)learningRate, T = (float)backgroundRatio, vT = (float)varThreshold;//学习速率、背景门限、方差门限  
75.     int K = nmixtures;//混合模型个数  
76.     MixData<float>* mptr = (MixData<float>*)bgmodel.data;  
77.       
78.     const float w0 = (float)defaultInitialWeight;//初始权值  
79.     const float sk0 = (float)(w0/(defaultNoiseSigma*2));//初始优先级  
80.     const float var0 = (float)(defaultNoiseSigma*defaultNoiseSigma*4);//初始方差  
81.     const float minVar = (float)(noiseSigma*noiseSigma);//最小方差  
82.       
83.     for( y = 0; y < rows; y++ )  
84.     {  
85.         const uchar* src = image.ptr<uchar>(y);  
86.         uchar* dst = fgmask.ptr<uchar>(y);  
87.           
88.         if( alpha > 0 )//如果学习速率为0，则退化为背景相减  
89.         {  
90.             for( x = 0; x < cols; x++, mptr += K )  
91.             {  
92.                 float wsum = 0;  
93.                 float pix = src[x];//每个像素  
94.                 int kHit = -1, kForeground = -1;//是否属于模型，是否属于前景  
95.                   
96.                 for( k = 0; k < K; k++ )//每个高斯模型  
97.                 {  
98.                     float w = mptr[k].weight;//当前模型的权值  
99.                     wsum += w;//权值累加  
100.                     if( w < FLT_EPSILON )  
101.                         break;  
102.                     float mu = mptr[k].mean;//当前模型的均值  
103.                     float var = mptr[k].var;//当前模型的方差  
104.                     float diff = pix - mu;//当前像素与模型均值之差  
105.                     float d2 = diff*diff;//平方  
106.                     if( d2 < vT*var )//是否小于方门限，即是否属于该模型  
107.                     {  
108.                         wsum -= w;//如果匹配，则把它减去，因为之后会更新它  
109.                         float dw = alpha*(1.f - w);  
110.                         mptr[k].weight = w + dw;//增加权值  
111.                         //注意源文章中涉及概率的部分多进行了简化，将概率变为1  
112.                         mptr[k].mean = mu + alpha*diff;//修正均值         
113.                         var = max(var + alpha*(d2 - var), minVar);//开始时方差清零0，所以这里使用噪声方差作为默认方差，否则使用上一次方差  
114.                         mptr[k].var = var;//修正方差  
115.                         mptr[k].sortKey = w/sqrt(var);//重新计算优先级，貌似这里不对，应该使用更新后的mptr[k].weight而不是w  
116.                           
117.                         for( k1 = k-1; k1 >= 0; k1-- )//从匹配的第k个模型开始向前比较，如果更新后的单高斯模型优先级超过他前面的那个，则交换顺序  
118.                         {  
119.                             if( mptr[k1].sortKey >= mptr[k1+1].sortKey )//如果优先级没有发生改变，则停止比较  
120.                                 break;  
121.                             std::swap( mptr[k1], mptr[k1+1] );//交换它们的顺序，始终保证优先级最大的在前面  
122.                         }  
123.                           
124.                         kHit = k1+1;//记录属于哪个模型  
125.                         break;  
126.                     }  
127.                 }  
128.                   
129.                 if( kHit < 0 ) // no appropriate gaussian mixture found at all, remove the weakest mixture and create a new one  
130.                                 //当前像素不属于任何一个模型  
131.                 {  
132.                     //初始化一个新模型  
133.                     kHit = k = min(k, K-1);//有两种情况，当最开始的初始化时，k并不是等于K-1的  
134.                     wsum += w0 - mptr[k].weight;//从权值总和中减去原来的那个模型，并加上默认权值  
135.                     mptr[k].weight = w0;//初始化权值  
136.                     mptr[k].mean = pix;//初始化均值  
137.                     mptr[k].var = var0; //初始化方差  
138.                     mptr[k].sortKey = sk0;//初始化权值  
139.                 }  
140.                 else  
141.                     for( ; k < K; k++ )  
142.                         wsum += mptr[k].weight;//求出剩下的总权值  
143.                   
144.                 float wscale = 1.f/wsum;//归一化  
145.                 wsum = 0;  
146.                 for( k = 0; k < K; k++ )  
147.                 {  
148.                     wsum += mptr[k].weight *= wscale;  
149.                     mptr[k].sortKey *= wscale;//计算归一化权值  
150.                     if( wsum > T && kForeground < 0 )  
151.                         kForeground = k+1;//第几个模型之后就判为前景了  
152.                 }  
153.                   
154.                 dst[x] = (uchar)(-(kHit >= kForeground));//判决：(ucahr)(-true) = 255;(uchar)(-(false)) = 0;  
155.             }  
156.         }  
157.         else//如果学习速率小于等于0，则没有背景更新过程，其他过程类似  
158.         {  
159.             for( x = 0; x < cols; x++, mptr += K )  
160.             {  
161.                 float pix = src[x];  
162.                 int kHit = -1, kForeground = -1;  
163.                   
164.                 for( k = 0; k < K; k++ )  
165.                 {  
166.                     if( mptr[k].weight < FLT_EPSILON )  
167.                         break;  
168.                     float mu = mptr[k].mean;  
169.                     float var = mptr[k].var;  
170.                     float diff = pix - mu;  
171.                     float d2 = diff*diff;  
172.                     if( d2 < vT*var )  
173.                     {  
174.                         kHit = k;  
175.                         break;  
176.                     }  
177.                 }  
178.                   
179.                 if( kHit >= 0 )  
180.                 {  
181.                     float wsum = 0;  
182.                     for( k = 0; k < K; k++ )  
183.                     {  
184.                         wsum += mptr[k].weight;  
185.                         if( wsum > T )  
186.                         {  
187.                             kForeground = k+1;  
188.                             break;  
189.                         }  
190.                     }  
191.                 }  
192.                   
193.                 dst[x] = (uchar)(kHit < 0 || kHit >= kForeground ? 255 : 0);  
194.             }  
195.         }  
196.     }  
197. }  
198.   
199.       
200. static void process8uC3( const Mat& image, Mat& fgmask, double learningRate,  
201.                          Mat& bgmodel, int nmixtures, double backgroundRatio,  
202.                          double varThreshold, double noiseSigma )  
203. {  
204.     int x, y, k, k1, rows = image.rows, cols = image.cols;  
205.     float alpha = (float)learningRate, T = (float)backgroundRatio, vT = (float)varThreshold;  
206.     int K = nmixtures;  
207.       
208.     const float w0 = (float)defaultInitialWeight;  
209.     const float sk0 = (float)(w0/(defaultNoiseSigma*2*sqrt(3.)));  
210.     const float var0 = (float)(defaultNoiseSigma*defaultNoiseSigma*4);  
211.     const float minVar = (float)(noiseSigma*noiseSigma);  
212.     MixData<Vec3f>* mptr = (MixData<Vec3f>*)bgmodel.data;  
213.       
214.     for( y = 0; y < rows; y++ )  
215.     {  
216.         const uchar* src = image.ptr<uchar>(y);  
217.         uchar* dst = fgmask.ptr<uchar>(y);  
218.           
219.         if( alpha > 0 )  
220.         {  
221.             for( x = 0; x < cols; x++, mptr += K )  
222.             {  
223.                 float wsum = 0;  
224.                 Vec3f pix(src[x*3], src[x*3+1], src[x*3+2]);  
225.                 int kHit = -1, kForeground = -1;  
226.                   
227.                 for( k = 0; k < K; k++ )  
228.                 {  
229.                     float w = mptr[k].weight;  
230.                     wsum += w;  
231.                     if( w < FLT_EPSILON )  
232.                         break;  
233.                     Vec3f mu = mptr[k].mean;  
234.                     Vec3f var = mptr[k].var;  
235.                     Vec3f diff = pix - mu;  
236.                     float d2 = diff.dot(diff);  
237.                     if( d2 < vT*(var[0] + var[1] + var[2]) )  
238.                     {  
239.                         wsum -= w;  
240.                         float dw = alpha*(1.f - w);  
241.                         mptr[k].weight = w + dw;  
242.                         mptr[k].mean = mu + alpha*diff;  
243.                         var = Vec3f(max(var[0] + alpha*(diff[0]*diff[0] - var[0]), minVar),  
244.                                     max(var[1] + alpha*(diff[1]*diff[1] - var[1]), minVar),  
245.                                     max(var[2] + alpha*(diff[2]*diff[2] - var[2]), minVar));  
246.                         mptr[k].var = var;  
247.                         mptr[k].sortKey = w/sqrt(var[0] + var[1] + var[2]);  
248.                           
249.                         for( k1 = k-1; k1 >= 0; k1-- )  
250.                         {  
251.                             if( mptr[k1].sortKey >= mptr[k1+1].sortKey )  
252.                                 break;  
253.                             std::swap( mptr[k1], mptr[k1+1] );  
254.                         }  
255.                           
256.                         kHit = k1+1;  
257.                         break;  
258.                     }  
259.                 }  
260.                   
261.                 if( kHit < 0 ) // no appropriate gaussian mixture found at all, remove the weakest mixture and create a new one  
262.                 {  
263.                     kHit = k = min(k, K-1);  
264.                     wsum += w0 - mptr[k].weight;  
265.                     mptr[k].weight = w0;  
266.                     mptr[k].mean = pix;  
267.                     mptr[k].var = Vec3f(var0, var0, var0);  
268.                     mptr[k].sortKey = sk0;  
269.                 }  
270.                 else  
271.                     for( ; k < K; k++ )  
272.                         wsum += mptr[k].weight;  
273.               
274.                 float wscale = 1.f/wsum;  
275.                 wsum = 0;  
276.                 for( k = 0; k < K; k++ )  
277.                 {  
278.                     wsum += mptr[k].weight *= wscale;  
279.                     mptr[k].sortKey *= wscale;  
280.                     if( wsum > T && kForeground < 0 )  
281.                         kForeground = k+1;  
282.                 }  
283.                   
284.                 dst[x] = (uchar)(-(kHit >= kForeground));  
285.             }  
286.         }  
287.         else  
288.         {  
289.             for( x = 0; x < cols; x++, mptr += K )  
290.             {  
291.                 Vec3f pix(src[x*3], src[x*3+1], src[x*3+2]);  
292.                 int kHit = -1, kForeground = -1;  
293.                   
294.                 for( k = 0; k < K; k++ )  
295.                 {  
296.                     if( mptr[k].weight < FLT_EPSILON )  
297.                         break;  
298.                     Vec3f mu = mptr[k].mean;  
299.                     Vec3f var = mptr[k].var;  
300.                     Vec3f diff = pix - mu;  
301.                     float d2 = diff.dot(diff);  
302.                     if( d2 < vT*(var[0] + var[1] + var[2]) )  
303.                     {  
304.                         kHit = k;  
305.                         break;  
306.                     }  
307.                 }  
308.    
309.                 if( kHit >= 0 )  
310.                 {  
311.                     float wsum = 0;  
312.                     for( k = 0; k < K; k++ )  
313.                     {  
314.                         wsum += mptr[k].weight;  
315.                         if( wsum > T )  
316.                         {  
317.                             kForeground = k+1;  
318.                             break;  
319.                         }  
320.                     }  
321.                 }  
322.                   
323.                 dst[x] = (uchar)(kHit < 0 || kHit >= kForeground ? 255 : 0);  
324.             }  
325.         }  
326.     }  
327. }  
328.   
329. void BackgroundSubtractorMOG::operator()(InputArray _image, OutputArray _fgmask, double learningRate)  
330. {  
331.     Mat image = _image.getMat();  
332.     bool needToInitialize = nframes == 0 || learningRate >= 1 || image.size() != frameSize || image.type() != frameType;//是否需要初始化  
333.       
334.     if( needToInitialize )  
335.         initialize(image.size(), image.type());//初始化  
336.       
337.     CV_Assert( image.depth() == CV_8U );  
338.     _fgmask.create( image.size(), CV_8U );  
339.     Mat fgmask = _fgmask.getMat();  
340.       
341.     ++nframes;  
342.     learningRate = learningRate >= 0 && nframes > 1 ? learningRate : 1./min( nframes, history );  
343.     CV_Assert(learningRate >= 0);  
344.       
345.     if( image.type() == CV_8UC1 )//处理灰度图像  
346.         process8uC1( image, fgmask, learningRate, bgmodel, nmixtures, backgroundRatio, varThreshold, noiseSigma );  
347.     else if( image.type() == CV_8UC3 )//处理彩色图像  
348.         process8uC3( image, fgmask, learningRate, bgmodel, nmixtures, backgroundRatio, varThreshold, noiseSigma );  
349.     else  
350.         CV_Error( CV_StsUnsupportedFormat, "Only 1- and 3-channel 8-bit images are supported in BackgroundSubtractorMOG" );  
351. }  
352.       
353. }  

static const int defaultNMixtures = 5;//默认混合模型个数
static const int defaultHistory = 200;//默认历史帧数
static const double defaultBackgroundRatio = 0.7;//默认背景门限
static const double defaultVarThreshold = 2.5*2.5;//默认方差门限
static const double defaultNoiseSigma = 30*0.5;//默认噪声方差
static const double defaultInitialWeight = 0.05;//默认初始权值
 //不带参数的构造函数，使用默认值  
BackgroundSubtractorMOG::BackgroundSubtractorMOG()
{
    frameSize = Size(0,0);
    frameType = 0;
    
    nframes = 0;
    nmixtures = defaultNMixtures;
    history = defaultHistory;
    varThreshold = defaultVarThreshold;
    backgroundRatio = defaultBackgroundRatio;
    noiseSigma = defaultNoiseSigma;
}
//带参数的构造函数，使用参数传进来的值    
BackgroundSubtractorMOG::BackgroundSubtractorMOG(int _history, int _nmixtures,
                                                 double _backgroundRatio,
                                                 double _noiseSigma)
{
    frameSize = Size(0,0);
    frameType = 0;
    
    nframes = 0;
    nmixtures = min(_nmixtures > 0 ? _nmixtures : defaultNMixtures, 8);//不能超过8个，否则就用默认的
    history = _history > 0 ? _history : defaultHistory;//不能小于0，否则就用默认的
    varThreshold = defaultVarThreshold;//门限使用默认的
    backgroundRatio = min(_backgroundRatio > 0 ? _backgroundRatio : 0.95, 1.);//背景门限必须大于0，小于1，否则使用0.95
    noiseSigma = _noiseSigma <= 0 ? defaultNoiseSigma : _noiseSigma;//噪声门限大于0
}
    
BackgroundSubtractorMOG::~BackgroundSubtractorMOG()
{
}


void BackgroundSubtractorMOG::initialize(Size _frameSize, int _frameType)
{
    frameSize = _frameSize;
    frameType = _frameType;
    nframes = 0;
    
    int nchannels = CV_MAT_CN(frameType);
    CV_Assert( CV_MAT_DEPTH(frameType) == CV_8U );
    
    // for each gaussian mixture of each pixel bg model we store ...
    // the mixture sort key (w/sum_of_variances), the mixture weight (w),
    // the mean (nchannels values) and
    // the diagonal covariance matrix (another nchannels values)
    bgmodel.create( 1, frameSize.height*frameSize.width*nmixtures*(2 + 2*nchannels), CV_32F );//初始化一个1行*XX列的矩阵
	//XX是这样计算的：图像的行*列*混合模型的个数*（1（优先级） + 1（权值） + 2（均值 + 方差） * 通道数）
    bgmodel = Scalar::all(0);//清零
}

//设为模版，就是考虑到了彩色图像与灰度图像两种情况    
template<typename VT> struct MixData
{
    float sortKey;
    float weight;
    VT mean;
    VT var;
};

    
static void process8uC1( const Mat& image, Mat& fgmask, double learningRate,
                         Mat& bgmodel, int nmixtures, double backgroundRatio,
                         double varThreshold, double noiseSigma )
{
    int x, y, k, k1, rows = image.rows, cols = image.cols;
    float alpha = (float)learningRate, T = (float)backgroundRatio, vT = (float)varThreshold;//学习速率、背景门限、方差门限
    int K = nmixtures;//混合模型个数
    MixData<float>* mptr = (MixData<float>*)bgmodel.data;
    
    const float w0 = (float)defaultInitialWeight;//初始权值
    const float sk0 = (float)(w0/(defaultNoiseSigma*2));//初始优先级
    const float var0 = (float)(defaultNoiseSigma*defaultNoiseSigma*4);//初始方差
    const float minVar = (float)(noiseSigma*noiseSigma);//最小方差
    
    for( y = 0; y < rows; y++ )
    {
        const uchar* src = image.ptr<uchar>(y);
        uchar* dst = fgmask.ptr<uchar>(y);
        
        if( alpha > 0 )//如果学习速率为0，则退化为背景相减
        {
            for( x = 0; x < cols; x++, mptr += K )
            {
                float wsum = 0;
                float pix = src[x];//每个像素
                int kHit = -1, kForeground = -1;//是否属于模型，是否属于前景
                
                for( k = 0; k < K; k++ )//每个高斯模型
                {
                    float w = mptr[k].weight;//当前模型的权值
                    wsum += w;//权值累加
                    if( w < FLT_EPSILON )
                        break;
                    float mu = mptr[k].mean;//当前模型的均值
                    float var = mptr[k].var;//当前模型的方差
                    float diff = pix - mu;//当前像素与模型均值之差
                    float d2 = diff*diff;//平方
                    if( d2 < vT*var )//是否小于方门限，即是否属于该模型
                    {
                        wsum -= w;//如果匹配，则把它减去，因为之后会更新它
                        float dw = alpha*(1.f - w);
                        mptr[k].weight = w + dw;//增加权值
						//注意源文章中涉及概率的部分多进行了简化，将概率变为1
                        mptr[k].mean = mu + alpha*diff;//修正均值		
                        var = max(var + alpha*(d2 - var), minVar);//开始时方差清零0，所以这里使用噪声方差作为默认方差，否则使用上一次方差
                        mptr[k].var = var;//修正方差
                        mptr[k].sortKey = w/sqrt(var);//重新计算优先级，貌似这里不对，应该使用更新后的mptr[k].weight而不是w
                        
                        for( k1 = k-1; k1 >= 0; k1-- )//从匹配的第k个模型开始向前比较，如果更新后的单高斯模型优先级超过他前面的那个，则交换顺序
                        {
                            if( mptr[k1].sortKey >= mptr[k1+1].sortKey )//如果优先级没有发生改变，则停止比较
                                break;
                            std::swap( mptr[k1], mptr[k1+1] );//交换它们的顺序，始终保证优先级最大的在前面
                        }
                        
                        kHit = k1+1;//记录属于哪个模型
                        break;
                    }
                }
                
                if( kHit < 0 ) // no appropriate gaussian mixture found at all, remove the weakest mixture and create a new one
								//当前像素不属于任何一个模型
                {
					//初始化一个新模型
                    kHit = k = min(k, K-1);//有两种情况，当最开始的初始化时，k并不是等于K-1的
                    wsum += w0 - mptr[k].weight;//从权值总和中减去原来的那个模型，并加上默认权值
                    mptr[k].weight = w0;//初始化权值
                    mptr[k].mean = pix;//初始化均值
                    mptr[k].var = var0;	//初始化方差
                    mptr[k].sortKey = sk0;//初始化权值
                }
                else
                    for( ; k < K; k++ )
                        wsum += mptr[k].weight;//求出剩下的总权值
                
                float wscale = 1.f/wsum;//归一化
                wsum = 0;
                for( k = 0; k < K; k++ )
                {
                    wsum += mptr[k].weight *= wscale;
                    mptr[k].sortKey *= wscale;//计算归一化权值
                    if( wsum > T && kForeground < 0 )
                        kForeground = k+1;//第几个模型之后就判为前景了
                }
                
                dst[x] = (uchar)(-(kHit >= kForeground));//判决：(ucahr)(-true) = 255;(uchar)(-(false)) = 0;
            }
        }
        else//如果学习速率小于等于0，则没有背景更新过程，其他过程类似
        {
            for( x = 0; x < cols; x++, mptr += K )
            {
                float pix = src[x];
                int kHit = -1, kForeground = -1;
                
                for( k = 0; k < K; k++ )
                {
                    if( mptr[k].weight < FLT_EPSILON )
                        break;
                    float mu = mptr[k].mean;
                    float var = mptr[k].var;
                    float diff = pix - mu;
                    float d2 = diff*diff;
                    if( d2 < vT*var )
                    {
                        kHit = k;
                        break;
                    }
                }
                
                if( kHit >= 0 )
                {
                    float wsum = 0;
                    for( k = 0; k < K; k++ )
                    {
                        wsum += mptr[k].weight;
                        if( wsum > T )
                        {
                            kForeground = k+1;
                            break;
                        }
                    }
                }
                
                dst[x] = (uchar)(kHit < 0 || kHit >= kForeground ? 255 : 0);
            }
        }
    }
}

    
static void process8uC3( const Mat& image, Mat& fgmask, double learningRate,
                         Mat& bgmodel, int nmixtures, double backgroundRatio,
                         double varThreshold, double noiseSigma )
{
    int x, y, k, k1, rows = image.rows, cols = image.cols;
    float alpha = (float)learningRate, T = (float)backgroundRatio, vT = (float)varThreshold;
    int K = nmixtures;
    
    const float w0 = (float)defaultInitialWeight;
    const float sk0 = (float)(w0/(defaultNoiseSigma*2*sqrt(3.)));
    const float var0 = (float)(defaultNoiseSigma*defaultNoiseSigma*4);
    const float minVar = (float)(noiseSigma*noiseSigma);
    MixData<Vec3f>* mptr = (MixData<Vec3f>*)bgmodel.data;
    
    for( y = 0; y < rows; y++ )
    {
        const uchar* src = image.ptr<uchar>(y);
        uchar* dst = fgmask.ptr<uchar>(y);
        
        if( alpha > 0 )
        {
            for( x = 0; x < cols; x++, mptr += K )
            {
                float wsum = 0;
                Vec3f pix(src[x*3], src[x*3+1], src[x*3+2]);
                int kHit = -1, kForeground = -1;
                
                for( k = 0; k < K; k++ )
                {
                    float w = mptr[k].weight;
                    wsum += w;
                    if( w < FLT_EPSILON )
                        break;
                    Vec3f mu = mptr[k].mean;
                    Vec3f var = mptr[k].var;
                    Vec3f diff = pix - mu;
                    float d2 = diff.dot(diff);
                    if( d2 < vT*(var[0] + var[1] + var[2]) )
                    {
                        wsum -= w;
                        float dw = alpha*(1.f - w);
                        mptr[k].weight = w + dw;
                        mptr[k].mean = mu + alpha*diff;
                        var = Vec3f(max(var[0] + alpha*(diff[0]*diff[0] - var[0]), minVar),
                                    max(var[1] + alpha*(diff[1]*diff[1] - var[1]), minVar),
                                    max(var[2] + alpha*(diff[2]*diff[2] - var[2]), minVar));
                        mptr[k].var = var;
                        mptr[k].sortKey = w/sqrt(var[0] + var[1] + var[2]);
                        
                        for( k1 = k-1; k1 >= 0; k1-- )
                        {
                            if( mptr[k1].sortKey >= mptr[k1+1].sortKey )
                                break;
                            std::swap( mptr[k1], mptr[k1+1] );
                        }
                        
                        kHit = k1+1;
                        break;
                    }
                }
                
                if( kHit < 0 ) // no appropriate gaussian mixture found at all, remove the weakest mixture and create a new one
                {
                    kHit = k = min(k, K-1);
                    wsum += w0 - mptr[k].weight;
                    mptr[k].weight = w0;
                    mptr[k].mean = pix;
                    mptr[k].var = Vec3f(var0, var0, var0);
                    mptr[k].sortKey = sk0;
                }
                else
                    for( ; k < K; k++ )
                        wsum += mptr[k].weight;
            
                float wscale = 1.f/wsum;
                wsum = 0;
                for( k = 0; k < K; k++ )
                {
                    wsum += mptr[k].weight *= wscale;
                    mptr[k].sortKey *= wscale;
                    if( wsum > T && kForeground < 0 )
                        kForeground = k+1;
                }
                
                dst[x] = (uchar)(-(kHit >= kForeground));
            }
        }
        else
        {
            for( x = 0; x < cols; x++, mptr += K )
            {
                Vec3f pix(src[x*3], src[x*3+1], src[x*3+2]);
                int kHit = -1, kForeground = -1;
                
                for( k = 0; k < K; k++ )
                {
                    if( mptr[k].weight < FLT_EPSILON )
                        break;
                    Vec3f mu = mptr[k].mean;
                    Vec3f var = mptr[k].var;
                    Vec3f diff = pix - mu;
                    float d2 = diff.dot(diff);
                    if( d2 < vT*(var[0] + var[1] + var[2]) )
                    {
                        kHit = k;
                        break;
                    }
                }
 
                if( kHit >= 0 )
                {
                    float wsum = 0;
                    for( k = 0; k < K; k++ )
                    {
                        wsum += mptr[k].weight;
                        if( wsum > T )
                        {
                            kForeground = k+1;
                            break;
                        }
                    }
                }
                
                dst[x] = (uchar)(kHit < 0 || kHit >= kForeground ? 255 : 0);
            }
        }
    }
}

void BackgroundSubtractorMOG::operator()(InputArray _image, OutputArray _fgmask, double learningRate)
{
    Mat image = _image.getMat();
    bool needToInitialize = nframes == 0 || learningRate >= 1 || image.size() != frameSize || image.type() != frameType;//是否需要初始化
    
    if( needToInitialize )
        initialize(image.size(), image.type());//初始化
    
    CV_Assert( image.depth() == CV_8U );
    _fgmask.create( image.size(), CV_8U );
    Mat fgmask = _fgmask.getMat();
    
    ++nframes;
    learningRate = learningRate >= 0 && nframes > 1 ? learningRate : 1./min( nframes, history );
    CV_Assert(learningRate >= 0);
    
    if( image.type() == CV_8UC1 )//处理灰度图像
        process8uC1( image, fgmask, learningRate, bgmodel, nmixtures, backgroundRatio, varThreshold, noiseSigma );
    else if( image.type() == CV_8UC3 )//处理彩色图像
        process8uC3( image, fgmask, learningRate, bgmodel, nmixtures, backgroundRatio, varThreshold, noiseSigma );
    else
        CV_Error( CV_StsUnsupportedFormat, "Only 1- and 3-channel 8-bit images are supported in BackgroundSubtractorMOG" );
}
    
}

其中处理3通道彩色图像与处理单通道灰度图像类似，我就没有进行注释了。

其中有几点需要注意：

1.在高斯混合模型中需要使用概率更新参数的地方，程序中都简化成为了1处理，否则计算一个正态分布的概率还是挺花时间的。（程序作者在注释中也指出了他不是完全按照论文写成的，而是做了一些小的修改）。但是除了将概率换成1，其他地方还是严格按照公式的，大家可以仔细推导一下，就会看出其中的差异。

2.作者原文中是如果没有一个高斯模型与该像素点匹配，则去掉一个一个概率最小的，而用当前像素初始化的分布来替代他，而在这里变成了去掉优先级最小的。

3.程序中为了避免多次做循环，把一些步骤拆开做了，比如归一化权值需要先求出总权值，调整权值后的排序之类的，计算背景模型个数等等。减少了遍历的次数。其中的巧妙之处也不得不佩服作者的良苦用心。

3.就是似乎更新优先级的计算有点小问题，也可能是我理解不对。

4.在初始化时，可以使用多种方式，大家一看程序就明白了。

最后附上一个小的示例程序，教你如何使用高斯混合模型：

[cpp] view plain copy

print ?

1. int main()  
2. {  
3.     VideoCapture capture("D:/videos/shadow/use3.MPG");  
4.     if( !capture.isOpened() )  
5.     {  
6.         cout<<"读取视频失败"<<endl;  
7.         return -1;  
8.     }  
9.     //获取整个帧数  
10.     long totalFrameNumber = capture.get(CV_CAP_PROP_FRAME_COUNT);  
11.     cout<<"整个视频共"<<totalFrameNumber<<"帧"<<endl;  
12.   
13.     //设置开始帧()  
14.     long frameToStart = 200;  
15.     capture.set( CV_CAP_PROP_POS_FRAMES,frameToStart);  
16.     cout<<"从第"<<frameToStart<<"帧开始读"<<endl;  
17.   
18.     //设置结束帧  
19.     int frameToStop = 650;  
20.   
21.     if(frameToStop < frameToStart)  
22.     {  
23.         cout<<"结束帧小于开始帧，程序错误，即将退出！"<<endl;  
24.         return -1;  
25.     }  
26.     else  
27.     {  
28.         cout<<"结束帧为：第"<<frameToStop<<"帧"<<endl;  
29.     }  
30.   
31.     double rate = capture.get(CV_CAP_PROP_FPS);  
32.     int delay = 1000/rate;  
33.   
34.     Mat frame;  
35.     //前景图片  
36.     Mat foreground;  
37.   
38.   
39.     //使用默认参数调用混合高斯模型  
40.     BackgroundSubtractorMOG mog;  
41.     bool stop(false);  
42.     //currentFrame是在循环体中控制读取到指定的帧后循环结束的变量  
43.     long currentFrame = frameToStart;  
44.     while( !stop )  
45.     {  
46.         if( !capture.read(frame) )  
47.         {  
48.             cout<<"从视频中读取图像失败或者读完整个视频"<<endl;  
49.             return -2;  
50.         }  
51.         cvtColor(frame,frame,CV_RGB2GRAY);  
52.         imshow("输入视频",frame);  
53.         //参数为：输入图像、输出图像、学习速率  
54.         mog(frame,foreground,0.01);  
55.   
56.   
57.         imshow("前景",foreground);  
58.   
59.         //按ESC键退出，按其他键会停止在当前帧  
60.   
61.         int c = waitKey(delay);  
62.   
63.         if ( (char)c == 27 || currentFrame >= frameToStop)  
64.         {  
65.             stop = true;  
66.         }  
67.         if ( c >= 0)  
68.         {  
69.             waitKey(0);  
70.         }  
71.         currentFrame++;  
72.   
73.     }  
74.   
75.     waitKey(0);  
76. }  

int main()
{
	VideoCapture capture("D:/videos/shadow/use3.MPG");
	if( !capture.isOpened() )
	{
		cout<<"读取视频失败"<<endl;
		return -1;
	}
	//获取整个帧数
	long totalFrameNumber = capture.get(CV_CAP_PROP_FRAME_COUNT);
	cout<<"整个视频共"<<totalFrameNumber<<"帧"<<endl;

	//设置开始帧()
	long frameToStart = 200;
	capture.set( CV_CAP_PROP_POS_FRAMES,frameToStart);
	cout<<"从第"<<frameToStart<<"帧开始读"<<endl;

	//设置结束帧
	int frameToStop = 650;

	if(frameToStop < frameToStart)
	{
		cout<<"结束帧小于开始帧，程序错误，即将退出！"<<endl;
		return -1;
	}
	else
	{
		cout<<"结束帧为：第"<<frameToStop<<"帧"<<endl;
	}

	double rate = capture.get(CV_CAP_PROP_FPS);
	int delay = 1000/rate;

	Mat frame;
	//前景图片
	Mat foreground;


	//使用默认参数调用混合高斯模型
	BackgroundSubtractorMOG mog;
	bool stop(false);
	//currentFrame是在循环体中控制读取到指定的帧后循环结束的变量
	long currentFrame = frameToStart;
	while( !stop )
	{
		if( !capture.read(frame) )
		{
			cout<<"从视频中读取图像失败或者读完整个视频"<<endl;
			return -2;
		}
		cvtColor(frame,frame,CV_RGB2GRAY);
		imshow("输入视频",frame);
		//参数为：输入图像、输出图像、学习速率
		mog(frame,foreground,0.01);


		imshow("前景",foreground);

		//按ESC键退出，按其他键会停止在当前帧

		int c = waitKey(delay);

		if ( (char)c == 27 || currentFrame >= frameToStop)
		{
			stop = true;
		}
		if ( c >= 0)
		{
			waitKey(0);
		}
		currentFrame++;

	}

	waitKey(0);
}

就说这么多吧，虽然我

顶

17

踩

1

 

 

• 上一篇我的OpenCV学习笔记（24）：详细讨论OpenCV中的数据结构
• 下一篇请不要责备那些翻译书的人

  相关文章推荐

• • VC++深入详解（1）：MFC框架程序剖析
• • 【MFC】程序框架及基础知识
• • 【OpenCV】高斯混合背景建模
• • opencv BackgroundSubtractorMOG和BackgroundSubtractorMOG2的区别
• • opencv中BackgroundSubtractorMOG类中的背景提取

• • opencv中BackgroundSubtractorMOG问题及解决方法
• • BackgroundSubtractorMOG2前后背景分离
• • 背景提取—修改高斯混合模型BackgroundSubtractorMOG2中的参数及使用
• • opencv BackgroundSubtractorMOG2重要方法
• • 【OpenCV】高斯混合背景提取