• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 //  By downloading, copying, installing or using the software you agree to this license.
6 //  If you do not agree to this license, do not download, install,
7 //  copy or use the software.
8 //
9 //
10 //                          License Agreement
11 //                For Open Source Computer Vision Library
12 //
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Copyright (C) 2013, OpenCV Foundation, all rights reserved.
15 // Third party copyrights are property of their respective owners.
16 //
17 // Redistribution and use in source and binary forms, with or without modification,
18 // are permitted provided that the following conditions are met:
19 //
20 //   * Redistribution's of source code must retain the above copyright notice,
21 //     this list of conditions and the following disclaimer.
22 //
23 //   * Redistribution's in binary form must reproduce the above copyright notice,
24 //     this list of conditions and the following disclaimer in the documentation
25 //     and/or other materials provided with the distribution.
26 //
27 //   * The name of the copyright holders may not be used to endorse or promote products
28 //     derived from this software without specific prior written permission.
29 //
30 // This software is provided by the copyright holders and contributors "as is" and
31 // any express or implied warranties, including, but not limited to, the implied
32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
33 // In no event shall the Intel Corporation or contributors be liable for any direct,
34 // indirect, incidental, special, exemplary, or consequential damages
35 // (including, but not limited to, procurement of substitute goods or services;
36 // loss of use, data, or profits; or business interruption) however caused
37 // and on any theory of liability, whether in contract, strict liability,
38 // or tort (including negligence or otherwise) arising in any way out of
39 // the use of this software, even if advised of the possibility of such damage.
40 //
41 //M*/
42 
43 /*//Implementation of the Gaussian mixture model background subtraction from:
44 //
45 //"Improved adaptive Gausian mixture model for background subtraction"
46 //Z.Zivkovic
47 //International Conference Pattern Recognition, UK, August, 2004
48 //http://www.zoranz.net/Publications/zivkovic2004ICPR.pdf
49 //The code is very fast and performs also shadow detection.
50 //Number of Gausssian components is adapted per pixel.
51 //
52 // and
53 //
54 //"Efficient Adaptive Density Estimapion per Image Pixel for the Task of Background Subtraction"
55 //Z.Zivkovic, F. van der Heijden
56 //Pattern Recognition Letters, vol. 27, no. 7, pages 773-780, 2006.
57 //
58 //The algorithm similar to the standard Stauffer&Grimson algorithm with
59 //additional selection of the number of the Gaussian components based on:
60 //
61 //"Recursive unsupervised learning of finite mixture models "
62 //Z.Zivkovic, F.van der Heijden
63 //IEEE Trans. on Pattern Analysis and Machine Intelligence, vol.26, no.5, pages 651-656, 2004
64 //http://www.zoranz.net/Publications/zivkovic2004PAMI.pdf
65 //
66 //
67 //Example usage with as cpp class
68 // BackgroundSubtractorMOG2 bg_model;
69 //For each new image the model is updates using:
70 // bg_model(img, fgmask);
71 //
72 //Example usage as part of the CvBGStatModel:
73 // CvBGStatModel* bg_model = cvCreateGaussianBGModel2( first_frame );
74 //
75 // //update for each frame
76 // cvUpdateBGStatModel( tmp_frame, bg_model );//segmentation result is in bg_model->foreground
77 //
78 // //release at the program termination
79 // cvReleaseBGStatModel( &bg_model );
80 //
81 //Author: Z.Zivkovic, www.zoranz.net
82 //Date: 7-April-2011, Version:1.0
83 ///////////*/
84 
85 #include "precomp.hpp"
86 #include "opencl_kernels_video.hpp"
87 
88 namespace cv
89 {
90 
91 /*
92  Interface of Gaussian mixture algorithm from:
93 
94  "Improved adaptive Gausian mixture model for background subtraction"
95  Z.Zivkovic
96  International Conference Pattern Recognition, UK, August, 2004
97  http://www.zoranz.net/Publications/zivkovic2004ICPR.pdf
98 
99  Advantages:
100  -fast - number of Gausssian components is constantly adapted per pixel.
101  -performs also shadow detection (see bgfg_segm_test.cpp example)
102 
103 */
104 
105 // default parameters of gaussian background detection algorithm
106 static const int defaultHistory2 = 500; // Learning rate; alpha = 1/defaultHistory2
107 static const float defaultVarThreshold2 = 4.0f*4.0f;
108 static const int defaultNMixtures2 = 5; // maximal number of Gaussians in mixture
109 static const float defaultBackgroundRatio2 = 0.9f; // threshold sum of weights for background test
110 static const float defaultVarThresholdGen2 = 3.0f*3.0f;
111 static const float defaultVarInit2 = 15.0f; // initial variance for new components
112 static const float defaultVarMax2 = 5*defaultVarInit2;
113 static const float defaultVarMin2 = 4.0f;
114 
115 // additional parameters
116 static const float defaultfCT2 = 0.05f; // complexity reduction prior constant 0 - no reduction of number of components
117 static const unsigned char defaultnShadowDetection2 = (unsigned char)127; // value to use in the segmentation mask for shadows, set 0 not to do shadow detection
118 static const float defaultfTau = 0.5f; // Tau - shadow threshold, see the paper for explanation
119 
120 
121 class BackgroundSubtractorMOG2Impl : public BackgroundSubtractorMOG2
122 {
123 public:
124     //! the default constructor
BackgroundSubtractorMOG2Impl()125     BackgroundSubtractorMOG2Impl()
126     {
127         frameSize = Size(0,0);
128         frameType = 0;
129 
130         nframes = 0;
131         history = defaultHistory2;
132         varThreshold = defaultVarThreshold2;
133         bShadowDetection = 1;
134 
135         nmixtures = defaultNMixtures2;
136         backgroundRatio = defaultBackgroundRatio2;
137         fVarInit = defaultVarInit2;
138         fVarMax  = defaultVarMax2;
139         fVarMin = defaultVarMin2;
140 
141         varThresholdGen = defaultVarThresholdGen2;
142         fCT = defaultfCT2;
143         nShadowDetection =  defaultnShadowDetection2;
144         fTau = defaultfTau;
145 
146         opencl_ON = true;
147     }
148     //! the full constructor that takes the length of the history,
149     // the number of gaussian mixtures, the background ratio parameter and the noise strength
BackgroundSubtractorMOG2Impl(int _history,float _varThreshold,bool _bShadowDetection=true)150     BackgroundSubtractorMOG2Impl(int _history,  float _varThreshold, bool _bShadowDetection=true)
151     {
152         frameSize = Size(0,0);
153         frameType = 0;
154 
155         nframes = 0;
156         history = _history > 0 ? _history : defaultHistory2;
157         varThreshold = (_varThreshold>0)? _varThreshold : defaultVarThreshold2;
158         bShadowDetection = _bShadowDetection;
159 
160         nmixtures = defaultNMixtures2;
161         backgroundRatio = defaultBackgroundRatio2;
162         fVarInit = defaultVarInit2;
163         fVarMax  = defaultVarMax2;
164         fVarMin = defaultVarMin2;
165 
166         varThresholdGen = defaultVarThresholdGen2;
167         fCT = defaultfCT2;
168         nShadowDetection =  defaultnShadowDetection2;
169         fTau = defaultfTau;
170         name_ = "BackgroundSubtractor.MOG2";
171 
172         opencl_ON = true;
173     }
174     //! the destructor
~BackgroundSubtractorMOG2Impl()175     ~BackgroundSubtractorMOG2Impl() {}
176     //! the update operator
177     void apply(InputArray image, OutputArray fgmask, double learningRate=-1);
178 
179     //! computes a background image which are the mean of all background gaussians
180     virtual void getBackgroundImage(OutputArray backgroundImage) const;
181 
182     //! re-initiaization method
initialize(Size _frameSize,int _frameType)183     void initialize(Size _frameSize, int _frameType)
184     {
185         frameSize = _frameSize;
186         frameType = _frameType;
187         nframes = 0;
188 
189         int nchannels = CV_MAT_CN(frameType);
190         CV_Assert( nchannels <= CV_CN_MAX );
191         CV_Assert( nmixtures <= 255);
192 
193         if (ocl::useOpenCL() && opencl_ON)
194         {
195             create_ocl_apply_kernel();
196             kernel_getBg.create("getBackgroundImage2_kernel", ocl::video::bgfg_mog2_oclsrc, format( "-D CN=%d -D NMIXTURES=%d", nchannels, nmixtures));
197 
198             if (kernel_apply.empty() || kernel_getBg.empty())
199                 opencl_ON = false;
200         }
201         else opencl_ON = false;
202 
203         if (opencl_ON)
204         {
205             u_weight.create(frameSize.height * nmixtures, frameSize.width, CV_32FC1);
206             u_weight.setTo(Scalar::all(0));
207 
208             u_variance.create(frameSize.height * nmixtures, frameSize.width, CV_32FC1);
209             u_variance.setTo(Scalar::all(0));
210 
211             if (nchannels==3)
212                 nchannels=4;
213             u_mean.create(frameSize.height * nmixtures, frameSize.width, CV_32FC(nchannels)); //4 channels
214             u_mean.setTo(Scalar::all(0));
215 
216             //make the array for keeping track of the used modes per pixel - all zeros at start
217             u_bgmodelUsedModes.create(frameSize, CV_8UC1);
218             u_bgmodelUsedModes.setTo(cv::Scalar::all(0));
219         }
220         else
221         {
222             // for each gaussian mixture of each pixel bg model we store ...
223             // the mixture weight (w),
224             // the mean (nchannels values) and
225             // the covariance
226             bgmodel.create( 1, frameSize.height*frameSize.width*nmixtures*(2 + nchannels), CV_32F );
227             //make the array for keeping track of the used modes per pixel - all zeros at start
228             bgmodelUsedModes.create(frameSize,CV_8U);
229             bgmodelUsedModes = Scalar::all(0);
230         }
231     }
232 
getHistory() const233     virtual int getHistory() const { return history; }
setHistory(int _nframes)234     virtual void setHistory(int _nframes) { history = _nframes; }
235 
getNMixtures() const236     virtual int getNMixtures() const { return nmixtures; }
setNMixtures(int nmix)237     virtual void setNMixtures(int nmix) { nmixtures = nmix; }
238 
getBackgroundRatio() const239     virtual double getBackgroundRatio() const { return backgroundRatio; }
setBackgroundRatio(double _backgroundRatio)240     virtual void setBackgroundRatio(double _backgroundRatio) { backgroundRatio = (float)_backgroundRatio; }
241 
getVarThreshold() const242     virtual double getVarThreshold() const { return varThreshold; }
setVarThreshold(double _varThreshold)243     virtual void setVarThreshold(double _varThreshold) { varThreshold = _varThreshold; }
244 
getVarThresholdGen() const245     virtual double getVarThresholdGen() const { return varThresholdGen; }
setVarThresholdGen(double _varThresholdGen)246     virtual void setVarThresholdGen(double _varThresholdGen) { varThresholdGen = (float)_varThresholdGen; }
247 
getVarInit() const248     virtual double getVarInit() const { return fVarInit; }
setVarInit(double varInit)249     virtual void setVarInit(double varInit) { fVarInit = (float)varInit; }
250 
getVarMin() const251     virtual double getVarMin() const { return fVarMin; }
setVarMin(double varMin)252     virtual void setVarMin(double varMin) { fVarMin = (float)varMin; }
253 
getVarMax() const254     virtual double getVarMax() const { return fVarMax; }
setVarMax(double varMax)255     virtual void setVarMax(double varMax) { fVarMax = (float)varMax; }
256 
getComplexityReductionThreshold() const257     virtual double getComplexityReductionThreshold() const { return fCT; }
setComplexityReductionThreshold(double ct)258     virtual void setComplexityReductionThreshold(double ct) { fCT = (float)ct; }
259 
getDetectShadows() const260     virtual bool getDetectShadows() const { return bShadowDetection; }
setDetectShadows(bool detectshadows)261     virtual void setDetectShadows(bool detectshadows)
262     {
263         if ((bShadowDetection && detectshadows) || (!bShadowDetection && !detectshadows))
264             return;
265         bShadowDetection = detectshadows;
266         if (!kernel_apply.empty())
267         {
268             create_ocl_apply_kernel();
269             CV_Assert( !kernel_apply.empty() );
270         }
271     }
272 
getShadowValue() const273     virtual int getShadowValue() const { return nShadowDetection; }
setShadowValue(int value)274     virtual void setShadowValue(int value) { nShadowDetection = (uchar)value; }
275 
getShadowThreshold() const276     virtual double getShadowThreshold() const { return fTau; }
setShadowThreshold(double value)277     virtual void setShadowThreshold(double value) { fTau = (float)value; }
278 
write(FileStorage & fs) const279     virtual void write(FileStorage& fs) const
280     {
281         fs << "name" << name_
282         << "history" << history
283         << "nmixtures" << nmixtures
284         << "backgroundRatio" << backgroundRatio
285         << "varThreshold" << varThreshold
286         << "varThresholdGen" << varThresholdGen
287         << "varInit" << fVarInit
288         << "varMin" << fVarMin
289         << "varMax" << fVarMax
290         << "complexityReductionThreshold" << fCT
291         << "detectShadows" << (int)bShadowDetection
292         << "shadowValue" << (int)nShadowDetection
293         << "shadowThreshold" << fTau;
294     }
295 
read(const FileNode & fn)296     virtual void read(const FileNode& fn)
297     {
298         CV_Assert( (String)fn["name"] == name_ );
299         history = (int)fn["history"];
300         nmixtures = (int)fn["nmixtures"];
301         backgroundRatio = (float)fn["backgroundRatio"];
302         varThreshold = (double)fn["varThreshold"];
303         varThresholdGen = (float)fn["varThresholdGen"];
304         fVarInit = (float)fn["varInit"];
305         fVarMin = (float)fn["varMin"];
306         fVarMax = (float)fn["varMax"];
307         fCT = (float)fn["complexityReductionThreshold"];
308         bShadowDetection = (int)fn["detectShadows"] != 0;
309         nShadowDetection = saturate_cast<uchar>((int)fn["shadowValue"]);
310         fTau = (float)fn["shadowThreshold"];
311     }
312 
313 protected:
314     Size frameSize;
315     int frameType;
316     Mat bgmodel;
317     Mat bgmodelUsedModes;//keep track of number of modes per pixel
318 
319     //for OCL
320 
321     mutable bool opencl_ON;
322 
323     UMat u_weight;
324     UMat u_variance;
325     UMat u_mean;
326     UMat u_bgmodelUsedModes;
327 
328     mutable ocl::Kernel kernel_apply;
329     mutable ocl::Kernel kernel_getBg;
330 
331     int nframes;
332     int history;
333     int nmixtures;
334     //! here it is the maximum allowed number of mixture components.
335     //! Actual number is determined dynamically per pixel
336     double varThreshold;
337     // threshold on the squared Mahalanobis distance to decide if it is well described
338     // by the background model or not. Related to Cthr from the paper.
339     // This does not influence the update of the background. A typical value could be 4 sigma
340     // and that is varThreshold=4*4=16; Corresponds to Tb in the paper.
341 
342     /////////////////////////
343     // less important parameters - things you might change but be carefull
344     ////////////////////////
345     float backgroundRatio;
346     // corresponds to fTB=1-cf from the paper
347     // TB - threshold when the component becomes significant enough to be included into
348     // the background model. It is the TB=1-cf from the paper. So I use cf=0.1 => TB=0.
349     // For alpha=0.001 it means that the mode should exist for approximately 105 frames before
350     // it is considered foreground
351     // float noiseSigma;
352     float varThresholdGen;
353     //correspondts to Tg - threshold on the squared Mahalan. dist. to decide
354     //when a sample is close to the existing components. If it is not close
355     //to any a new component will be generated. I use 3 sigma => Tg=3*3=9.
356     //Smaller Tg leads to more generated components and higher Tg might make
357     //lead to small number of components but they can grow too large
358     float fVarInit;
359     float fVarMin;
360     float fVarMax;
361     //initial variance  for the newly generated components.
362     //It will will influence the speed of adaptation. A good guess should be made.
363     //A simple way is to estimate the typical standard deviation from the images.
364     //I used here 10 as a reasonable value
365     // min and max can be used to further control the variance
366     float fCT;//CT - complexity reduction prior
367     //this is related to the number of samples needed to accept that a component
368     //actually exists. We use CT=0.05 of all the samples. By setting CT=0 you get
369     //the standard Stauffer&Grimson algorithm (maybe not exact but very similar)
370 
371     //shadow detection parameters
372     bool bShadowDetection;//default 1 - do shadow detection
373     unsigned char nShadowDetection;//do shadow detection - insert this value as the detection result - 127 default value
374     float fTau;
375     // Tau - shadow threshold. The shadow is detected if the pixel is darker
376     //version of the background. Tau is a threshold on how much darker the shadow can be.
377     //Tau= 0.5 means that if pixel is more than 2 times darker then it is not shadow
378     //See: Prati,Mikic,Trivedi,Cucchiarra,"Detecting Moving Shadows...",IEEE PAMI,2003.
379 
380     String name_;
381 
382     bool ocl_getBackgroundImage(OutputArray backgroundImage) const;
383     bool ocl_apply(InputArray _image, OutputArray _fgmask, double learningRate=-1);
384     void create_ocl_apply_kernel();
385 };
386 
387 struct GaussBGStatModel2Params
388 {
389     //image info
390     int nWidth;
391     int nHeight;
392     int nND;//number of data dimensions (image channels)
393 
394     bool bPostFiltering;//defult 1 - do postfiltering - will make shadow detection results also give value 255
395     double  minArea; // for postfiltering
396 
397     bool bInit;//default 1, faster updates at start
398 
399     /////////////////////////
400     //very important parameters - things you will change
401     ////////////////////////
402     float fAlphaT;
403     //alpha - speed of update - if the time interval you want to average over is T
404     //set alpha=1/T. It is also usefull at start to make T slowly increase
405     //from 1 until the desired T
406     float fTb;
407     //Tb - threshold on the squared Mahalan. dist. to decide if it is well described
408     //by the background model or not. Related to Cthr from the paper.
409     //This does not influence the update of the background. A typical value could be 4 sigma
410     //and that is Tb=4*4=16;
411 
412     /////////////////////////
413     //less important parameters - things you might change but be carefull
414     ////////////////////////
415     float fTg;
416     //Tg - threshold on the squared Mahalan. dist. to decide
417     //when a sample is close to the existing components. If it is not close
418     //to any a new component will be generated. I use 3 sigma => Tg=3*3=9.
419     //Smaller Tg leads to more generated components and higher Tg might make
420     //lead to small number of components but they can grow too large
421     float fTB;//1-cf from the paper
422     //TB - threshold when the component becomes significant enough to be included into
423     //the background model. It is the TB=1-cf from the paper. So I use cf=0.1 => TB=0.
424     //For alpha=0.001 it means that the mode should exist for approximately 105 frames before
425     //it is considered foreground
426     float fVarInit;
427     float fVarMax;
428     float fVarMin;
429     //initial standard deviation  for the newly generated components.
430     //It will will influence the speed of adaptation. A good guess should be made.
431     //A simple way is to estimate the typical standard deviation from the images.
432     //I used here 10 as a reasonable value
433     float fCT;//CT - complexity reduction prior
434     //this is related to the number of samples needed to accept that a component
435     //actually exists. We use CT=0.05 of all the samples. By setting CT=0 you get
436     //the standard Stauffer&Grimson algorithm (maybe not exact but very similar)
437 
438     //even less important parameters
439     int nM;//max number of modes - const - 4 is usually enough
440 
441     //shadow detection parameters
442     bool bShadowDetection;//default 1 - do shadow detection
443     unsigned char nShadowDetection;//do shadow detection - insert this value as the detection result
444     float fTau;
445     // Tau - shadow threshold. The shadow is detected if the pixel is darker
446     //version of the background. Tau is a threshold on how much darker the shadow can be.
447     //Tau= 0.5 means that if pixel is more than 2 times darker then it is not shadow
448     //See: Prati,Mikic,Trivedi,Cucchiarra,"Detecting Moving Shadows...",IEEE PAMI,2003.
449 };
450 
451 struct GMM
452 {
453     float weight;
454     float variance;
455 };
456 
457 // shadow detection performed per pixel
458 // should work for rgb data, could be usefull for gray scale and depth data as well
459 // See: Prati,Mikic,Trivedi,Cucchiarra,"Detecting Moving Shadows...",IEEE PAMI,2003.
460 CV_INLINE bool
detectShadowGMM(const float * data,int nchannels,int nmodes,const GMM * gmm,const float * mean,float Tb,float TB,float tau)461 detectShadowGMM(const float* data, int nchannels, int nmodes,
462                 const GMM* gmm, const float* mean,
463                 float Tb, float TB, float tau)
464 {
465     float tWeight = 0;
466 
467     // check all the components  marked as background:
468     for( int mode = 0; mode < nmodes; mode++, mean += nchannels )
469     {
470         GMM g = gmm[mode];
471 
472         float numerator = 0.0f;
473         float denominator = 0.0f;
474         for( int c = 0; c < nchannels; c++ )
475         {
476             numerator   += data[c] * mean[c];
477             denominator += mean[c] * mean[c];
478         }
479 
480         // no division by zero allowed
481         if( denominator == 0 )
482             return false;
483 
484         // if tau < a < 1 then also check the color distortion
485         if( numerator <= denominator && numerator >= tau*denominator )
486         {
487             float a = numerator / denominator;
488             float dist2a = 0.0f;
489 
490             for( int c = 0; c < nchannels; c++ )
491             {
492                 float dD= a*mean[c] - data[c];
493                 dist2a += dD*dD;
494             }
495 
496             if (dist2a < Tb*g.variance*a*a)
497                 return true;
498         };
499 
500         tWeight += g.weight;
501         if( tWeight > TB )
502             return false;
503     };
504     return false;
505 }
506 
507 //update GMM - the base update function performed per pixel
508 //
509 //"Efficient Adaptive Density Estimapion per Image Pixel for the Task of Background Subtraction"
510 //Z.Zivkovic, F. van der Heijden
511 //Pattern Recognition Letters, vol. 27, no. 7, pages 773-780, 2006.
512 //
513 //The algorithm similar to the standard Stauffer&Grimson algorithm with
514 //additional selection of the number of the Gaussian components based on:
515 //
516 //"Recursive unsupervised learning of finite mixture models "
517 //Z.Zivkovic, F.van der Heijden
518 //IEEE Trans. on Pattern Analysis and Machine Intelligence, vol.26, no.5, pages 651-656, 2004
519 //http://www.zoranz.net/Publications/zivkovic2004PAMI.pdf
520 
521 class MOG2Invoker : public ParallelLoopBody
522 {
523 public:
MOG2Invoker(const Mat & _src,Mat & _dst,GMM * _gmm,float * _mean,uchar * _modesUsed,int _nmixtures,float _alphaT,float _Tb,float _TB,float _Tg,float _varInit,float _varMin,float _varMax,float _prune,float _tau,bool _detectShadows,uchar _shadowVal)524     MOG2Invoker(const Mat& _src, Mat& _dst,
525                 GMM* _gmm, float* _mean,
526                 uchar* _modesUsed,
527                 int _nmixtures, float _alphaT,
528                 float _Tb, float _TB, float _Tg,
529                 float _varInit, float _varMin, float _varMax,
530                 float _prune, float _tau, bool _detectShadows,
531                 uchar _shadowVal)
532     {
533         src = &_src;
534         dst = &_dst;
535         gmm0 = _gmm;
536         mean0 = _mean;
537         modesUsed0 = _modesUsed;
538         nmixtures = _nmixtures;
539         alphaT = _alphaT;
540         Tb = _Tb;
541         TB = _TB;
542         Tg = _Tg;
543         varInit = _varInit;
544         varMin = MIN(_varMin, _varMax);
545         varMax = MAX(_varMin, _varMax);
546         prune = _prune;
547         tau = _tau;
548         detectShadows = _detectShadows;
549         shadowVal = _shadowVal;
550     }
551 
operator ()(const Range & range) const552     void operator()(const Range& range) const
553     {
554         int y0 = range.start, y1 = range.end;
555         int ncols = src->cols, nchannels = src->channels();
556         AutoBuffer<float> buf(src->cols*nchannels);
557         float alpha1 = 1.f - alphaT;
558         float dData[CV_CN_MAX];
559 
560         for( int y = y0; y < y1; y++ )
561         {
562             const float* data = buf;
563             if( src->depth() != CV_32F )
564                 src->row(y).convertTo(Mat(1, ncols, CV_32FC(nchannels), (void*)data), CV_32F);
565             else
566                 data = src->ptr<float>(y);
567 
568             float* mean = mean0 + ncols*nmixtures*nchannels*y;
569             GMM* gmm = gmm0 + ncols*nmixtures*y;
570             uchar* modesUsed = modesUsed0 + ncols*y;
571             uchar* mask = dst->ptr(y);
572 
573             for( int x = 0; x < ncols; x++, data += nchannels, gmm += nmixtures, mean += nmixtures*nchannels )
574             {
575                 //calculate distances to the modes (+ sort)
576                 //here we need to go in descending order!!!
577                 bool background = false;//return value -> true - the pixel classified as background
578 
579                 //internal:
580                 bool fitsPDF = false;//if it remains zero a new GMM mode will be added
581                 int nmodes = modesUsed[x], nNewModes = nmodes;//current number of modes in GMM
582                 float totalWeight = 0.f;
583 
584                 float* mean_m = mean;
585 
586                 //////
587                 //go through all modes
588                 for( int mode = 0; mode < nmodes; mode++, mean_m += nchannels )
589                 {
590                     float weight = alpha1*gmm[mode].weight + prune;//need only weight if fit is found
591                     int swap_count = 0;
592                     ////
593                     //fit not found yet
594                     if( !fitsPDF )
595                     {
596                         //check if it belongs to some of the remaining modes
597                         float var = gmm[mode].variance;
598 
599                         //calculate difference and distance
600                         float dist2;
601 
602                         if( nchannels == 3 )
603                         {
604                             dData[0] = mean_m[0] - data[0];
605                             dData[1] = mean_m[1] - data[1];
606                             dData[2] = mean_m[2] - data[2];
607                             dist2 = dData[0]*dData[0] + dData[1]*dData[1] + dData[2]*dData[2];
608                         }
609                         else
610                         {
611                             dist2 = 0.f;
612                             for( int c = 0; c < nchannels; c++ )
613                             {
614                                 dData[c] = mean_m[c] - data[c];
615                                 dist2 += dData[c]*dData[c];
616                             }
617                         }
618 
619                         //background? - Tb - usually larger than Tg
620                         if( totalWeight < TB && dist2 < Tb*var )
621                             background = true;
622 
623                         //check fit
624                         if( dist2 < Tg*var )
625                         {
626                             /////
627                             //belongs to the mode
628                             fitsPDF = true;
629 
630                             //update distribution
631 
632                             //update weight
633                             weight += alphaT;
634                             float k = alphaT/weight;
635 
636                             //update mean
637                             for( int c = 0; c < nchannels; c++ )
638                                 mean_m[c] -= k*dData[c];
639 
640                             //update variance
641                             float varnew = var + k*(dist2-var);
642                             //limit the variance
643                             varnew = MAX(varnew, varMin);
644                             varnew = MIN(varnew, varMax);
645                             gmm[mode].variance = varnew;
646 
647                             //sort
648                             //all other weights are at the same place and
649                             //only the matched (iModes) is higher -> just find the new place for it
650                             for( int i = mode; i > 0; i-- )
651                             {
652                                 //check one up
653                                 if( weight < gmm[i-1].weight )
654                                     break;
655 
656                                 swap_count++;
657                                 //swap one up
658                                 std::swap(gmm[i], gmm[i-1]);
659                                 for( int c = 0; c < nchannels; c++ )
660                                     std::swap(mean[i*nchannels + c], mean[(i-1)*nchannels + c]);
661                             }
662                             //belongs to the mode - bFitsPDF becomes 1
663                             /////
664                         }
665                     }//!bFitsPDF)
666 
667                     //check prune
668                     if( weight < -prune )
669                     {
670                         weight = 0.0;
671                         nmodes--;
672                     }
673 
674                     gmm[mode-swap_count].weight = weight;//update weight by the calculated value
675                     totalWeight += weight;
676                 }
677                 //go through all modes
678                 //////
679 
680                 //renormalize weights
681                 totalWeight = 1.f/totalWeight;
682                 for( int mode = 0; mode < nmodes; mode++ )
683                 {
684                     gmm[mode].weight *= totalWeight;
685                 }
686 
687                 nmodes = nNewModes;
688 
689                 //make new mode if needed and exit
690                 if( !fitsPDF && alphaT > 0.f )
691                 {
692                     // replace the weakest or add a new one
693                     int mode = nmodes == nmixtures ? nmixtures-1 : nmodes++;
694 
695                     if (nmodes==1)
696                         gmm[mode].weight = 1.f;
697                     else
698                     {
699                         gmm[mode].weight = alphaT;
700 
701                         // renormalize all other weights
702                         for( int i = 0; i < nmodes-1; i++ )
703                             gmm[i].weight *= alpha1;
704                     }
705 
706                     // init
707                     for( int c = 0; c < nchannels; c++ )
708                         mean[mode*nchannels + c] = data[c];
709 
710                     gmm[mode].variance = varInit;
711 
712                     //sort
713                     //find the new place for it
714                     for( int i = nmodes - 1; i > 0; i-- )
715                     {
716                         // check one up
717                         if( alphaT < gmm[i-1].weight )
718                             break;
719 
720                         // swap one up
721                         std::swap(gmm[i], gmm[i-1]);
722                         for( int c = 0; c < nchannels; c++ )
723                             std::swap(mean[i*nchannels + c], mean[(i-1)*nchannels + c]);
724                     }
725                 }
726 
727                 //set the number of modes
728                 modesUsed[x] = uchar(nmodes);
729                 mask[x] = background ? 0 :
730                     detectShadows && detectShadowGMM(data, nchannels, nmodes, gmm, mean, Tb, TB, tau) ?
731                     shadowVal : 255;
732             }
733         }
734     }
735 
736     const Mat* src;
737     Mat* dst;
738     GMM* gmm0;
739     float* mean0;
740     uchar* modesUsed0;
741 
742     int nmixtures;
743     float alphaT, Tb, TB, Tg;
744     float varInit, varMin, varMax, prune, tau;
745 
746     bool detectShadows;
747     uchar shadowVal;
748 };
749 
750 #ifdef HAVE_OPENCL
751 
ocl_apply(InputArray _image,OutputArray _fgmask,double learningRate)752 bool BackgroundSubtractorMOG2Impl::ocl_apply(InputArray _image, OutputArray _fgmask, double learningRate)
753 {
754     ++nframes;
755     learningRate = learningRate >= 0 && nframes > 1 ? learningRate : 1./std::min( 2*nframes, history );
756     CV_Assert(learningRate >= 0);
757 
758     _fgmask.create(_image.size(), CV_8U);
759     UMat fgmask = _fgmask.getUMat();
760 
761     const double alpha1 = 1.0f - learningRate;
762 
763     UMat frame = _image.getUMat();
764 
765     float varMax = MAX(fVarMin, fVarMax);
766     float varMin = MIN(fVarMin, fVarMax);
767 
768     int idxArg = 0;
769     idxArg = kernel_apply.set(idxArg, ocl::KernelArg::ReadOnly(frame));
770     idxArg = kernel_apply.set(idxArg, ocl::KernelArg::PtrReadWrite(u_bgmodelUsedModes));
771     idxArg = kernel_apply.set(idxArg, ocl::KernelArg::PtrReadWrite(u_weight));
772     idxArg = kernel_apply.set(idxArg, ocl::KernelArg::PtrReadWrite(u_mean));
773     idxArg = kernel_apply.set(idxArg, ocl::KernelArg::PtrReadWrite(u_variance));
774     idxArg = kernel_apply.set(idxArg, ocl::KernelArg::WriteOnlyNoSize(fgmask));
775 
776     idxArg = kernel_apply.set(idxArg, (float)learningRate);        //alphaT
777     idxArg = kernel_apply.set(idxArg, (float)alpha1);
778     idxArg = kernel_apply.set(idxArg, (float)(-learningRate*fCT));   //prune
779 
780     idxArg = kernel_apply.set(idxArg, (float)varThreshold); //c_Tb
781     idxArg = kernel_apply.set(idxArg, backgroundRatio);     //c_TB
782     idxArg = kernel_apply.set(idxArg, varThresholdGen);     //c_Tg
783     idxArg = kernel_apply.set(idxArg, varMin);
784     idxArg = kernel_apply.set(idxArg, varMax);
785     idxArg = kernel_apply.set(idxArg, fVarInit);
786     idxArg = kernel_apply.set(idxArg, fTau);
787     if (bShadowDetection)
788         kernel_apply.set(idxArg, nShadowDetection);
789 
790     size_t globalsize[] = {frame.cols, frame.rows, 1};
791     return kernel_apply.run(2, globalsize, NULL, true);
792 }
793 
ocl_getBackgroundImage(OutputArray _backgroundImage) const794 bool BackgroundSubtractorMOG2Impl::ocl_getBackgroundImage(OutputArray _backgroundImage) const
795 {
796     CV_Assert(frameType == CV_8UC1 || frameType == CV_8UC3);
797 
798     _backgroundImage.create(frameSize, frameType);
799     UMat dst = _backgroundImage.getUMat();
800 
801     int idxArg = 0;
802     idxArg = kernel_getBg.set(idxArg, ocl::KernelArg::PtrReadOnly(u_bgmodelUsedModes));
803     idxArg = kernel_getBg.set(idxArg, ocl::KernelArg::PtrReadOnly(u_weight));
804     idxArg = kernel_getBg.set(idxArg, ocl::KernelArg::PtrReadOnly(u_mean));
805     idxArg = kernel_getBg.set(idxArg, ocl::KernelArg::WriteOnly(dst));
806     kernel_getBg.set(idxArg, backgroundRatio);
807 
808     size_t globalsize[2] = {u_bgmodelUsedModes.cols, u_bgmodelUsedModes.rows};
809 
810     return kernel_getBg.run(2, globalsize, NULL, false);
811 }
812 
813 #endif
814 
create_ocl_apply_kernel()815 void BackgroundSubtractorMOG2Impl::create_ocl_apply_kernel()
816 {
817     int nchannels = CV_MAT_CN(frameType);
818     String opts = format("-D CN=%d -D NMIXTURES=%d%s", nchannels, nmixtures, bShadowDetection ? " -D SHADOW_DETECT" : "");
819     kernel_apply.create("mog2_kernel", ocl::video::bgfg_mog2_oclsrc, opts);
820 }
821 
apply(InputArray _image,OutputArray _fgmask,double learningRate)822 void BackgroundSubtractorMOG2Impl::apply(InputArray _image, OutputArray _fgmask, double learningRate)
823 {
824     bool needToInitialize = nframes == 0 || learningRate >= 1 || _image.size() != frameSize || _image.type() != frameType;
825 
826     if( needToInitialize )
827         initialize(_image.size(), _image.type());
828 
829     if (opencl_ON)
830     {
831         CV_OCL_RUN(opencl_ON, ocl_apply(_image, _fgmask, learningRate))
832 
833         opencl_ON = false;
834         initialize(_image.size(), _image.type());
835     }
836 
837     Mat image = _image.getMat();
838     _fgmask.create( image.size(), CV_8U );
839     Mat fgmask = _fgmask.getMat();
840 
841     ++nframes;
842     learningRate = learningRate >= 0 && nframes > 1 ? learningRate : 1./std::min( 2*nframes, history );
843     CV_Assert(learningRate >= 0);
844 
845     parallel_for_(Range(0, image.rows),
846                   MOG2Invoker(image, fgmask,
847                               bgmodel.ptr<GMM>(),
848                               (float*)(bgmodel.ptr() + sizeof(GMM)*nmixtures*image.rows*image.cols),
849                               bgmodelUsedModes.ptr(), nmixtures, (float)learningRate,
850                               (float)varThreshold,
851                               backgroundRatio, varThresholdGen,
852                               fVarInit, fVarMin, fVarMax, float(-learningRate*fCT), fTau,
853                               bShadowDetection, nShadowDetection),
854                               image.total()/(double)(1 << 16));
855 }
856 
getBackgroundImage(OutputArray backgroundImage) const857 void BackgroundSubtractorMOG2Impl::getBackgroundImage(OutputArray backgroundImage) const
858 {
859     if (opencl_ON)
860     {
861         CV_OCL_RUN(opencl_ON, ocl_getBackgroundImage(backgroundImage))
862 
863         opencl_ON = false;
864         return;
865     }
866 
867     int nchannels = CV_MAT_CN(frameType);
868     CV_Assert(nchannels == 1 || nchannels == 3);
869     Mat meanBackground(frameSize, CV_MAKETYPE(CV_8U, nchannels), Scalar::all(0));
870     int firstGaussianIdx = 0;
871     const GMM* gmm = bgmodel.ptr<GMM>();
872     const float* mean = reinterpret_cast<const float*>(gmm + frameSize.width*frameSize.height*nmixtures);
873     std::vector<float> meanVal(nchannels, 0.f);
874     for(int row=0; row<meanBackground.rows; row++)
875     {
876         for(int col=0; col<meanBackground.cols; col++)
877         {
878             int nmodes = bgmodelUsedModes.at<uchar>(row, col);
879             float totalWeight = 0.f;
880             for(int gaussianIdx = firstGaussianIdx; gaussianIdx < firstGaussianIdx + nmodes; gaussianIdx++)
881             {
882                 GMM gaussian = gmm[gaussianIdx];
883                 size_t meanPosition = gaussianIdx*nchannels;
884                 for(int chn = 0; chn < nchannels; chn++)
885                 {
886                     meanVal[chn] += gaussian.weight * mean[meanPosition + chn];
887                 }
888                 totalWeight += gaussian.weight;
889 
890                 if(totalWeight > backgroundRatio)
891                     break;
892             }
893             float invWeight = 1.f/totalWeight;
894             switch(nchannels)
895             {
896             case 1:
897                 meanBackground.at<uchar>(row, col) = (uchar)(meanVal[0] * invWeight);
898                 meanVal[0] = 0.f;
899                 break;
900             case 3:
901                 Vec3f& meanVec = *reinterpret_cast<Vec3f*>(&meanVal[0]);
902                 meanBackground.at<Vec3b>(row, col) = Vec3b(meanVec * invWeight);
903                 meanVec = 0.f;
904                 break;
905             }
906             firstGaussianIdx += nmixtures;
907         }
908     }
909     meanBackground.copyTo(backgroundImage);
910 }
911 
createBackgroundSubtractorMOG2(int _history,double _varThreshold,bool _bShadowDetection)912 Ptr<BackgroundSubtractorMOG2> createBackgroundSubtractorMOG2(int _history, double _varThreshold,
913                                                              bool _bShadowDetection)
914 {
915     return makePtr<BackgroundSubtractorMOG2Impl>(_history, (float)_varThreshold, _bShadowDetection);
916 }
917 
918 }
919 
920 /* End of file. */
921