opencv/src/cvaux/cvbgfg_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 "_cvaux.h"
  42 #include <float.h>
  43
  44 // to make sure we can use these short names
  45 #undef K
  46 #undef L
  47 #undef T
  48
  49 // This is based on the "An Improved Adaptive Background Mixture Model for
  50 // Real-time Tracking with Shadow Detection" by P. KaewTraKulPong and R. Bowden
  51 // http://personal.ee.surrey.ac.uk/Personal/R.Bowden/publications/avbs01/avbs01.pdf
  52 //
  53 // The windowing method is used, but not the shadow detection. I make some of my
  54 // own modifications which make more sense. There are some errors in some of their
  55 // equations.
  56 //
  57
  58 namespace cv
  59 {
  60
  61 BackgroundSubtractor::~BackgroundSubtractor() {}
  62 void BackgroundSubtractor::operator()(const Mat&, Mat&, double)
  63 {
  64 }
  65
  66 static const int defaultNMixtures = CV_BGFG_MOG_NGAUSSIANS;
  67 static const int defaultHistory = CV_BGFG_MOG_WINDOW_SIZE;
  68 static const double defaultBackgroundRatio = CV_BGFG_MOG_BACKGROUND_THRESHOLD;
  69 static const double defaultVarThreshold = CV_BGFG_MOG_STD_THRESHOLD*CV_BGFG_MOG_STD_THRESHOLD;
  70 static const double defaultNoiseSigma = CV_BGFG_MOG_SIGMA_INIT*0.5;
  71
  72 BackgroundSubtractorMOG::BackgroundSubtractorMOG()
  73 {
  74     frameSize = Size(0,0);
  75     frameType = 0;
  76
  77     nframes = 0;
  78     nmixtures = defaultNMixtures;
  79     history = defaultHistory;
  80     varThreshold = defaultVarThreshold;
  81     backgroundRatio = defaultBackgroundRatio;
  82     noiseSigma = defaultNoiseSigma;
  83 }
  84
  85 BackgroundSubtractorMOG::BackgroundSubtractorMOG(int _history, int _nmixtures,
  86                                                  double _backgroundRatio,
  87                                                  double _noiseSigma)
  88 {
  89     frameSize = Size(0,0);
  90     frameType = 0;
  91
  92     nframes = 0;
  93     nmixtures = min(_nmixtures > 0 ? _nmixtures : defaultNMixtures, 8);
  94     history = _history > 0 ? _history : defaultHistory;
  95     varThreshold = defaultVarThreshold;
  96     backgroundRatio = min(_backgroundRatio > 0 ? _backgroundRatio : 0.95, 1.);
  97     noiseSigma = _noiseSigma <= 0 ? defaultNoiseSigma : _noiseSigma;
  98 }
  99
 100 BackgroundSubtractorMOG::~BackgroundSubtractorMOG()
 101 {
 102 }
 103
 104
 105 void BackgroundSubtractorMOG::initialize(Size _frameSize, int _frameType)
 106 {
 107     frameSize = _frameSize;
 108     frameType = _frameType;
 109     nframes = 0;
 110
 111     int nchannels = CV_MAT_CN(frameType);
 112     CV_Assert( CV_MAT_DEPTH(frameType) == CV_8U );
 113
 114     // for each gaussian mixture of each pixel bg model we store ...
 115     // the mixture sort key (w/sum_of_variances), the mixture weight (w),
 116     // the mean (nchannels values) and
 117     // the diagonal covariance matrix (another nchannels values)
 118     bgmodel.create( 1, frameSize.height*frameSize.width*nmixtures*(2 + 2*nchannels), CV_32F );
 119     bgmodel = Scalar::all(0);
 120 }
 121
 122
 123 template<typename VT> struct MixData
 124 {
 125     float sortKey;
 126     float weight;
 127     VT mean;
 128     VT var;
 129 };
 130
 131
 132 static void process8uC1( BackgroundSubtractorMOG& obj, const Mat& image, Mat& fgmask, double learningRate )
 133 {
 134     int x, y, k, k1, rows = image.rows, cols = image.cols;
 135     float alpha = (float)learningRate, T = (float)obj.backgroundRatio, vT = (float)obj.varThreshold;
 136     int K = obj.nmixtures;
 137     MixData<float>* mptr = (MixData<float>*)obj.bgmodel.data;
 138
 139     const float w0 = (float)CV_BGFG_MOG_WEIGHT_INIT;
 140     const float sk0 = (float)(w0/CV_BGFG_MOG_SIGMA_INIT);
 141     const float var0 = (float)(CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT);
 142     const float minVar = (float)(obj.noiseSigma*obj.noiseSigma);
 143
 144     for( y = 0; y < rows; y++ )
 145     {
 146         const uchar* src = image.ptr<uchar>(y);
 147         uchar* dst = fgmask.ptr<uchar>(y);
 148
 149         if( alpha > 0 )
 150         {
 151             for( x = 0; x < cols; x++, mptr += K )
 152             {
 153                 float wsum = 0;
 154                 float pix = src[x];
 155                 int kHit = -1, kForeground = -1;
 156
 157                 for( k = 0; k < K; k++ )
 158                 {
 159                     float w = mptr[k].weight;
 160                     wsum += w;
 161                     if( w < FLT_EPSILON )
 162                         break;
 163                     float mu = mptr[k].mean;
 164                     float var = mptr[k].var;
 165                     float diff = pix - mu;
 166                     float d2 = diff*diff;
 167                     if( d2 < vT*var )
 168                     {
 169                         wsum -= w;
 170                         float dw = alpha*(1.f - w);
 171                         mptr[k].weight = w + dw;
 172                         mptr[k].mean = mu + alpha*diff;
 173                         var = max(var + alpha*(d2 - var), minVar);
 174                         mptr[k].var = var;
 175                         mptr[k].sortKey = w/sqrt(var);
 176
 177                         for( k1 = k-1; k1 >= 0; k1-- )
 178                         {
 179                             if( mptr[k1].sortKey >= mptr[k1+1].sortKey )
 180                                 break;
 181                             std::swap( mptr[k1], mptr[k1+1] );
 182                         }
 183
 184                         kHit = k1+1;
 185                         break;
 186                     }
 187                 }
 188
 189                 if( kHit < 0 ) // no appropriate gaussian mixture found at all, remove the weakest mixture and create a new one
 190                 {
 191                     kHit = k = min(k, K-1);
 192                     wsum += w0 - mptr[k].weight;
 193                     mptr[k].weight = w0;
 194                     mptr[k].mean = pix;
 195                     mptr[k].var = var0;
 196                     mptr[k].sortKey = sk0;
 197                 }
 198                 else
 199                     for( ; k < K; k++ )
 200                         wsum += mptr[k].weight;
 201
 202                 float wscale = 1.f/wsum;
 203                 wsum = 0;
 204                 for( k = 0; k < K; k++ )
 205                 {
 206                     wsum += mptr[k].weight *= wscale;
 207                     mptr[k].sortKey *= wscale;
 208                     if( wsum > T && kForeground < 0 )
 209                         kForeground = k+1;
 210                 }
 211
 212                 dst[x] = (uchar)(-(kHit >= kForeground));
 213             }
 214         }
 215         else
 216         {
 217             for( x = 0; x < cols; x++, mptr += K )
 218             {
 219                 float pix = src[x];
 220                 int kHit = -1, kForeground = -1;
 221
 222                 for( k = 0; k < K; k++ )
 223                 {
 224                     if( mptr[k].weight < FLT_EPSILON )
 225                         break;
 226                     float mu = mptr[k].mean;
 227                     float var = mptr[k].var;
 228                     float diff = pix - mu;
 229                     float d2 = diff*diff;
 230                     if( d2 < vT*var )
 231                     {
 232                         kHit = k;
 233                         break;
 234                     }
 235                 }
 236
 237                 if( kHit >= 0 )
 238                 {
 239                     float wsum = 0;
 240                     for( k = 0; k < K; k++ )
 241                     {
 242                         wsum += mptr[k].weight;
 243                         if( wsum > T )
 244                         {
 245                             kForeground = k+1;
 246                             break;
 247                         }
 248                     }
 249                 }
 250
 251                 dst[x] = (uchar)(kHit < 0 || kHit >= kForeground ? 255 : 0);
 252             }
 253         }
 254     }
 255 }
 256
 257 static void process8uC3( BackgroundSubtractorMOG& obj, const Mat& image, Mat& fgmask, double learningRate )
 258 {
 259     int x, y, k, k1, rows = image.rows, cols = image.cols;
 260     float alpha = (float)learningRate, T = (float)obj.backgroundRatio, vT = (float)obj.varThreshold;
 261     int K = obj.nmixtures;
 262
 263     const float w0 = (float)CV_BGFG_MOG_WEIGHT_INIT;
 264     const float sk0 = (float)(w0/CV_BGFG_MOG_SIGMA_INIT*sqrt(3.));
 265     const float var0 = (float)(CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT);
 266     const float minVar = (float)(obj.noiseSigma*obj.noiseSigma);
 267     MixData<Vec3f>* mptr = (MixData<Vec3f>*)obj.bgmodel.data;
 268
 269     for( y = 0; y < rows; y++ )
 270     {
 271         const uchar* src = image.ptr<uchar>(y);
 272         uchar* dst = fgmask.ptr<uchar>(y);
 273
 274         if( alpha > 0 )
 275         {
 276             for( x = 0; x < cols; x++, mptr += K )
 277             {
 278                 float wsum = 0;
 279                 Vec3f pix(src[x*3], src[x*3+1], src[x*3+2]);
 280                 int kHit = -1, kForeground = -1;
 281
 282                 for( k = 0; k < K; k++ )
 283                 {
 284                     float w = mptr[k].weight;
 285                     wsum += w;
 286                     if( w < FLT_EPSILON )
 287                         break;
 288                     Vec3f mu = mptr[k].mean;
 289                     Vec3f var = mptr[k].var;
 290                     Vec3f diff = pix - mu;
 291                     float d2 = diff.dot(diff);
 292                     if( d2 < vT*(var[0] + var[1] + var[2]) )
 293                     {
 294                         wsum -= w;
 295                         float dw = alpha*(1.f - w);
 296                         mptr[k].weight = w + dw;
 297                         mptr[k].mean = mu + alpha*diff;
 298                         var = Vec3f(max(var[0] + alpha*(diff[0]*diff[0] - var[0]), minVar),
 299                                     max(var[1] + alpha*(diff[1]*diff[1] - var[1]), minVar),
 300                                     max(var[2] + alpha*(diff[2]*diff[2] - var[2]), minVar));
 301                         mptr[k].var = var;
 302                         mptr[k].sortKey = w/sqrt(var[0] + var[1] + var[2]);
 303
 304                         for( k1 = k-1; k1 >= 0; k1-- )
 305                         {
 306                             if( mptr[k1].sortKey >= mptr[k1+1].sortKey )
 307                                 break;
 308                             std::swap( mptr[k1], mptr[k1+1] );
 309                         }
 310
 311                         kHit = k1+1;
 312                         break;
 313                     }
 314                 }
 315
 316                 if( kHit < 0 ) // no appropriate gaussian mixture found at all, remove the weakest mixture and create a new one
 317                 {
 318                     kHit = k = min(k, K-1);
 319                     wsum += w0 - mptr[k].weight;
 320                     mptr[k].weight = w0;
 321                     mptr[k].mean = pix;
 322                     mptr[k].var = Vec3f(var0, var0, var0);
 323                     mptr[k].sortKey = sk0;
 324                 }
 325                 else
 326                     for( ; k < K; k++ )
 327                         wsum += mptr[k].weight;
 328
 329                 float wscale = 1.f/wsum;
 330                 wsum = 0;
 331                 for( k = 0; k < K; k++ )
 332                 {
 333                     wsum += mptr[k].weight *= wscale;
 334                     mptr[k].sortKey *= wscale;
 335                     if( wsum > T && kForeground < 0 )
 336                         kForeground = k+1;
 337                 }
 338
 339                 dst[x] = (uchar)(-(kHit >= kForeground));
 340             }
 341         }
 342         else
 343         {
 344             for( x = 0; x < cols; x++, mptr += K )
 345             {
 346                 Vec3f pix(src[x*3], src[x*3+1], src[x*3+2]);
 347                 int kHit = -1, kForeground = -1;
 348
 349                 for( k = 0; k < K; k++ )
 350                 {
 351                     if( mptr[k].weight < FLT_EPSILON )
 352                         break;
 353                     Vec3f mu = mptr[k].mean;
 354                     Vec3f var = mptr[k].var;
 355                     Vec3f diff = pix - mu;
 356                     float d2 = diff.dot(diff);
 357                     if( d2 < vT*(var[0] + var[1] + var[2]) )
 358                     {
 359                         kHit = k;
 360                         break;
 361                     }
 362                 }
 363
 364                 if( kHit >= 0 )
 365                 {
 366                     float wsum = 0;
 367                     for( k = 0; k < K; k++ )
 368                     {
 369                         wsum += mptr[k].weight;
 370                         if( wsum > T )
 371                         {
 372                             kForeground = k+1;
 373                             break;
 374                         }
 375                     }
 376                 }
 377
 378                 dst[x] = (uchar)(kHit < 0 || kHit >= kForeground ? 255 : 0);
 379             }
 380         }
 381     }
 382 }
 383
 384 void BackgroundSubtractorMOG::operator()(const Mat& image, Mat& fgmask, double learningRate)
 385 {
 386     bool needToInitialize = nframes == 0 || learningRate >= 1 || image.size() != frameSize || image.type() != frameType;
 387
 388     if( needToInitialize )
 389         initialize(image.size(), image.type());
 390
 391     CV_Assert( image.depth() == CV_8U );
 392     fgmask.create( image.size(), CV_8U );
 393
 394     ++nframes;
 395     learningRate = learningRate >= 0 && nframes > 1 ? learningRate : 1./min( nframes, history );
 396     CV_Assert(learningRate >= 0);
 397
 398     if( image.type() == CV_8UC1 )
 399         process8uC1( *this, image, fgmask, learningRate );
 400     else if( image.type() == CV_8UC3 )
 401         process8uC3( *this, image, fgmask, learningRate );
 402     else
 403         CV_Error( CV_StsUnsupportedFormat, "Only 1- and 3-channel 8-bit images are supported in BackgroundSubtractorMOG" );
 404 }
 405
 406 }
 407
 408
 409 static void CV_CDECL
 410 icvReleaseGaussianBGModel( CvGaussBGModel** bg_model )
 411 {
 412     if( !bg_model )
 413         CV_Error( CV_StsNullPtr, "" );
 414
 415     if( *bg_model )
 416     {
 417         delete (cv::Mat*)((*bg_model)->g_point);
 418         cvReleaseImage( &(*bg_model)->background );
 419         cvReleaseImage( &(*bg_model)->foreground );
 420         cvReleaseMemStorage(&(*bg_model)->storage);
 421         memset( *bg_model, 0, sizeof(**bg_model) );
 422         delete *bg_model;
 423         *bg_model = 0;
 424     }
 425 }
 426
 427
 428 static int CV_CDECL
 429 icvUpdateGaussianBGModel( IplImage* curr_frame, CvGaussBGModel*  bg_model, double learningRate )
 430 {
 431     int region_count = 0;
 432
 433     cv::Mat image = cv::cvarrToMat(curr_frame), mask = cv::cvarrToMat(bg_model->foreground);
 434
 435     cv::BackgroundSubtractorMOG mog;
 436     mog.bgmodel = *(cv::Mat*)bg_model->g_point;
 437     mog.frameSize = mog.bgmodel.data ? cv::Size(cvGetSize(curr_frame)) : cv::Size();
 438     mog.frameType = image.type();
 439
 440     mog.nframes = bg_model->countFrames;
 441     mog.history = bg_model->params.win_size;
 442     mog.nmixtures = bg_model->params.n_gauss;
 443     mog.varThreshold = bg_model->params.std_threshold;
 444     mog.backgroundRatio = bg_model->params.bg_threshold;
 445
 446     mog(image, mask, learningRate);
 447
 448     bg_model->countFrames = mog.nframes;
 449     if( ((cv::Mat*)bg_model->g_point)->data != mog.bgmodel.data )
 450         *((cv::Mat*)bg_model->g_point) = mog.bgmodel;
 451
 452     //foreground filtering
 453
 454     //filter small regions
 455     cvClearMemStorage(bg_model->storage);
 456
 457     //cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_OPEN, 1 );
 458     //cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_CLOSE, 1 );
 459
 460     /*
 461     CvSeq *first_seq = NULL, *prev_seq = NULL, *seq = NULL;
 462     cvFindContours( bg_model->foreground, bg_model->storage, &first_seq, sizeof(CvContour), CV_RETR_LIST );
 463     for( seq = first_seq; seq; seq = seq->h_next )
 464     {
 465         CvContour* cnt = (CvContour*)seq;
 466         if( cnt->rect.width * cnt->rect.height < bg_model->params.minArea )
 467         {
 468             //delete small contour
 469             prev_seq = seq->h_prev;
 470             if( prev_seq )
 471             {
 472                 prev_seq->h_next = seq->h_next;
 473                 if( seq->h_next ) seq->h_next->h_prev = prev_seq;
 474             }
 475             else
 476             {
 477                 first_seq = seq->h_next;
 478                 if( seq->h_next ) seq->h_next->h_prev = NULL;
 479             }
 480         }
 481         else
 482         {
 483             region_count++;
 484         }
 485     }
 486     bg_model->foreground_regions = first_seq;
 487     cvZero(bg_model->foreground);
 488     cvDrawContours(bg_model->foreground, first_seq, CV_RGB(0, 0, 255), CV_RGB(0, 0, 255), 10, -1);*/
 489     CvMat _mask = mask;
 490     cvCopy(&_mask, bg_model->foreground);
 491
 492     return region_count;
 493 }
 494
 495 CV_IMPL CvBGStatModel*
 496 cvCreateGaussianBGModel( IplImage* first_frame, CvGaussBGStatModelParams* parameters )
 497 {
 498     CvGaussBGStatModelParams params;
 499
 500     CV_Assert( CV_IS_IMAGE(first_frame) );
 501
 502     //init parameters
 503     if( parameters == NULL )
 504     {                        /* These constants are defined in cvaux/include/cvaux.h: */
 505         params.win_size      = CV_BGFG_MOG_WINDOW_SIZE;
 506         params.bg_threshold  = CV_BGFG_MOG_BACKGROUND_THRESHOLD;
 507
 508         params.std_threshold = CV_BGFG_MOG_STD_THRESHOLD;
 509         params.weight_init   = CV_BGFG_MOG_WEIGHT_INIT;
 510
 511         params.variance_init = CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT;
 512         params.minArea       = CV_BGFG_MOG_MINAREA;
 513         params.n_gauss       = CV_BGFG_MOG_NGAUSSIANS;
 514     }
 515     else
 516         params = *parameters;
 517
 518     CvGaussBGModel* bg_model = new CvGaussBGModel;
 519     memset( bg_model, 0, sizeof(*bg_model) );
 520     bg_model->type = CV_BG_MODEL_MOG;
 521     bg_model->release = (CvReleaseBGStatModel)icvReleaseGaussianBGModel;
 522     bg_model->update = (CvUpdateBGStatModel)icvUpdateGaussianBGModel;
 523
 524     bg_model->params = params;
 525
 526     //prepare storages
 527     bg_model->g_point = (CvGaussBGPoint*)new cv::Mat();
 528
 529     bg_model->background = cvCreateImage(cvSize(first_frame->width,
 530         first_frame->height), IPL_DEPTH_8U, first_frame->nChannels);
 531     bg_model->foreground = cvCreateImage(cvSize(first_frame->width,
 532         first_frame->height), IPL_DEPTH_8U, 1);
 533
 534     bg_model->storage = cvCreateMemStorage();
 535
 536     bg_model->countFrames = 0;
 537
 538     icvUpdateGaussianBGModel( first_frame, bg_model, 1 );
 539
 540     return (CvBGStatModel*)bg_model;
 541 }
 542
 543 /* End of file. */
 544