modules/legacy/src/bgfg_gaussmix.cpp

   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 //                        Intel License Agreement
  11 //
  12 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  13 // Third party copyrights are property of their respective owners.
  14 //
  15 // Redistribution and use in source and binary forms, with or without modification,
  16 // are permitted provided that the following conditions are met:
  17 //
  18 //   * Redistribution's of source code must retain the above copyright notice,
  19 //     this list of conditions and the following disclaimer.
  20 //
  21 //   * Redistribution's in binary form must reproduce the above copyright notice,
  22 //     this list of conditions and the following disclaimer in the documentation
  23 //     and/or other materials provided with the distribution.
  24 //
  25 //   * The name of Intel Corporation may not be used to endorse or promote products
  26 //     derived from this software without specific prior written permission.
  27 //
  28 // This software is provided by the copyright holders and contributors "as is" and
  29 // any express or implied warranties, including, but not limited to, the implied
  30 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  31 // In no event shall the Intel Corporation or contributors be liable for any direct,
  32 // indirect, incidental, special, exemplary, or consequential damages
  33 // (including, but not limited to, procurement of substitute goods or services;
  34 // loss of use, data, or profits; or business interruption) however caused
  35 // and on any theory of liability, whether in contract, strict liability,
  36 // or tort (including negligence or otherwise) arising in any way out of
  37 // the use of this software, even if advised of the possibility of such damage.
  38 //
  39 //M*/
  40
  41 #include "precomp.hpp"
  42
  43 /////////////////////////////////////// MOG model //////////////////////////////////////////
  44
  45 static void CV_CDECL
  46 icvReleaseGaussianBGModel( CvGaussBGModel** bg_model )
  47 {
  48     if( !bg_model )
  49         CV_Error( CV_StsNullPtr, "" );
  50
  51     if( *bg_model )
  52     {
  53         delete (cv::Ptr<cv::BackgroundSubtractor>*)((*bg_model)->mog);
  54         cvReleaseImage( &(*bg_model)->background );
  55         cvReleaseImage( &(*bg_model)->foreground );
  56         memset( *bg_model, 0, sizeof(**bg_model) );
  57         delete *bg_model;
  58         *bg_model = 0;
  59     }
  60 }
  61
  62
  63 static int CV_CDECL
  64 icvUpdateGaussianBGModel( IplImage* curr_frame, CvGaussBGModel*  bg_model, double learningRate )
  65 {
  66     cv::Mat image = cv::cvarrToMat(curr_frame), mask = cv::cvarrToMat(bg_model->foreground);
  67
  68     cv::Ptr<cv::BackgroundSubtractor>* mog = (cv::Ptr<cv::BackgroundSubtractor>*)(bg_model->mog);
  69     CV_Assert(mog != 0);
  70
  71     (*mog)->apply(image, mask, learningRate);
  72     bg_model->countFrames++;
  73
  74     return 0;
  75 }
  76
  77 CV_IMPL CvBGStatModel*
  78 cvCreateGaussianBGModel( IplImage* first_frame, CvGaussBGStatModelParams* parameters )
  79 {
  80     CvGaussBGStatModelParams params;
  81
  82     CV_Assert( CV_IS_IMAGE(first_frame) );
  83
  84     //init parameters
  85     if( parameters == NULL )
  86     {                        // These constants are defined in cvaux/include/cvaux.h
  87         params.win_size      = CV_BGFG_MOG_WINDOW_SIZE;
  88         params.bg_threshold  = CV_BGFG_MOG_BACKGROUND_THRESHOLD;
  89
  90         params.std_threshold = CV_BGFG_MOG_STD_THRESHOLD;
  91         params.weight_init   = CV_BGFG_MOG_WEIGHT_INIT;
  92
  93         params.variance_init = CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT;
  94         params.minArea       = CV_BGFG_MOG_MINAREA;
  95         params.n_gauss       = CV_BGFG_MOG_NGAUSSIANS;
  96     }
  97     else
  98         params = *parameters;
  99
 100     CvGaussBGModel* bg_model = new CvGaussBGModel;
 101     memset( bg_model, 0, sizeof(*bg_model) );
 102     bg_model->type = CV_BG_MODEL_MOG;
 103     bg_model->release = (CvReleaseBGStatModel)icvReleaseGaussianBGModel;
 104     bg_model->update = (CvUpdateBGStatModel)icvUpdateGaussianBGModel;
 105
 106     bg_model->params = params;
 107
 108     cv::Ptr<cv::BackgroundSubtractor> mog = cv::createBackgroundSubtractorMOG(params.win_size, params.n_gauss,
 109                                                                               params.bg_threshold);
 110     cv::Ptr<cv::BackgroundSubtractor>* pmog = new cv::Ptr<cv::BackgroundSubtractor>;
 111     *pmog = mog;
 112     bg_model->mog = pmog;
 113
 114     CvSize sz = cvGetSize(first_frame);
 115     bg_model->background = cvCreateImage(sz, IPL_DEPTH_8U, first_frame->nChannels);
 116     bg_model->foreground = cvCreateImage(sz, IPL_DEPTH_8U, 1);
 117
 118     bg_model->countFrames = 0;
 119
 120     icvUpdateGaussianBGModel( first_frame, bg_model, 1 );
 121
 122     return (CvBGStatModel*)bg_model;
 123 }
 124
 125
 126 //////////////////////////////////////////// MOG2 //////////////////////////////////////////////
 127
 128 /*M///////////////////////////////////////////////////////////////////////////////////////
 129  //
 130  //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
 131  //
 132  //  By downloading, copying, installing or using the software you agree to this license.
 133  //  If you do not agree to this license, do not download, install,
 134  //  copy or use the software.
 135  //
 136  //
 137  //                        Intel License Agreement
 138  //
 139  // Copyright (C) 2000, Intel Corporation, all rights reserved.
 140  // Third party copyrights are property of their respective owners.
 141  //
 142  // Redistribution and use in source and binary forms, with or without modification,
 143  // are permitted provided that the following conditions are met:
 144  //
 145  //   * Redistribution's of source code must retain the above copyright notice,
 146  //     this list of conditions and the following disclaimer.
 147  //
 148  //   * Redistribution's in binary form must reproduce the above copyright notice,
 149  //     this list of conditions and the following disclaimer in the documentation
 150  //     and/or other materials provided with the distribution.
 151  //
 152  //   * The name of Intel Corporation may not be used to endorse or promote products
 153  //     derived from this software without specific prior written permission.
 154  //
 155  // This software is provided by the copyright holders and contributors "as is" and
 156  // any express or implied warranties, including, but not limited to, the implied
 157  // warranties of merchantability and fitness for a particular purpose are disclaimed.
 158  // In no event shall the Intel Corporation or contributors be liable for any direct,
 159  // indirect, incidental, special, exemplary, or consequential damages
 160  // (including, but not limited to, procurement of substitute goods or services;
 161  // loss of use, data, or profits; or business interruption) however caused
 162  // and on any theory of liability, whether in contract, strict liability,
 163  // or tort (including negligence or otherwise) arising in any way out of
 164  // the use of this software, even if advised of the possibility of such damage.
 165  //
 166  //M*/
 167
 168 /*//Implementation of the Gaussian mixture model background subtraction from:
 169  //
 170  //"Improved adaptive Gausian mixture model for background subtraction"
 171  //Z.Zivkovic
 172  //International Conference Pattern Recognition, UK, August, 2004
 173  //http://www.zoranz.net/Publications/zivkovic2004ICPR.pdf
 174  //The code is very fast and performs also shadow detection.
 175  //Number of Gausssian components is adapted per pixel.
 176  //
 177  // and
 178  //
 179  //"Efficient Adaptive Density Estimapion per Image Pixel for the Task of Background Subtraction"
 180  //Z.Zivkovic, F. van der Heijden
 181  //Pattern Recognition Letters, vol. 27, no. 7, pages 773-780, 2006.
 182  //
 183  //The algorithm similar to the standard Stauffer&Grimson algorithm with
 184  //additional selection of the number of the Gaussian components based on:
 185  //
 186  //"Recursive unsupervised learning of finite mixture models "
 187  //Z.Zivkovic, F.van der Heijden
 188  //IEEE Trans. on Pattern Analysis and Machine Intelligence, vol.26, no.5, pages 651-656, 2004
 189  //http://www.zoranz.net/Publications/zivkovic2004PAMI.pdf
 190  //
 191  //
 192  //Example usage with as cpp class
 193  // BackgroundSubtractorMOG2 bg_model;
 194  //For each new image the model is updates using:
 195  // bg_model(img, fgmask);
 196  //
 197  //Example usage as part of the CvBGStatModel:
 198  // CvBGStatModel* bg_model = cvCreateGaussianBGModel2( first_frame );
 199  //
 200  // //update for each frame
 201  // cvUpdateBGStatModel( tmp_frame, bg_model );//segmentation result is in bg_model->foreground
 202  //
 203  // //release at the program termination
 204  // cvReleaseBGStatModel( &bg_model );
 205  //
 206  //Author: Z.Zivkovic, www.zoranz.net
 207  //Date: 7-April-2011, Version:1.0
 208  ///////////*/
 209
 210 #include "precomp.hpp"
 211
 212
 213 /*
 214  Interface of Gaussian mixture algorithm from:
 215
 216  "Improved adaptive Gausian mixture model for background subtraction"
 217  Z.Zivkovic
 218  International Conference Pattern Recognition, UK, August, 2004
 219  http://www.zoranz.net/Publications/zivkovic2004ICPR.pdf
 220
 221  Advantages:
 222  -fast - number of Gausssian components is constantly adapted per pixel.
 223  -performs also shadow detection (see bgfg_segm_test.cpp example)
 224
 225  */
 226
 227
 228 #define CV_BG_MODEL_MOG2            3                 /* "Mixture of Gaussians 2".  */
 229
 230
 231 /* default parameters of gaussian background detection algorithm */
 232 #define CV_BGFG_MOG2_STD_THRESHOLD            4.0f     /* lambda=2.5 is 99% */
 233 #define CV_BGFG_MOG2_WINDOW_SIZE              500      /* Learning rate; alpha = 1/CV_GBG_WINDOW_SIZE */
 234 #define CV_BGFG_MOG2_BACKGROUND_THRESHOLD     0.9f     /* threshold sum of weights for background test */
 235 #define CV_BGFG_MOG2_STD_THRESHOLD_GENERATE   3.0f     /* lambda=2.5 is 99% */
 236 #define CV_BGFG_MOG2_NGAUSSIANS               5        /* = K = number of Gaussians in mixture */
 237 #define CV_BGFG_MOG2_VAR_INIT                 15.0f    /* initial variance for new components*/
 238 #define CV_BGFG_MOG2_VAR_MIN                    4.0f
 239 #define CV_BGFG_MOG2_VAR_MAX                      5*CV_BGFG_MOG2_VAR_INIT
 240 #define CV_BGFG_MOG2_MINAREA                  15.0f    /* for postfiltering */
 241
 242 /* additional parameters */
 243 #define CV_BGFG_MOG2_CT                               0.05f     /* complexity reduction prior constant 0 - no reduction of number of components*/
 244 #define CV_BGFG_MOG2_SHADOW_VALUE             127       /* value to use in the segmentation mask for shadows, sot 0 not to do shadow detection*/
 245 #define CV_BGFG_MOG2_SHADOW_TAU               0.5f      /* Tau - shadow threshold, see the paper for explanation*/
 246
 247 typedef struct CvGaussBGStatModel2Params
 248 {
 249     //image info
 250     int nWidth;
 251     int nHeight;
 252     int nND;//number of data dimensions (image channels)
 253
 254     bool bPostFiltering;//defult 1 - do postfiltering - will make shadow detection results also give value 255
 255     double  minArea; // for postfiltering
 256
 257     bool bInit;//default 1, faster updates at start
 258
 259     /////////////////////////
 260     //very important parameters - things you will change
 261     ////////////////////////
 262     float fAlphaT;
 263     //alpha - speed of update - if the time interval you want to average over is T
 264     //set alpha=1/T. It is also usefull at start to make T slowly increase
 265     //from 1 until the desired T
 266     float fTb;
 267     //Tb - threshold on the squared Mahalan. dist. to decide if it is well described
 268     //by the background model or not. Related to Cthr from the paper.
 269     //This does not influence the update of the background. A typical value could be 4 sigma
 270     //and that is Tb=4*4=16;
 271
 272     /////////////////////////
 273     //less important parameters - things you might change but be carefull
 274     ////////////////////////
 275     float fTg;
 276     //Tg - threshold on the squared Mahalan. dist. to decide
 277     //when a sample is close to the existing components. If it is not close
 278     //to any a new component will be generated. I use 3 sigma => Tg=3*3=9.
 279     //Smaller Tg leads to more generated components and higher Tg might make
 280     //lead to small number of components but they can grow too large
 281     float fTB;//1-cf from the paper
 282     //TB - threshold when the component becomes significant enough to be included into
 283     //the background model. It is the TB=1-cf from the paper. So I use cf=0.1 => TB=0.
 284     //For alpha=0.001 it means that the mode should exist for approximately 105 frames before
 285     //it is considered foreground
 286     float fVarInit;
 287     float fVarMax;
 288     float fVarMin;
 289     //initial standard deviation  for the newly generated components.
 290     //It will will influence the speed of adaptation. A good guess should be made.
 291     //A simple way is to estimate the typical standard deviation from the images.
 292     //I used here 10 as a reasonable value
 293     float fCT;//CT - complexity reduction prior
 294     //this is related to the number of samples needed to accept that a component
 295     //actually exists. We use CT=0.05 of all the samples. By setting CT=0 you get
 296     //the standard Stauffer&Grimson algorithm (maybe not exact but very similar)
 297
 298     //even less important parameters
 299     int nM;//max number of modes - const - 4 is usually enough
 300
 301     //shadow detection parameters
 302     bool bShadowDetection;//default 1 - do shadow detection
 303     unsigned char nShadowDetection;//do shadow detection - insert this value as the detection result
 304     float fTau;
 305     // Tau - shadow threshold. The shadow is detected if the pixel is darker
 306     //version of the background. Tau is a threshold on how much darker the shadow can be.
 307     //Tau= 0.5 means that if pixel is more than 2 times darker then it is not shadow
 308     //See: Prati,Mikic,Trivedi,Cucchiarra,"Detecting Moving Shadows...",IEEE PAMI,2003.
 309 } CvGaussBGStatModel2Params;
 310
 311 #define CV_BGFG_MOG2_NDMAX 3
 312
 313 typedef struct CvPBGMMGaussian
 314 {
 315     float weight;
 316     float mean[CV_BGFG_MOG2_NDMAX];
 317     float variance;
 318 }CvPBGMMGaussian;
 319
 320 typedef struct CvGaussBGStatModel2Data
 321 {
 322     CvPBGMMGaussian* rGMM; //array for the mixture of Gaussians
 323     unsigned char* rnUsedModes;//number of Gaussian components per pixel (maximum 255)
 324 } CvGaussBGStatModel2Data;
 325
 326
 327 /*
 328  //only foreground image is updated
 329  //no filtering included
 330  typedef struct CvGaussBGModel2
 331  {
 332  CV_BG_STAT_MODEL_FIELDS();
 333  CvGaussBGStatModel2Params params;
 334  CvGaussBGStatModel2Data   data;
 335  int                       countFrames;
 336  } CvGaussBGModel2;
 337
 338  CVAPI(CvBGStatModel*) cvCreateGaussianBGModel2( IplImage* first_frame,
 339  CvGaussBGStatModel2Params* params CV_DEFAULT(NULL) );
 340  */
 341
 342 //shadow detection performed per pixel
 343 // should work for rgb data, could be usefull for gray scale and depth data as well
 344 //  See: Prati,Mikic,Trivedi,Cucchiarra,"Detecting Moving Shadows...",IEEE PAMI,2003.
 345 CV_INLINE int _icvRemoveShadowGMM(float* data, int nD,
 346                                   unsigned char nModes,
 347                                   CvPBGMMGaussian* pGMM,
 348                                   float m_fTb,
 349                                   float m_fTB,
 350                                   float m_fTau)
 351 {
 352     float tWeight = 0;
 353     float numerator, denominator;
 354     // check all the components  marked as background:
 355     for (int iModes=0;iModes<nModes;iModes++)
 356     {
 357
 358         CvPBGMMGaussian g=pGMM[iModes];
 359
 360         numerator = 0.0f;
 361         denominator = 0.0f;
 362         for (int iD=0;iD<nD;iD++)
 363         {
 364             numerator   += data[iD]  * g.mean[iD];
 365             denominator += g.mean[iD]* g.mean[iD];
 366         }
 367
 368         // no division by zero allowed
 369         if (denominator == 0)
 370         {
 371             return 0;
 372         };
 373         float a = numerator / denominator;
 374
 375         // if tau < a < 1 then also check the color distortion
 376         if ((a <= 1) && (a >= m_fTau))
 377         {
 378
 379             float dist2a=0.0f;
 380
 381             for (int iD=0;iD<nD;iD++)
 382             {
 383                 float dD= a*g.mean[iD] - data[iD];
 384                 dist2a += (dD*dD);
 385             }
 386
 387             if (dist2a<m_fTb*g.variance*a*a)
 388             {
 389                 return 2;
 390             }
 391         };
 392
 393         tWeight += g.weight;
 394         if (tWeight > m_fTB)
 395         {
 396             return 0;
 397         };
 398     };
 399     return 0;
 400 }
 401
 402 //update GMM - the base update function performed per pixel
 403 //
 404 //"Efficient Adaptive Density Estimapion per Image Pixel for the Task of Background Subtraction"
 405 //Z.Zivkovic, F. van der Heijden
 406 //Pattern Recognition Letters, vol. 27, no. 7, pages 773-780, 2006.
 407 //
 408 //The algorithm similar to the standard Stauffer&Grimson algorithm with
 409 //additional selection of the number of the Gaussian components based on:
 410 //
 411 //"Recursive unsupervised learning of finite mixture models "
 412 //Z.Zivkovic, F.van der Heijden
 413 //IEEE Trans. on Pattern Analysis and Machine Intelligence, vol.26, no.5, pages 651-656, 2004
 414 //http://www.zoranz.net/Publications/zivkovic2004PAMI.pdf
 415
 416 #if defined(__GNUC__) && (__GNUC__ == 4) && (__GNUC_MINOR__ == 8)
 417 # pragma GCC diagnostic push
 418 # pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
 419 #endif
 420
 421 CV_INLINE int _icvUpdateGMM(float* data, int nD,
 422                             unsigned char* pModesUsed,
 423                             CvPBGMMGaussian* pGMM,
 424                             int m_nM,
 425                             float m_fAlphaT,
 426                             float m_fTb,
 427                             float m_fTB,
 428                             float m_fTg,
 429                             float m_fVarInit,
 430                             float m_fVarMax,
 431                             float m_fVarMin,
 432                             float m_fPrune)
 433 {
 434     //calculate distances to the modes (+ sort)
 435     //here we need to go in descending order!!!
 436     bool bBackground=0;//return value -> true - the pixel classified as background
 437
 438     //internal:
 439     bool bFitsPDF=0;//if it remains zero a new GMM mode will be added
 440     float m_fOneMinAlpha=1-m_fAlphaT;
 441     unsigned char nModes=*pModesUsed;//current number of modes in GMM
 442     float totalWeight=0.0f;
 443
 444     //////
 445     //go through all modes
 446     int iMode=0;
 447     CvPBGMMGaussian* pGauss=pGMM;
 448     for (;iMode<nModes;iMode++,pGauss++)
 449     {
 450         float weight = pGauss->weight;//need only weight if fit is found
 451         weight=m_fOneMinAlpha*weight+m_fPrune;
 452
 453         ////
 454         //fit not found yet
 455         if (!bFitsPDF)
 456         {
 457             //check if it belongs to some of the remaining modes
 458             float var=pGauss->variance;
 459
 460             //calculate difference and distance
 461             float dist2=0.0f;
 462 #if (CV_BGFG_MOG2_NDMAX==1)
 463             float dData=pGauss->mean[0]-data[0];
 464             dist2=dData*dData;
 465 #else
 466             float dData[CV_BGFG_MOG2_NDMAX];
 467
 468             for (int iD=0;iD<nD;iD++)
 469             {
 470                 dData[iD]=pGauss->mean[iD]-data[iD];
 471                 dist2+=dData[iD]*dData[iD];
 472             }
 473 #endif
 474             //background? - m_fTb - usually larger than m_fTg
 475             if ((totalWeight<m_fTB)&&(dist2<m_fTb*var))
 476                 bBackground=1;
 477
 478             //check fit
 479             if (dist2<m_fTg*var)
 480             {
 481                 /////
 482                 //belongs to the mode - bFitsPDF becomes 1
 483                 bFitsPDF=1;
 484
 485                 //update distribution
 486
 487                 //update weight
 488                 weight+=m_fAlphaT;
 489
 490                 float k = m_fAlphaT/weight;
 491
 492                 //update mean
 493 #if (CV_BGFG_MOG2_NDMAX==1)
 494                 pGauss->mean[0]-=k*dData;
 495 #else
 496                 for (int iD=0;iD<nD;iD++)
 497                 {
 498                     pGauss->mean[iD]-=k*dData[iD];
 499                 }
 500 #endif
 501
 502                 //update variance
 503                 float varnew = var + k*(dist2-var);
 504                 //limit the variance
 505                 pGauss->variance = MIN(m_fVarMax,MAX(varnew,m_fVarMin));
 506
 507                 //sort
 508                 //all other weights are at the same place and
 509                 //only the matched (iModes) is higher -> just find the new place for it
 510                 for (int iLocal = iMode;iLocal>0;iLocal--)
 511                 {
 512                     //check one up
 513                     if (weight < (pGMM[iLocal-1].weight))
 514                     {
 515                         break;
 516                     }
 517                     else
 518                     {
 519                         //swap one up
 520                         CvPBGMMGaussian temp = pGMM[iLocal];
 521                         pGMM[iLocal] = pGMM[iLocal-1];
 522                         pGMM[iLocal-1] = temp;
 523                         pGauss--;
 524                     }
 525                 }
 526                 //belongs to the mode - bFitsPDF becomes 1
 527                 /////
 528             }
 529         }//!bFitsPDF)
 530
 531         //check prune
 532         if (weight<-m_fPrune)
 533         {
 534             weight=0.0;
 535             nModes--;
 536         }
 537
 538         pGauss->weight=weight;//update weight by the calculated value
 539         totalWeight+=weight;
 540     }
 541     //go through all modes
 542     //////
 543
 544     //renormalize weights
 545     for (iMode = 0; iMode < nModes; iMode++)
 546     {
 547         pGMM[iMode].weight = pGMM[iMode].weight/totalWeight;
 548     }
 549
 550     //make new mode if needed and exit
 551     if (!bFitsPDF)
 552     {
 553         if (nModes==m_nM)
 554         {
 555             //replace the weakest
 556             pGauss=pGMM+m_nM-1;
 557         }
 558         else
 559         {
 560             //add a new one
 561             pGauss=pGMM+nModes;
 562             nModes++;
 563         }
 564
 565         if (nModes==1)
 566         {
 567             pGauss->weight=1;
 568         }
 569         else
 570         {
 571             pGauss->weight=m_fAlphaT;
 572
 573             //renormalize all weights
 574             for (iMode = 0; iMode < nModes-1; iMode++)
 575             {
 576                 pGMM[iMode].weight *=m_fOneMinAlpha;
 577             }
 578         }
 579
 580         //init
 581         memcpy(pGauss->mean,data,nD*sizeof(float));
 582         pGauss->variance=m_fVarInit;
 583
 584         //sort
 585         //find the new place for it
 586         for (int iLocal = nModes-1;iLocal>0;iLocal--)
 587         {
 588             //check one up
 589             if (m_fAlphaT < (pGMM[iLocal-1].weight))
 590             {
 591                 break;
 592             }
 593             else
 594             {
 595                 //swap one up
 596                 CvPBGMMGaussian temp = pGMM[iLocal];
 597                 pGMM[iLocal] = pGMM[iLocal-1];
 598                 pGMM[iLocal-1] = temp;
 599             }
 600         }
 601     }
 602
 603     //set the number of modes
 604     *pModesUsed=nModes;
 605
 606     return bBackground;
 607 }
 608
 609 #if defined(__GNUC__) && (__GNUC__ == 4) && (__GNUC_MINOR__ == 8)
 610 # pragma GCC diagnostic pop
 611 #endif
 612
 613 // a bit more efficient implementation for common case of 3 channel (rgb) images
 614 CV_INLINE int _icvUpdateGMM_C3(float r,float g, float b,
 615                                unsigned char* pModesUsed,
 616                                CvPBGMMGaussian* pGMM,
 617                                int m_nM,
 618                                float m_fAlphaT,
 619                                float m_fTb,
 620                                float m_fTB,
 621                                float m_fTg,
 622                                float m_fVarInit,
 623                                float m_fVarMax,
 624                                float m_fVarMin,
 625                                float m_fPrune)
 626 {
 627     //calculate distances to the modes (+ sort)
 628     //here we need to go in descending order!!!
 629     bool bBackground=0;//return value -> true - the pixel classified as background
 630
 631     //internal:
 632     bool bFitsPDF=0;//if it remains zero a new GMM mode will be added
 633     float m_fOneMinAlpha=1-m_fAlphaT;
 634     unsigned char nModes=*pModesUsed;//current number of modes in GMM
 635     float totalWeight=0.0f;
 636
 637     //////
 638     //go through all modes
 639     int iMode=0;
 640     CvPBGMMGaussian* pGauss=pGMM;
 641     for (;iMode<nModes;iMode++,pGauss++)
 642     {
 643         float weight = pGauss->weight;//need only weight if fit is found
 644         weight=m_fOneMinAlpha*weight+m_fPrune;
 645
 646         ////
 647         //fit not found yet
 648         if (!bFitsPDF)
 649         {
 650             //check if it belongs to some of the remaining modes
 651             float var=pGauss->variance;
 652
 653             //calculate difference and distance
 654             float muR = pGauss->mean[0];
 655             float muG = pGauss->mean[1];
 656             float muB = pGauss->mean[2];
 657
 658             float dR=muR - r;
 659             float dG=muG - g;
 660             float dB=muB - b;
 661
 662             float dist2=(dR*dR+dG*dG+dB*dB);
 663
 664             //background? - m_fTb - usually larger than m_fTg
 665             if ((totalWeight<m_fTB)&&(dist2<m_fTb*var))
 666                 bBackground=1;
 667
 668             //check fit
 669             if (dist2<m_fTg*var)
 670             {
 671                 /////
 672                 //belongs to the mode - bFitsPDF becomes 1
 673                 bFitsPDF=1;
 674
 675                 //update distribution
 676
 677                 //update weight
 678                 weight+=m_fAlphaT;
 679
 680                 float k = m_fAlphaT/weight;
 681
 682                 //update mean
 683                 pGauss->mean[0] = muR - k*(dR);
 684                 pGauss->mean[1] = muG - k*(dG);
 685                 pGauss->mean[2] = muB - k*(dB);
 686
 687                 //update variance
 688                 float varnew = var + k*(dist2-var);
 689                 //limit the variance
 690                 pGauss->variance = MIN(m_fVarMax,MAX(varnew,m_fVarMin));
 691
 692                 //sort
 693                 //all other weights are at the same place and
 694                 //only the matched (iModes) is higher -> just find the new place for it
 695                 for (int iLocal = iMode;iLocal>0;iLocal--)
 696                 {
 697                     //check one up
 698                     if (weight < (pGMM[iLocal-1].weight))
 699                     {
 700                         break;
 701                     }
 702                     else
 703                     {
 704                         //swap one up
 705                         CvPBGMMGaussian temp = pGMM[iLocal];
 706                         pGMM[iLocal] = pGMM[iLocal-1];
 707                         pGMM[iLocal-1] = temp;
 708                         pGauss--;
 709                     }
 710                 }
 711                 //belongs to the mode - bFitsPDF becomes 1
 712                 /////
 713             }
 714
 715         }//!bFitsPDF)
 716
 717         //check prunning
 718         if (weight<-m_fPrune)
 719         {
 720             weight=0.0;
 721             nModes--;
 722         }
 723
 724         pGauss->weight=weight;
 725         totalWeight+=weight;
 726     }
 727     //go through all modes
 728     //////
 729
 730     //renormalize weights
 731     for (iMode = 0; iMode < nModes; iMode++)
 732     {
 733         pGMM[iMode].weight = pGMM[iMode].weight/totalWeight;
 734     }
 735
 736     //make new mode if needed and exit
 737     if (!bFitsPDF)
 738     {
 739         if (nModes==m_nM)
 740         {
 741             //replace the weakest
 742             pGauss=pGMM+m_nM-1;
 743         }
 744         else
 745         {
 746             //add a new one
 747             pGauss=pGMM+nModes;
 748             nModes++;
 749         }
 750
 751         if (nModes==1)
 752         {
 753             pGauss->weight=1;
 754         }
 755         else
 756         {
 757             pGauss->weight=m_fAlphaT;
 758
 759             //renormalize all weights
 760             for (iMode = 0; iMode < nModes-1; iMode++)
 761             {
 762                 pGMM[iMode].weight *=m_fOneMinAlpha;
 763             }
 764         }
 765
 766         //init
 767         pGauss->mean[0]=r;
 768         pGauss->mean[1]=g;
 769         pGauss->mean[2]=b;
 770
 771         pGauss->variance=m_fVarInit;
 772
 773         //sort
 774         //find the new place for it
 775         for (int iLocal = nModes-1;iLocal>0;iLocal--)
 776         {
 777             //check one up
 778             if (m_fAlphaT < (pGMM[iLocal-1].weight))
 779             {
 780                 break;
 781             }
 782             else
 783             {
 784                 //swap one up
 785                 CvPBGMMGaussian temp = pGMM[iLocal];
 786                 pGMM[iLocal] = pGMM[iLocal-1];
 787                 pGMM[iLocal-1] = temp;
 788             }
 789         }
 790     }
 791
 792     //set the number of modes
 793     *pModesUsed=nModes;
 794
 795     return bBackground;
 796 }
 797
 798 //the main function to update the background model
 799 static void icvUpdatePixelBackgroundGMM2( const CvArr* srcarr, CvArr* dstarr ,
 800                                   CvPBGMMGaussian *pGMM,
 801                                   unsigned char *pUsedModes,
 802                                   //CvGaussBGStatModel2Params* pGMMPar,
 803                                   int nM,
 804                                   float fTb,
 805                                   float fTB,
 806                                   float fTg,
 807                                   float fVarInit,
 808                                   float fVarMax,
 809                                   float fVarMin,
 810                                   float fCT,
 811                                   float fTau,
 812                                   bool bShadowDetection,
 813                                   unsigned char  nShadowDetection,
 814                                   float alpha)
 815 {
 816     CvMat sstub, *src = cvGetMat(srcarr, &sstub);
 817     CvMat dstub, *dst = cvGetMat(dstarr, &dstub);
 818     CvSize size = cvGetMatSize(src);
 819     int nD=CV_MAT_CN(src->type);
 820
 821     //reshape if possible
 822     if( CV_IS_MAT_CONT(src->type & dst->type) )
 823     {
 824         size.width *= size.height;
 825         size.height = 1;
 826     }
 827
 828     int x, y;
 829     float data[CV_BGFG_MOG2_NDMAX];
 830     float prune=-alpha*fCT;
 831
 832     //general nD
 833
 834     if (nD!=3)
 835     {
 836         switch (CV_MAT_DEPTH(src->type))
 837         {
 838             case CV_8U:
 839                 for( y = 0; y < size.height; y++ )
 840                 {
 841                     uchar* sptr = src->data.ptr + src->step*y;
 842                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
 843                     for( x = 0; x < size.width; x++,
 844                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
 845                     {
 846                         //convert data
 847                         for (int iD=0;iD<nD;iD++) data[iD]=float(sptr[iD]);
 848                         //update GMM model
 849                         int result = _icvUpdateGMM(data,nD,pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
 850                         //detect shadows in the foreground
 851                         if (bShadowDetection)
 852                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
 853                         //generate output
 854                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
 855                     }
 856                 }
 857                 break;
 858             case CV_16S:
 859                 for( y = 0; y < size.height; y++ )
 860                 {
 861                     short* sptr = src->data.s + src->step*y;
 862                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
 863                     for( x = 0; x < size.width; x++,
 864                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
 865                     {
 866                         //convert data
 867                         for (int iD=0;iD<nD;iD++) data[iD]=float(sptr[iD]);
 868                         //update GMM model
 869                         int result = _icvUpdateGMM(data,nD,pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
 870                         //detect shadows in the foreground
 871                         if (bShadowDetection)
 872                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
 873                         //generate output
 874                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
 875                     }
 876                 }
 877                 break;
 878             case CV_16U:
 879                 for( y = 0; y < size.height; y++ )
 880                 {
 881                     unsigned short* sptr = (unsigned short*) (src->data.s + src->step*y);
 882                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
 883                     for( x = 0; x < size.width; x++,
 884                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
 885                     {
 886                         //convert data
 887                         for (int iD=0;iD<nD;iD++) data[iD]=float(sptr[iD]);
 888                         //update GMM model
 889                         int result = _icvUpdateGMM(data,nD,pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
 890                         //detect shadows in the foreground
 891                         if (bShadowDetection)
 892                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
 893                         //generate output
 894                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
 895                     }
 896                 }
 897                 break;
 898             case CV_32S:
 899                 for( y = 0; y < size.height; y++ )
 900                 {
 901                     int* sptr = src->data.i + src->step*y;
 902                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
 903                     for( x = 0; x < size.width; x++,
 904                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
 905                     {
 906                         //convert data
 907                         for (int iD=0;iD<nD;iD++) data[iD]=float(sptr[iD]);
 908                         //update GMM model
 909                         int result = _icvUpdateGMM(data,nD,pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
 910                         //detect shadows in the foreground
 911                         if (bShadowDetection)
 912                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
 913                         //generate output
 914                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
 915                     }
 916                 }
 917                 break;
 918             case CV_32F:
 919                 for( y = 0; y < size.height; y++ )
 920                 {
 921                     float* sptr = src->data.fl + src->step*y;
 922                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
 923                     for( x = 0; x < size.width; x++,
 924                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
 925                     {
 926                         //update GMM model
 927                         int result = _icvUpdateGMM(sptr,nD,pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
 928                         //detect shadows in the foreground
 929                         if (bShadowDetection)
 930                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
 931                         //generate output
 932                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
 933                     }
 934                 }
 935                 break;
 936             case CV_64F:
 937                 for( y = 0; y < size.height; y++ )
 938                 {
 939                     double* sptr = src->data.db + src->step*y;
 940                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
 941                     for( x = 0; x < size.width; x++,
 942                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
 943                     {
 944                         //convert data
 945                         for (int iD=0;iD<nD;iD++) data[iD]=float(sptr[iD]);
 946                         //update GMM model
 947                         int result = _icvUpdateGMM(data,nD,pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
 948                         //detect shadows in the foreground
 949                         if (bShadowDetection)
 950                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
 951                         //generate output
 952                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
 953                     }
 954                 }
 955                 break;
 956         }
 957     }else ///if (nD==3) - a bit faster
 958     {
 959         switch (CV_MAT_DEPTH(src->type))
 960         {
 961             case CV_8U:
 962                 for( y = 0; y < size.height; y++ )
 963                 {
 964                     uchar* sptr = src->data.ptr + src->step*y;
 965                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
 966                     for( x = 0; x < size.width; x++,
 967                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
 968                     {
 969                         //convert data
 970                         data[0]=float(sptr[0]),data[1]=float(sptr[1]),data[2]=float(sptr[2]);
 971                         //update GMM model
 972                         int result = _icvUpdateGMM_C3(data[0],data[1],data[2],pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
 973                         //detect shadows in the foreground
 974                         if (bShadowDetection)
 975                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
 976                         //generate output
 977                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
 978                     }
 979                 }
 980                 break;
 981             case CV_16S:
 982                 for( y = 0; y < size.height; y++ )
 983                 {
 984                     short* sptr = src->data.s + src->step*y;
 985                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
 986                     for( x = 0; x < size.width; x++,
 987                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
 988                     {
 989                         //convert data
 990                         data[0]=float(sptr[0]),data[1]=float(sptr[1]),data[2]=float(sptr[2]);
 991                         //update GMM model
 992                         int result = _icvUpdateGMM_C3(data[0],data[1],data[2],pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
 993                         //detect shadows in the foreground
 994                         if (bShadowDetection)
 995                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
 996                         //generate output
 997                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
 998                     }
 999                 }
1000                 break;
1001             case CV_16U:
1002                 for( y = 0; y < size.height; y++ )
1003                 {
1004                     unsigned short* sptr = (unsigned short*) src->data.s + src->step*y;
1005                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
1006                     for( x = 0; x < size.width; x++,
1007                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
1008                     {
1009                         //convert data
1010                         data[0]=float(sptr[0]),data[1]=float(sptr[1]),data[2]=float(sptr[2]);
1011                         //update GMM model
1012                         int result = _icvUpdateGMM_C3(data[0],data[1],data[2],pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
1013                         //detect shadows in the foreground
1014                         if (bShadowDetection)
1015                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
1016                         //generate output
1017                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
1018                     }
1019                 }
1020                 break;
1021             case CV_32S:
1022                 for( y = 0; y < size.height; y++ )
1023                 {
1024                     int* sptr = src->data.i + src->step*y;
1025                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
1026                     for( x = 0; x < size.width; x++,
1027                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
1028                     {
1029                         //convert data
1030                         data[0]=float(sptr[0]),data[1]=float(sptr[1]),data[2]=float(sptr[2]);
1031                         //update GMM model
1032                         int result = _icvUpdateGMM_C3(data[0],data[1],data[2],pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
1033                         //detect shadows in the foreground
1034                         if (bShadowDetection)
1035                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
1036                         //generate output
1037                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
1038                     }
1039                 }
1040                 break;
1041             case CV_32F:
1042                 for( y = 0; y < size.height; y++ )
1043                 {
1044                     float* sptr = src->data.fl + src->step*y;
1045                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
1046                     for( x = 0; x < size.width; x++,
1047                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
1048                     {
1049                         //update GMM model
1050                         int result = _icvUpdateGMM_C3(sptr[0],sptr[1],sptr[2],pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
1051                         //detect shadows in the foreground
1052                         if (bShadowDetection)
1053                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
1054                         //generate output
1055                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
1056                     }
1057                 }
1058                 break;
1059             case CV_64F:
1060                 for( y = 0; y < size.height; y++ )
1061                 {
1062                     double* sptr = src->data.db + src->step*y;
1063                     uchar* pDataOutput = dst->data.ptr + dst->step*y;
1064                     for( x = 0; x < size.width; x++,
1065                         pGMM+=nM,pUsedModes++,pDataOutput++,sptr+=nD)
1066                     {
1067                         //convert data
1068                         data[0]=float(sptr[0]),data[1]=float(sptr[1]),data[2]=float(sptr[2]);
1069                         //update GMM model
1070                         int result = _icvUpdateGMM_C3(data[0],data[1],data[2],pUsedModes,pGMM,nM,alpha, fTb, fTB, fTg, fVarInit, fVarMax, fVarMin,prune);
1071                         //detect shadows in the foreground
1072                         if (bShadowDetection)
1073                             if (result==0) result= _icvRemoveShadowGMM(data,nD,(*pUsedModes),pGMM,fTb,fTB,fTau);
1074                         //generate output
1075                         (* pDataOutput)= (result==1) ? 0 : (result==2) ? (nShadowDetection) : 255;
1076                     }
1077                 }
1078                 break;
1079         }
1080     }//a bit faster for nD=3;
1081 }
1082
1083
1084 //only foreground image is updated
1085 //no filtering included
1086 typedef struct CvGaussBGModel2
1087 {
1088     CV_BG_STAT_MODEL_FIELDS();
1089     CvGaussBGStatModel2Params params;
1090     CvGaussBGStatModel2Data   data;
1091     int                       countFrames;
1092 } CvGaussBGModel2;
1093
1094 CVAPI(CvBGStatModel*) cvCreateGaussianBGModel2( IplImage* first_frame,
1095                                                CvGaussBGStatModel2Params* params CV_DEFAULT(NULL) );
1096
1097 //////////////////////////////////////////////
1098 //implementation as part of the CvBGStatModel
1099 static void CV_CDECL icvReleaseGaussianBGModel2( CvGaussBGModel2** bg_model );
1100 static int CV_CDECL icvUpdateGaussianBGModel2( IplImage* curr_frame, CvGaussBGModel2*  bg_model );
1101
1102
1103 CV_IMPL CvBGStatModel*
1104 cvCreateGaussianBGModel2( IplImage* first_frame, CvGaussBGStatModel2Params* parameters )
1105 {
1106     CvGaussBGModel2* bg_model = 0;
1107     int w,h;
1108
1109     CV_FUNCNAME( "cvCreateGaussianBGModel2" );
1110
1111     __BEGIN__;
1112
1113     CvGaussBGStatModel2Params params;
1114
1115     if( !CV_IS_IMAGE(first_frame) )
1116         CV_ERROR( CV_StsBadArg, "Invalid or NULL first_frame parameter" );
1117
1118     if( first_frame->nChannels>CV_BGFG_MOG2_NDMAX )
1119         CV_ERROR( CV_StsBadArg, "Maxumum number of channels in the image is excedded (change CV_BGFG_MOG2_MAXBANDS constant)!" );
1120
1121
1122     CV_CALL( bg_model = (CvGaussBGModel2*)cvAlloc( sizeof(*bg_model) ));
1123     memset( bg_model, 0, sizeof(*bg_model) );
1124     bg_model->type    = CV_BG_MODEL_MOG2;
1125     bg_model->release = (CvReleaseBGStatModel) icvReleaseGaussianBGModel2;
1126     bg_model->update  = (CvUpdateBGStatModel)  icvUpdateGaussianBGModel2;
1127
1128     //init parameters
1129     if( parameters == NULL )
1130     {
1131         memset(&params, 0, sizeof(params));
1132
1133         // These constants are defined in cvaux/include/cvaux.h
1134         params.bShadowDetection = 1;
1135         params.bPostFiltering=0;
1136         params.minArea=CV_BGFG_MOG2_MINAREA;
1137
1138         //set parameters
1139         // K - max number of Gaussians per pixel
1140         params.nM = CV_BGFG_MOG2_NGAUSSIANS;//4;
1141         // Tb - the threshold - n var
1142         //pGMM->fTb = 4*4;
1143         params.fTb = CV_BGFG_MOG2_STD_THRESHOLD*CV_BGFG_MOG2_STD_THRESHOLD;
1144         // Tbf - the threshold
1145         //pGMM->fTB = 0.9f;//1-cf from the paper
1146         params.fTB = CV_BGFG_MOG2_BACKGROUND_THRESHOLD;
1147         // Tgenerate - the threshold
1148         params.fTg = CV_BGFG_MOG2_STD_THRESHOLD_GENERATE*CV_BGFG_MOG2_STD_THRESHOLD_GENERATE;//update the mode or generate new
1149         //pGMM->fSigma= 11.0f;//sigma for the new mode
1150         params.fVarInit = CV_BGFG_MOG2_VAR_INIT;
1151         params.fVarMax = CV_BGFG_MOG2_VAR_MAX;
1152         params.fVarMin = CV_BGFG_MOG2_VAR_MIN;
1153         // alpha - the learning factor
1154         params.fAlphaT = 1.0f/CV_BGFG_MOG2_WINDOW_SIZE;//0.003f;
1155         // complexity reduction prior constant
1156         params.fCT = CV_BGFG_MOG2_CT;//0.05f;
1157
1158         //shadow
1159         // Shadow detection
1160         params.nShadowDetection = (unsigned char)CV_BGFG_MOG2_SHADOW_VALUE;//value 0 to turn off
1161         params.fTau = CV_BGFG_MOG2_SHADOW_TAU;//0.5f;// Tau - shadow threshold
1162     }
1163     else
1164     {
1165         params = *parameters;
1166     }
1167
1168     bg_model->params = params;
1169
1170     //image data
1171     w = first_frame->width;
1172     h = first_frame->height;
1173
1174     bg_model->params.nWidth = w;
1175     bg_model->params.nHeight = h;
1176
1177     bg_model->params.nND = first_frame->nChannels;
1178
1179
1180     //allocate GMM data
1181
1182     //GMM for each pixel
1183     bg_model->data.rGMM = (CvPBGMMGaussian*) malloc(w*h * params.nM * sizeof(CvPBGMMGaussian));
1184     //used modes per pixel
1185     bg_model->data.rnUsedModes = (unsigned char* ) malloc(w*h);
1186     memset(bg_model->data.rnUsedModes,0,w*h);//no modes used
1187
1188     //prepare storages
1189     CV_CALL( bg_model->background = cvCreateImage(cvSize(w,h), IPL_DEPTH_8U, first_frame->nChannels));
1190     CV_CALL( bg_model->foreground = cvCreateImage(cvSize(w,h), IPL_DEPTH_8U, 1));
1191
1192     //for eventual filtering
1193     CV_CALL( bg_model->storage = cvCreateMemStorage());
1194
1195     bg_model->countFrames = 0;
1196
1197     __END__;
1198
1199     if( cvGetErrStatus() < 0 )
1200     {
1201         CvBGStatModel* base_ptr = (CvBGStatModel*)bg_model;
1202
1203         if( bg_model && bg_model->release )
1204             bg_model->release( &base_ptr );
1205         else
1206             cvFree( &bg_model );
1207         bg_model = 0;
1208     }
1209
1210     return (CvBGStatModel*)bg_model;
1211 }
1212
1213
1214 static void CV_CDECL
1215 icvReleaseGaussianBGModel2( CvGaussBGModel2** _bg_model )
1216 {
1217     CV_FUNCNAME( "icvReleaseGaussianBGModel2" );
1218
1219     __BEGIN__;
1220
1221     if( !_bg_model )
1222         CV_ERROR( CV_StsNullPtr, "" );
1223
1224     if( *_bg_model )
1225     {
1226         CvGaussBGModel2* bg_model = *_bg_model;
1227
1228         free (bg_model->data.rGMM);
1229         free (bg_model->data.rnUsedModes);
1230
1231         cvReleaseImage( &bg_model->background );
1232         cvReleaseImage( &bg_model->foreground );
1233         cvReleaseMemStorage(&bg_model->storage);
1234         memset( bg_model, 0, sizeof(*bg_model) );
1235         cvFree( _bg_model );
1236     }
1237
1238     __END__;
1239 }
1240
1241
1242 static int CV_CDECL
1243 icvUpdateGaussianBGModel2( IplImage* curr_frame, CvGaussBGModel2*  bg_model )
1244 {
1245     //checks
1246     if ((curr_frame->height!=bg_model->params.nHeight)||(curr_frame->width!=bg_model->params.nWidth)||(curr_frame->nChannels!=bg_model->params.nND))
1247         CV_Error( CV_StsBadSize, "the image not the same size as the reserved GMM background model");
1248
1249     float alpha=bg_model->params.fAlphaT;
1250     bg_model->countFrames++;
1251
1252     //faster initial updates - increase value of alpha
1253     if (bg_model->params.bInit){
1254         float alphaInit=(1.0f/(2*bg_model->countFrames+1));
1255         if (alphaInit>alpha)
1256         {
1257             alpha = alphaInit;
1258         }
1259         else
1260         {
1261             bg_model->params.bInit = 0;
1262         }
1263     }
1264
1265     //update background
1266     //icvUpdatePixelBackgroundGMM2( curr_frame, bg_model->foreground, bg_model->data.rGMM,bg_model->data.rnUsedModes,&(bg_model->params),alpha);
1267     icvUpdatePixelBackgroundGMM2( curr_frame, bg_model->foreground, bg_model->data.rGMM,bg_model->data.rnUsedModes,
1268                                  bg_model->params.nM,
1269                                  bg_model->params.fTb,
1270                                  bg_model->params.fTB,
1271                                  bg_model->params.fTg,
1272                                  bg_model->params.fVarInit,
1273                                  bg_model->params.fVarMax,
1274                                  bg_model->params.fVarMin,
1275                                  bg_model->params.fCT,
1276                                  bg_model->params.fTau,
1277                                  bg_model->params.bShadowDetection,
1278                                  bg_model->params.nShadowDetection,
1279                                  alpha);
1280
1281     //foreground filtering
1282     if (bg_model->params.bPostFiltering==1)
1283     {
1284         int region_count = 0;
1285         CvSeq *first_seq = NULL, *prev_seq = NULL, *seq = NULL;
1286
1287
1288         //filter small regions
1289         cvClearMemStorage(bg_model->storage);
1290
1291         cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_OPEN, 1 );
1292         cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_CLOSE, 1 );
1293
1294         cvFindContours( bg_model->foreground, bg_model->storage, &first_seq, sizeof(CvContour), CV_RETR_LIST );
1295         for( seq = first_seq; seq; seq = seq->h_next )
1296         {
1297             CvContour* cnt = (CvContour*)seq;
1298             if( cnt->rect.width * cnt->rect.height < bg_model->params.minArea )
1299             {
1300                 //delete small contour
1301                 prev_seq = seq->h_prev;
1302                 if( prev_seq )
1303                 {
1304                     prev_seq->h_next = seq->h_next;
1305                     if( seq->h_next ) seq->h_next->h_prev = prev_seq;
1306                 }
1307                 else
1308                 {
1309                     first_seq = seq->h_next;
1310                     if( seq->h_next ) seq->h_next->h_prev = NULL;
1311                 }
1312             }
1313             else
1314             {
1315                 region_count++;
1316             }
1317         }
1318         bg_model->foreground_regions = first_seq;
1319         cvZero(bg_model->foreground);
1320         cvDrawContours(bg_model->foreground, first_seq, CV_RGB(0, 0, 255), CV_RGB(0, 0, 255), 10, -1);
1321
1322         return region_count;
1323     }
1324
1325     return 1;
1326 }
1327
1328 /* End of file. */