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 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 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 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 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 233 virtual int getHistory() const { return history; } 234 virtual void setHistory(int _nframes) { history = _nframes; } 235 236 virtual int getNMixtures() const { return nmixtures; } 237 virtual void setNMixtures(int nmix) { nmixtures = nmix; } 238 239 virtual double getBackgroundRatio() const { return backgroundRatio; } 240 virtual void setBackgroundRatio(double _backgroundRatio) { backgroundRatio = (float)_backgroundRatio; } 241 242 virtual double getVarThreshold() const { return varThreshold; } 243 virtual void setVarThreshold(double _varThreshold) { varThreshold = _varThreshold; } 244 245 virtual double getVarThresholdGen() const { return varThresholdGen; } 246 virtual void setVarThresholdGen(double _varThresholdGen) { varThresholdGen = (float)_varThresholdGen; } 247 248 virtual double getVarInit() const { return fVarInit; } 249 virtual void setVarInit(double varInit) { fVarInit = (float)varInit; } 250 251 virtual double getVarMin() const { return fVarMin; } 252 virtual void setVarMin(double varMin) { fVarMin = (float)varMin; } 253 254 virtual double getVarMax() const { return fVarMax; } 255 virtual void setVarMax(double varMax) { fVarMax = (float)varMax; } 256 257 virtual double getComplexityReductionThreshold() const { return fCT; } 258 virtual void setComplexityReductionThreshold(double ct) { fCT = (float)ct; } 259 260 virtual bool getDetectShadows() const { return bShadowDetection; } 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 273 virtual int getShadowValue() const { return nShadowDetection; } 274 virtual void setShadowValue(int value) { nShadowDetection = (uchar)value; } 275 276 virtual double getShadowThreshold() const { return fTau; } 277 virtual void setShadowThreshold(double value) { fTau = (float)value; } 278 279 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 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 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: 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 552 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 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 794 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 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 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 857 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 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
