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 Gaussian 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 Gaussian 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 CV_FINAL : 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 #ifdef HAVE_OPENCL
146 opencl_ON = true;
147 #endif
148 }
149 //! the full constructor that takes the length of the history,
150 // the number of gaussian mixtures, the background ratio parameter and the noise strength
BackgroundSubtractorMOG2Impl(int _history,float _varThreshold,bool _bShadowDetection=true)151 BackgroundSubtractorMOG2Impl(int _history, float _varThreshold, bool _bShadowDetection=true)
152 {
153 frameSize = Size(0,0);
154 frameType = 0;
155
156 nframes = 0;
157 history = _history > 0 ? _history : defaultHistory2;
158 varThreshold = (_varThreshold>0)? _varThreshold : defaultVarThreshold2;
159 bShadowDetection = _bShadowDetection;
160
161 nmixtures = defaultNMixtures2;
162 backgroundRatio = defaultBackgroundRatio2;
163 fVarInit = defaultVarInit2;
164 fVarMax = defaultVarMax2;
165 fVarMin = defaultVarMin2;
166
167 varThresholdGen = defaultVarThresholdGen2;
168 fCT = defaultfCT2;
169 nShadowDetection = defaultnShadowDetection2;
170 fTau = defaultfTau;
171 name_ = "BackgroundSubtractor.MOG2";
172 #ifdef HAVE_OPENCL
173 opencl_ON = true;
174 #endif
175 }
176 //! the destructor
~BackgroundSubtractorMOG2Impl()177 ~BackgroundSubtractorMOG2Impl() CV_OVERRIDE {}
178 //! the update operator
179 void apply(InputArray image, OutputArray fgmask, double learningRate) CV_OVERRIDE;
180
181 //! computes a background image which are the mean of all background gaussians
182 virtual void getBackgroundImage(OutputArray backgroundImage) const CV_OVERRIDE;
183
184 //! re-initialization method
initialize(Size _frameSize,int _frameType)185 void initialize(Size _frameSize, int _frameType)
186 {
187 frameSize = _frameSize;
188 frameType = _frameType;
189 nframes = 0;
190
191 int nchannels = CV_MAT_CN(frameType);
192 CV_Assert( nchannels <= CV_CN_MAX );
193 CV_Assert( nmixtures <= 255);
194
195 #ifdef HAVE_OPENCL
196 if (ocl::isOpenCLActivated() && opencl_ON)
197 {
198 create_ocl_apply_kernel();
199
200 bool isFloat = CV_MAKETYPE(CV_32F,nchannels) == frameType;
201 kernel_getBg.create("getBackgroundImage2_kernel", ocl::video::bgfg_mog2_oclsrc, format( "-D CN=%d -D FL=%d -D NMIXTURES=%d", nchannels, isFloat, nmixtures));
202
203 if (kernel_apply.empty() || kernel_getBg.empty())
204 opencl_ON = false;
205 }
206 else opencl_ON = false;
207
208 if (opencl_ON)
209 {
210 u_weight.create(frameSize.height * nmixtures, frameSize.width, CV_32FC1);
211 u_weight.setTo(Scalar::all(0));
212
213 u_variance.create(frameSize.height * nmixtures, frameSize.width, CV_32FC1);
214 u_variance.setTo(Scalar::all(0));
215
216 if (nchannels==3)
217 nchannels=4;
218 u_mean.create(frameSize.height * nmixtures, frameSize.width, CV_32FC(nchannels)); //4 channels
219 u_mean.setTo(Scalar::all(0));
220
221 //make the array for keeping track of the used modes per pixel - all zeros at start
222 u_bgmodelUsedModes.create(frameSize, CV_8UC1);
223 u_bgmodelUsedModes.setTo(cv::Scalar::all(0));
224 }
225 else
226 #endif
227 {
228 // for each gaussian mixture of each pixel bg model we store ...
229 // the mixture weight (w),
230 // the mean (nchannels values) and
231 // the covariance
232 bgmodel.create( 1, frameSize.height*frameSize.width*nmixtures*(2 + nchannels), CV_32F );
233 //make the array for keeping track of the used modes per pixel - all zeros at start
234 bgmodelUsedModes.create(frameSize,CV_8U);
235 bgmodelUsedModes = Scalar::all(0);
236 }
237 }
238
getHistory() const239 virtual int getHistory() const CV_OVERRIDE { return history; }
setHistory(int _nframes)240 virtual void setHistory(int _nframes) CV_OVERRIDE { history = _nframes; }
241
getNMixtures() const242 virtual int getNMixtures() const CV_OVERRIDE { return nmixtures; }
setNMixtures(int nmix)243 virtual void setNMixtures(int nmix) CV_OVERRIDE { nmixtures = nmix; }
244
getBackgroundRatio() const245 virtual double getBackgroundRatio() const CV_OVERRIDE { return backgroundRatio; }
setBackgroundRatio(double _backgroundRatio)246 virtual void setBackgroundRatio(double _backgroundRatio) CV_OVERRIDE { backgroundRatio = (float)_backgroundRatio; }
247
getVarThreshold() const248 virtual double getVarThreshold() const CV_OVERRIDE { return varThreshold; }
setVarThreshold(double _varThreshold)249 virtual void setVarThreshold(double _varThreshold) CV_OVERRIDE { varThreshold = _varThreshold; }
250
getVarThresholdGen() const251 virtual double getVarThresholdGen() const CV_OVERRIDE { return varThresholdGen; }
setVarThresholdGen(double _varThresholdGen)252 virtual void setVarThresholdGen(double _varThresholdGen) CV_OVERRIDE { varThresholdGen = (float)_varThresholdGen; }
253
getVarInit() const254 virtual double getVarInit() const CV_OVERRIDE { return fVarInit; }
setVarInit(double varInit)255 virtual void setVarInit(double varInit) CV_OVERRIDE { fVarInit = (float)varInit; }
256
getVarMin() const257 virtual double getVarMin() const CV_OVERRIDE { return fVarMin; }
setVarMin(double varMin)258 virtual void setVarMin(double varMin) CV_OVERRIDE { fVarMin = (float)varMin; }
259
getVarMax() const260 virtual double getVarMax() const CV_OVERRIDE { return fVarMax; }
setVarMax(double varMax)261 virtual void setVarMax(double varMax) CV_OVERRIDE { fVarMax = (float)varMax; }
262
getComplexityReductionThreshold() const263 virtual double getComplexityReductionThreshold() const CV_OVERRIDE { return fCT; }
setComplexityReductionThreshold(double ct)264 virtual void setComplexityReductionThreshold(double ct) CV_OVERRIDE { fCT = (float)ct; }
265
getDetectShadows() const266 virtual bool getDetectShadows() const CV_OVERRIDE { return bShadowDetection; }
setDetectShadows(bool detectshadows)267 virtual void setDetectShadows(bool detectshadows) CV_OVERRIDE
268 {
269 if (bShadowDetection == detectshadows)
270 return;
271 bShadowDetection = detectshadows;
272 #ifdef HAVE_OPENCL
273 if (!kernel_apply.empty())
274 {
275 create_ocl_apply_kernel();
276 CV_Assert( !kernel_apply.empty() );
277 }
278 #endif
279 }
280
getShadowValue() const281 virtual int getShadowValue() const CV_OVERRIDE { return nShadowDetection; }
setShadowValue(int value)282 virtual void setShadowValue(int value) CV_OVERRIDE { nShadowDetection = (uchar)value; }
283
getShadowThreshold() const284 virtual double getShadowThreshold() const CV_OVERRIDE { return fTau; }
setShadowThreshold(double value)285 virtual void setShadowThreshold(double value) CV_OVERRIDE { fTau = (float)value; }
286
write(FileStorage & fs) const287 virtual void write(FileStorage& fs) const CV_OVERRIDE
288 {
289 writeFormat(fs);
290 fs << "name" << name_
291 << "history" << history
292 << "nmixtures" << nmixtures
293 << "backgroundRatio" << backgroundRatio
294 << "varThreshold" << varThreshold
295 << "varThresholdGen" << varThresholdGen
296 << "varInit" << fVarInit
297 << "varMin" << fVarMin
298 << "varMax" << fVarMax
299 << "complexityReductionThreshold" << fCT
300 << "detectShadows" << (int)bShadowDetection
301 << "shadowValue" << (int)nShadowDetection
302 << "shadowThreshold" << fTau;
303 }
304
read(const FileNode & fn)305 virtual void read(const FileNode& fn) CV_OVERRIDE
306 {
307 CV_Assert( (String)fn["name"] == name_ );
308 history = (int)fn["history"];
309 nmixtures = (int)fn["nmixtures"];
310 backgroundRatio = (float)fn["backgroundRatio"];
311 varThreshold = (double)fn["varThreshold"];
312 varThresholdGen = (float)fn["varThresholdGen"];
313 fVarInit = (float)fn["varInit"];
314 fVarMin = (float)fn["varMin"];
315 fVarMax = (float)fn["varMax"];
316 fCT = (float)fn["complexityReductionThreshold"];
317 bShadowDetection = (int)fn["detectShadows"] != 0;
318 nShadowDetection = saturate_cast<uchar>((int)fn["shadowValue"]);
319 fTau = (float)fn["shadowThreshold"];
320 }
321
322 protected:
323 Size frameSize;
324 int frameType;
325 Mat bgmodel;
326 Mat bgmodelUsedModes;//keep track of number of modes per pixel
327
328 #ifdef HAVE_OPENCL
329 //for OCL
330
331 mutable bool opencl_ON;
332
333 UMat u_weight;
334 UMat u_variance;
335 UMat u_mean;
336 UMat u_bgmodelUsedModes;
337
338 mutable ocl::Kernel kernel_apply;
339 mutable ocl::Kernel kernel_getBg;
340 #endif
341
342 int nframes;
343 int history;
344 int nmixtures;
345 //! here it is the maximum allowed number of mixture components.
346 //! Actual number is determined dynamically per pixel
347 double varThreshold;
348 // threshold on the squared Mahalanobis distance to decide if it is well described
349 // by the background model or not. Related to Cthr from the paper.
350 // This does not influence the update of the background. A typical value could be 4 sigma
351 // and that is varThreshold=4*4=16; Corresponds to Tb in the paper.
352
353 /////////////////////////
354 // less important parameters - things you might change but be careful
355 ////////////////////////
356 float backgroundRatio;
357 // corresponds to fTB=1-cf from the paper
358 // TB - threshold when the component becomes significant enough to be included into
359 // the background model. It is the TB=1-cf from the paper. So I use cf=0.1 => TB=0.
360 // For alpha=0.001 it means that the mode should exist for approximately 105 frames before
361 // it is considered foreground
362 // float noiseSigma;
363 float varThresholdGen;
364 //correspondts to Tg - threshold on the squared Mahalan. dist. to decide
365 //when a sample is close to the existing components. If it is not close
366 //to any a new component will be generated. I use 3 sigma => Tg=3*3=9.
367 //Smaller Tg leads to more generated components and higher Tg might make
368 //lead to small number of components but they can grow too large
369 float fVarInit;
370 float fVarMin;
371 float fVarMax;
372 //initial variance for the newly generated components.
373 //It will will influence the speed of adaptation. A good guess should be made.
374 //A simple way is to estimate the typical standard deviation from the images.
375 //I used here 10 as a reasonable value
376 // min and max can be used to further control the variance
377 float fCT;//CT - complexity reduction prior
378 //this is related to the number of samples needed to accept that a component
379 //actually exists. We use CT=0.05 of all the samples. By setting CT=0 you get
380 //the standard Stauffer&Grimson algorithm (maybe not exact but very similar)
381
382 //shadow detection parameters
383 bool bShadowDetection;//default 1 - do shadow detection
384 unsigned char nShadowDetection;//do shadow detection - insert this value as the detection result - 127 default value
385 float fTau;
386 // Tau - shadow threshold. The shadow is detected if the pixel is darker
387 //version of the background. Tau is a threshold on how much darker the shadow can be.
388 //Tau= 0.5 means that if pixel is more than 2 times darker then it is not shadow
389 //See: Prati,Mikic,Trivedi,Cucchiara,"Detecting Moving Shadows...",IEEE PAMI,2003.
390
391 String name_;
392
393 template <typename T, int CN>
394 void getBackgroundImage_intern(OutputArray backgroundImage) const;
395
396 #ifdef HAVE_OPENCL
397 bool ocl_getBackgroundImage(OutputArray backgroundImage) const;
398 bool ocl_apply(InputArray _image, OutputArray _fgmask, double learningRate=-1);
399 void create_ocl_apply_kernel();
400 #endif
401 };
402
403 struct GaussBGStatModel2Params
404 {
405 //image info
406 int nWidth;
407 int nHeight;
408 int nND;//number of data dimensions (image channels)
409
410 bool bPostFiltering;//default 1 - do postfiltering - will make shadow detection results also give value 255
411 double minArea; // for postfiltering
412
413 bool bInit;//default 1, faster updates at start
414
415 /////////////////////////
416 //very important parameters - things you will change
417 ////////////////////////
418 float fAlphaT;
419 //alpha - speed of update - if the time interval you want to average over is T
420 //set alpha=1/T. It is also useful at start to make T slowly increase
421 //from 1 until the desired T
422 float fTb;
423 //Tb - threshold on the squared Mahalan. dist. to decide if it is well described
424 //by the background model or not. Related to Cthr from the paper.
425 //This does not influence the update of the background. A typical value could be 4 sigma
426 //and that is Tb=4*4=16;
427
428 /////////////////////////
429 //less important parameters - things you might change but be careful
430 ////////////////////////
431 float fTg;
432 //Tg - threshold on the squared Mahalan. dist. to decide
433 //when a sample is close to the existing components. If it is not close
434 //to any a new component will be generated. I use 3 sigma => Tg=3*3=9.
435 //Smaller Tg leads to more generated components and higher Tg might make
436 //lead to small number of components but they can grow too large
437 float fTB;//1-cf from the paper
438 //TB - threshold when the component becomes significant enough to be included into
439 //the background model. It is the TB=1-cf from the paper. So I use cf=0.1 => TB=0.
440 //For alpha=0.001 it means that the mode should exist for approximately 105 frames before
441 //it is considered foreground
442 float fVarInit;
443 float fVarMax;
444 float fVarMin;
445 //initial standard deviation for the newly generated components.
446 //It will will influence the speed of adaptation. A good guess should be made.
447 //A simple way is to estimate the typical standard deviation from the images.
448 //I used here 10 as a reasonable value
449 float fCT;//CT - complexity reduction prior
450 //this is related to the number of samples needed to accept that a component
451 //actually exists. We use CT=0.05 of all the samples. By setting CT=0 you get
452 //the standard Stauffer&Grimson algorithm (maybe not exact but very similar)
453
454 //even less important parameters
455 int nM;//max number of modes - const - 4 is usually enough
456
457 //shadow detection parameters
458 bool bShadowDetection;//default 1 - do shadow detection
459 unsigned char nShadowDetection;//do shadow detection - insert this value as the detection result
460 float fTau;
461 // Tau - shadow threshold. The shadow is detected if the pixel is darker
462 //version of the background. Tau is a threshold on how much darker the shadow can be.
463 //Tau= 0.5 means that if pixel is more than 2 times darker then it is not shadow
464 //See: Prati,Mikic,Trivedi,Cucchiara,"Detecting Moving Shadows...",IEEE PAMI,2003.
465 };
466
467 struct GMM
468 {
469 float weight;
470 float variance;
471 };
472
473 // shadow detection performed per pixel
474 // should work for rgb data, could be useful for gray scale and depth data as well
475 // See: Prati,Mikic,Trivedi,Cucchiara,"Detecting Moving Shadows...",IEEE PAMI,2003.
476 CV_INLINE bool
detectShadowGMM(const float * data,int nchannels,int nmodes,const GMM * gmm,const float * mean,float Tb,float TB,float tau)477 detectShadowGMM(const float* data, int nchannels, int nmodes,
478 const GMM* gmm, const float* mean,
479 float Tb, float TB, float tau)
480 {
481 float tWeight = 0;
482
483 // check all the components marked as background:
484 for( int mode = 0; mode < nmodes; mode++, mean += nchannels )
485 {
486 GMM g = gmm[mode];
487
488 float numerator = 0.0f;
489 float denominator = 0.0f;
490 for( int c = 0; c < nchannels; c++ )
491 {
492 numerator += data[c] * mean[c];
493 denominator += mean[c] * mean[c];
494 }
495
496 // no division by zero allowed
497 if( denominator == 0 )
498 return false;
499
500 // if tau < a < 1 then also check the color distortion
501 if( numerator <= denominator && numerator >= tau*denominator )
502 {
503 float a = numerator / denominator;
504 float dist2a = 0.0f;
505
506 for( int c = 0; c < nchannels; c++ )
507 {
508 float dD= a*mean[c] - data[c];
509 dist2a += dD*dD;
510 }
511
512 if (dist2a < Tb*g.variance*a*a)
513 return true;
514 };
515
516 tWeight += g.weight;
517 if( tWeight > TB )
518 return false;
519 };
520 return false;
521 }
522
523 //update GMM - the base update function performed per pixel
524 //
525 //"Efficient Adaptive Density Estimapion per Image Pixel for the Task of Background Subtraction"
526 //Z.Zivkovic, F. van der Heijden
527 //Pattern Recognition Letters, vol. 27, no. 7, pages 773-780, 2006.
528 //
529 //The algorithm similar to the standard Stauffer&Grimson algorithm with
530 //additional selection of the number of the Gaussian components based on:
531 //
532 //"Recursive unsupervised learning of finite mixture models "
533 //Z.Zivkovic, F.van der Heijden
534 //IEEE Trans. on Pattern Analysis and Machine Intelligence, vol.26, no.5, pages 651-656, 2004
535 //http://www.zoranz.net/Publications/zivkovic2004PAMI.pdf
536
537 class MOG2Invoker : public ParallelLoopBody
538 {
539 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)540 MOG2Invoker(const Mat& _src, Mat& _dst,
541 GMM* _gmm, float* _mean,
542 uchar* _modesUsed,
543 int _nmixtures, float _alphaT,
544 float _Tb, float _TB, float _Tg,
545 float _varInit, float _varMin, float _varMax,
546 float _prune, float _tau, bool _detectShadows,
547 uchar _shadowVal)
548 {
549 src = &_src;
550 dst = &_dst;
551 gmm0 = _gmm;
552 mean0 = _mean;
553 modesUsed0 = _modesUsed;
554 nmixtures = _nmixtures;
555 alphaT = _alphaT;
556 Tb = _Tb;
557 TB = _TB;
558 Tg = _Tg;
559 varInit = _varInit;
560 varMin = MIN(_varMin, _varMax);
561 varMax = MAX(_varMin, _varMax);
562 prune = _prune;
563 tau = _tau;
564 detectShadows = _detectShadows;
565 shadowVal = _shadowVal;
566 }
567
operator ()(const Range & range) const568 void operator()(const Range& range) const CV_OVERRIDE
569 {
570 int y0 = range.start, y1 = range.end;
571 int ncols = src->cols, nchannels = src->channels();
572 AutoBuffer<float> buf(src->cols*nchannels);
573 float alpha1 = 1.f - alphaT;
574 float dData[CV_CN_MAX];
575
576 for( int y = y0; y < y1; y++ )
577 {
578 const float* data = buf.data();
579 if( src->depth() != CV_32F )
580 src->row(y).convertTo(Mat(1, ncols, CV_32FC(nchannels), (void*)data), CV_32F);
581 else
582 data = src->ptr<float>(y);
583
584 float* mean = mean0 + ncols*nmixtures*nchannels*y;
585 GMM* gmm = gmm0 + ncols*nmixtures*y;
586 uchar* modesUsed = modesUsed0 + ncols*y;
587 uchar* mask = dst->ptr(y);
588
589 for( int x = 0; x < ncols; x++, data += nchannels, gmm += nmixtures, mean += nmixtures*nchannels )
590 {
591 //calculate distances to the modes (+ sort)
592 //here we need to go in descending order!!!
593 bool background = false;//return value -> true - the pixel classified as background
594
595 //internal:
596 bool fitsPDF = false;//if it remains zero a new GMM mode will be added
597 int nmodes = modesUsed[x];//current number of modes in GMM
598 float totalWeight = 0.f;
599
600 float* mean_m = mean;
601
602 //////
603 //go through all modes
604 for( int mode = 0; mode < nmodes; mode++, mean_m += nchannels )
605 {
606 float weight = alpha1*gmm[mode].weight + prune;//need only weight if fit is found
607 int swap_count = 0;
608 ////
609 //fit not found yet
610 if( !fitsPDF )
611 {
612 //check if it belongs to some of the remaining modes
613 float var = gmm[mode].variance;
614
615 //calculate difference and distance
616 float dist2;
617
618 if( nchannels == 3 )
619 {
620 dData[0] = mean_m[0] - data[0];
621 dData[1] = mean_m[1] - data[1];
622 dData[2] = mean_m[2] - data[2];
623 dist2 = dData[0]*dData[0] + dData[1]*dData[1] + dData[2]*dData[2];
624 }
625 else
626 {
627 dist2 = 0.f;
628 for( int c = 0; c < nchannels; c++ )
629 {
630 dData[c] = mean_m[c] - data[c];
631 dist2 += dData[c]*dData[c];
632 }
633 }
634
635 //background? - Tb - usually larger than Tg
636 if( totalWeight < TB && dist2 < Tb*var )
637 background = true;
638
639 //check fit
640 if( dist2 < Tg*var )
641 {
642 /////
643 //belongs to the mode
644 fitsPDF = true;
645
646 //update distribution
647
648 //update weight
649 weight += alphaT;
650 float k = alphaT/weight;
651
652 //update mean
653 for( int c = 0; c < nchannels; c++ )
654 mean_m[c] -= k*dData[c];
655
656 //update variance
657 float varnew = var + k*(dist2-var);
658 //limit the variance
659 varnew = MAX(varnew, varMin);
660 varnew = MIN(varnew, varMax);
661 gmm[mode].variance = varnew;
662
663 //sort
664 //all other weights are at the same place and
665 //only the matched (iModes) is higher -> just find the new place for it
666 for( int i = mode; i > 0; i-- )
667 {
668 //check one up
669 if( weight < gmm[i-1].weight )
670 break;
671
672 swap_count++;
673 //swap one up
674 std::swap(gmm[i], gmm[i-1]);
675 for( int c = 0; c < nchannels; c++ )
676 std::swap(mean[i*nchannels + c], mean[(i-1)*nchannels + c]);
677 }
678 //belongs to the mode - bFitsPDF becomes 1
679 /////
680 }
681 }//!bFitsPDF)
682
683 //check prune
684 if( weight < -prune )
685 {
686 weight = 0.0;
687 nmodes--;
688 }
689
690 gmm[mode-swap_count].weight = weight;//update weight by the calculated value
691 totalWeight += weight;
692 }
693 //go through all modes
694 //////
695
696 // Renormalize weights. In the special case that the pixel does
697 // not agree with any modes, set weights to zero (a new mode will be added below).
698 float invWeight = 0.f;
699 if (std::abs(totalWeight) > FLT_EPSILON) {
700 invWeight = 1.f/totalWeight;
701 }
702
703 for( int mode = 0; mode < nmodes; mode++ )
704 {
705 gmm[mode].weight *= invWeight;
706 }
707
708 //make new mode if needed and exit
709 if( !fitsPDF && alphaT > 0.f )
710 {
711 // replace the weakest or add a new one
712 int mode = nmodes == nmixtures ? nmixtures-1 : nmodes++;
713
714 if (nmodes==1)
715 gmm[mode].weight = 1.f;
716 else
717 {
718 gmm[mode].weight = alphaT;
719
720 // renormalize all other weights
721 for( int i = 0; i < nmodes-1; i++ )
722 gmm[i].weight *= alpha1;
723 }
724
725 // init
726 for( int c = 0; c < nchannels; c++ )
727 mean[mode*nchannels + c] = data[c];
728
729 gmm[mode].variance = varInit;
730
731 //sort
732 //find the new place for it
733 for( int i = nmodes - 1; i > 0; i-- )
734 {
735 // check one up
736 if( alphaT < gmm[i-1].weight )
737 break;
738
739 // swap one up
740 std::swap(gmm[i], gmm[i-1]);
741 for( int c = 0; c < nchannels; c++ )
742 std::swap(mean[i*nchannels + c], mean[(i-1)*nchannels + c]);
743 }
744 }
745
746 //set the number of modes
747 modesUsed[x] = uchar(nmodes);
748 mask[x] = background ? 0 :
749 detectShadows && detectShadowGMM(data, nchannels, nmodes, gmm, mean, Tb, TB, tau) ?
750 shadowVal : 255;
751 }
752 }
753 }
754
755 const Mat* src;
756 Mat* dst;
757 GMM* gmm0;
758 float* mean0;
759 uchar* modesUsed0;
760
761 int nmixtures;
762 float alphaT, Tb, TB, Tg;
763 float varInit, varMin, varMax, prune, tau;
764
765 bool detectShadows;
766 uchar shadowVal;
767 };
768
769 #ifdef HAVE_OPENCL
770
ocl_apply(InputArray _image,OutputArray _fgmask,double learningRate)771 bool BackgroundSubtractorMOG2Impl::ocl_apply(InputArray _image, OutputArray _fgmask, double learningRate)
772 {
773 bool needToInitialize = nframes == 0 || learningRate >= 1 || _image.size() != frameSize || _image.type() != frameType;
774
775 if( needToInitialize )
776 initialize(_image.size(), _image.type());
777
778 ++nframes;
779 learningRate = learningRate >= 0 && nframes > 1 ? learningRate : 1./std::min( 2*nframes, history );
780 CV_Assert(learningRate >= 0);
781
782 _fgmask.create(_image.size(), CV_8U);
783 UMat fgmask = _fgmask.getUMat();
784
785 const double alpha1 = 1.0f - learningRate;
786
787 UMat frame = _image.getUMat();
788
789 float varMax = MAX(fVarMin, fVarMax);
790 float varMin = MIN(fVarMin, fVarMax);
791
792 int idxArg = 0;
793 idxArg = kernel_apply.set(idxArg, ocl::KernelArg::ReadOnly(frame));
794 idxArg = kernel_apply.set(idxArg, ocl::KernelArg::PtrReadWrite(u_bgmodelUsedModes));
795 idxArg = kernel_apply.set(idxArg, ocl::KernelArg::PtrReadWrite(u_weight));
796 idxArg = kernel_apply.set(idxArg, ocl::KernelArg::PtrReadWrite(u_mean));
797 idxArg = kernel_apply.set(idxArg, ocl::KernelArg::PtrReadWrite(u_variance));
798 idxArg = kernel_apply.set(idxArg, ocl::KernelArg::WriteOnlyNoSize(fgmask));
799
800 idxArg = kernel_apply.set(idxArg, (float)learningRate); //alphaT
801 idxArg = kernel_apply.set(idxArg, (float)alpha1);
802 idxArg = kernel_apply.set(idxArg, (float)(-learningRate*fCT)); //prune
803
804 idxArg = kernel_apply.set(idxArg, (float)varThreshold); //c_Tb
805 idxArg = kernel_apply.set(idxArg, backgroundRatio); //c_TB
806 idxArg = kernel_apply.set(idxArg, varThresholdGen); //c_Tg
807 idxArg = kernel_apply.set(idxArg, varMin);
808 idxArg = kernel_apply.set(idxArg, varMax);
809 idxArg = kernel_apply.set(idxArg, fVarInit);
810 idxArg = kernel_apply.set(idxArg, fTau);
811 if (bShadowDetection)
812 kernel_apply.set(idxArg, nShadowDetection);
813
814 size_t globalsize[] = {(size_t)frame.cols, (size_t)frame.rows, 1};
815 return kernel_apply.run(2, globalsize, NULL, true);
816 }
817
ocl_getBackgroundImage(OutputArray _backgroundImage) const818 bool BackgroundSubtractorMOG2Impl::ocl_getBackgroundImage(OutputArray _backgroundImage) const
819 {
820 _backgroundImage.create(frameSize, frameType);
821 UMat dst = _backgroundImage.getUMat();
822
823 int idxArg = 0;
824 idxArg = kernel_getBg.set(idxArg, ocl::KernelArg::PtrReadOnly(u_bgmodelUsedModes));
825 idxArg = kernel_getBg.set(idxArg, ocl::KernelArg::PtrReadOnly(u_weight));
826 idxArg = kernel_getBg.set(idxArg, ocl::KernelArg::PtrReadOnly(u_mean));
827 idxArg = kernel_getBg.set(idxArg, ocl::KernelArg::WriteOnly(dst));
828 kernel_getBg.set(idxArg, backgroundRatio);
829
830 size_t globalsize[2] = {(size_t)u_bgmodelUsedModes.cols, (size_t)u_bgmodelUsedModes.rows};
831
832 return kernel_getBg.run(2, globalsize, NULL, false);
833 }
834
create_ocl_apply_kernel()835 void BackgroundSubtractorMOG2Impl::create_ocl_apply_kernel()
836 {
837 int nchannels = CV_MAT_CN(frameType);
838 bool isFloat = CV_MAKETYPE(CV_32F,nchannels) == frameType;
839 String opts = format("-D CN=%d -D FL=%d -D NMIXTURES=%d%s", nchannels, isFloat, nmixtures, bShadowDetection ? " -D SHADOW_DETECT" : "");
840 kernel_apply.create("mog2_kernel", ocl::video::bgfg_mog2_oclsrc, opts);
841 }
842
843 #endif
844
apply(InputArray _image,OutputArray _fgmask,double learningRate)845 void BackgroundSubtractorMOG2Impl::apply(InputArray _image, OutputArray _fgmask, double learningRate)
846 {
847 CV_INSTRUMENT_REGION();
848
849 #ifdef HAVE_OPENCL
850 if (opencl_ON)
851 {
852 CV_OCL_RUN(_fgmask.isUMat(), ocl_apply(_image, _fgmask, learningRate))
853
854 opencl_ON = false;
855 nframes = 0;
856 }
857 #endif
858
859 bool needToInitialize = nframes == 0 || learningRate >= 1 || _image.size() != frameSize || _image.type() != frameType;
860
861 if( needToInitialize )
862 initialize(_image.size(), _image.type());
863
864 Mat image = _image.getMat();
865 _fgmask.create( image.size(), CV_8U );
866 Mat fgmask = _fgmask.getMat();
867
868 ++nframes;
869 learningRate = learningRate >= 0 && nframes > 1 ? learningRate : 1./std::min( 2*nframes, history );
870 CV_Assert(learningRate >= 0);
871
872 parallel_for_(Range(0, image.rows),
873 MOG2Invoker(image, fgmask,
874 bgmodel.ptr<GMM>(),
875 (float*)(bgmodel.ptr() + sizeof(GMM)*nmixtures*image.rows*image.cols),
876 bgmodelUsedModes.ptr(), nmixtures, (float)learningRate,
877 (float)varThreshold,
878 backgroundRatio, varThresholdGen,
879 fVarInit, fVarMin, fVarMax, float(-learningRate*fCT), fTau,
880 bShadowDetection, nShadowDetection),
881 image.total()/(double)(1 << 16));
882 }
883
884 template <typename T, int CN>
getBackgroundImage_intern(OutputArray backgroundImage) const885 void BackgroundSubtractorMOG2Impl::getBackgroundImage_intern(OutputArray backgroundImage) const
886 {
887 CV_INSTRUMENT_REGION();
888
889 Mat meanBackground(frameSize, frameType, Scalar::all(0));
890 int firstGaussianIdx = 0;
891 const GMM* gmm = bgmodel.ptr<GMM>();
892 const float* mean = reinterpret_cast<const float*>(gmm + frameSize.width*frameSize.height*nmixtures);
893 Vec<float,CN> meanVal(0.f);
894 for(int row=0; row<meanBackground.rows; row++)
895 {
896 for(int col=0; col<meanBackground.cols; col++)
897 {
898 int nmodes = bgmodelUsedModes.at<uchar>(row, col);
899 float totalWeight = 0.f;
900 for(int gaussianIdx = firstGaussianIdx; gaussianIdx < firstGaussianIdx + nmodes; gaussianIdx++)
901 {
902 GMM gaussian = gmm[gaussianIdx];
903 size_t meanPosition = gaussianIdx*CN;
904 for(int chn = 0; chn < CN; chn++)
905 {
906 meanVal(chn) += gaussian.weight * mean[meanPosition + chn];
907 }
908 totalWeight += gaussian.weight;
909
910 if(totalWeight > backgroundRatio)
911 break;
912 }
913 float invWeight = 0.f;
914 if (std::abs(totalWeight) > FLT_EPSILON) {
915 invWeight = 1.f/totalWeight;
916 }
917
918 meanBackground.at<Vec<T,CN> >(row, col) = Vec<T,CN>(meanVal * invWeight);
919 meanVal = 0.f;
920
921 firstGaussianIdx += nmixtures;
922 }
923 }
924 meanBackground.copyTo(backgroundImage);
925 }
926
getBackgroundImage(OutputArray backgroundImage) const927 void BackgroundSubtractorMOG2Impl::getBackgroundImage(OutputArray backgroundImage) const
928 {
929 CV_Assert(frameType == CV_8UC1 || frameType == CV_8UC3 || frameType == CV_32FC1 || frameType == CV_32FC3);
930
931 #ifdef HAVE_OPENCL
932 if (opencl_ON)
933 {
934 CV_OCL_RUN(opencl_ON, ocl_getBackgroundImage(backgroundImage))
935
936 opencl_ON = false;
937 }
938 #endif
939
940 switch(frameType)
941 {
942 case CV_8UC1:
943 getBackgroundImage_intern<uchar,1>(backgroundImage);
944 break;
945 case CV_8UC3:
946 getBackgroundImage_intern<uchar,3>(backgroundImage);
947 break;
948 case CV_32FC1:
949 getBackgroundImage_intern<float,1>(backgroundImage);
950 break;
951 case CV_32FC3:
952 getBackgroundImage_intern<float,3>(backgroundImage);
953 break;
954 }
955 }
956
createBackgroundSubtractorMOG2(int _history,double _varThreshold,bool _bShadowDetection)957 Ptr<BackgroundSubtractorMOG2> createBackgroundSubtractorMOG2(int _history, double _varThreshold,
958 bool _bShadowDetection)
959 {
960 return makePtr<BackgroundSubtractorMOG2Impl>(_history, (float)_varThreshold, _bShadowDetection);
961 }
962
963 }
964
965 /* End of file. */
966