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
  71 BackgroundSubtractorMOG::BackgroundSubtractorMOG()
  72 {
  73     frameSize = Size(0,0);
  74     frameType = 0;
  75
  76     nframes = 0;
  77     nmixtures = defaultNMixtures;
  78     history = defaultHistory;
  79     varThreshold = defaultVarThreshold;
  80     backgroundRatio = defaultBackgroundRatio;
  81 }
  82
  83 BackgroundSubtractorMOG::BackgroundSubtractorMOG(int _history, int _nmixtures, double _backgroundRatio)
  84 {
  85     frameSize = Size(0,0);
  86     frameType = 0;
  87
  88     nframes = 0;
  89     nmixtures = min(_nmixtures > 0 ? _nmixtures : defaultNMixtures, 8);
  90     history = _history > 0 ? _history : defaultHistory;
  91     varThreshold = defaultVarThreshold;
  92     backgroundRatio = min(_backgroundRatio > 0 ? _backgroundRatio : 0.95, 1.);
  93 }
  94
  95 BackgroundSubtractorMOG::~BackgroundSubtractorMOG()
  96 {
  97 }
  98
  99
 100 void BackgroundSubtractorMOG::initialize(Size _frameSize, int _frameType)
 101 {
 102     frameSize = _frameSize;
 103     frameType = _frameType;
 104     nframes = 0;
 105
 106     int nchannels = CV_MAT_CN(frameType);
 107     CV_Assert( CV_MAT_DEPTH(frameType) == CV_8U );
 108
 109     // for each gaussian mixture of each pixel bg model we store ...
 110     // the mixture sort key (w/sum_of_variances), the mixture weight (w),
 111     // the mean (nchannels values) and
 112     // the diagonal covariance matrix (another nchannels values)
 113     bgmodel.create( frameSize.height, frameSize.width*nmixtures*(2 + 2*nchannels), CV_32F );
 114     const float w0 = (float)CV_BGFG_MOG_WEIGHT_INIT;
 115     const float var0 = (float)(CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT);
 116     const float sk0 = (float)(CV_BGFG_MOG_WEIGHT_INIT/(CV_BGFG_MOG_SIGMA_INIT*sqrt((double)nchannels)));
 117
 118     for( int y = 0; y < frameSize.height; y++ )
 119     {
 120         float* mptr = bgmodel.ptr<float>(y);
 121         for( int x = 0; x < frameSize.width; x++ )
 122         {
 123             for( int k = 0; k < nmixtures; k++ )
 124             {
 125                 mptr[0] = sk0;
 126                 mptr[1] = w0;
 127                 mptr += 2;
 128                 for( int c = 0; c < nchannels; c++ )
 129                 {
 130                     mptr[c] = 0;
 131                     mptr[c + nchannels] = var0;
 132                 }
 133                 mptr += nchannels*2;
 134             }
 135         }
 136     }
 137 }
 138
 139
 140 template<typename VT> struct MixData
 141 {
 142     float sortKey;
 143     float weight;
 144     VT mean;
 145     VT var;
 146 };
 147
 148
 149 static void process8uC1( BackgroundSubtractorMOG& obj, const Mat& image, Mat& fgmask, double learningRate )
 150 {
 151     int x, y, k, k1, rows = image.rows, cols = image.cols;
 152     float alpha = (float)learningRate, T = (float)obj.backgroundRatio, vT = (float)obj.varThreshold;
 153     int K = obj.nmixtures;
 154
 155     const float w0 = (float)CV_BGFG_MOG_WEIGHT_INIT;
 156     const float sk0 = (float)(CV_BGFG_MOG_WEIGHT_INIT/CV_BGFG_MOG_SIGMA_INIT);
 157     const float var0 = (float)(CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT);
 158
 159     for( y = 0; y < rows; y++ )
 160     {
 161         const uchar* src = image.ptr<uchar>(y);
 162         uchar* dst = fgmask.ptr<uchar>(y);
 163         MixData<float>* mptr = (MixData<float>*)obj.bgmodel.ptr(y);
 164
 165         for( x = 0; x < cols; x++, mptr += K )
 166         {
 167             float wsum = 0, dw = 0;
 168             float pix = src[x];
 169             for( k = 0; k < K; k++ )
 170             {
 171                 float w = mptr[k].weight;
 172                 float mu = mptr[k].mean;
 173                 float var = mptr[k].var;
 174                 float diff = pix - mu, d2 = diff*diff;
 175                 if( d2 < vT*var )
 176                 {
 177                     dw = alpha*(1.f - w);
 178                     mptr[k].weight = w + dw;
 179                     mptr[k].mean = mu + alpha*diff;
 180                     mptr[k].var = var = max(var + alpha*(d2 - var), FLT_EPSILON);
 181                     mptr[k].sortKey = w/sqrt(var);
 182
 183                     for( k1 = k-1; k1 >= 0; k1-- )
 184                     {
 185                         if( mptr[k1].sortKey > mptr[k1+1].sortKey )
 186                             break;
 187                         std::swap( mptr[k1], mptr[k1+1] );
 188                     }
 189                     break;
 190                 }
 191                 wsum += w;
 192             }
 193
 194             dst[x] = (uchar)(-(wsum >= T));
 195             wsum += dw;
 196
 197             if( k == K ) // no appropriate gaussian mixture found at all, remove the weakest mixture and create a new one
 198             {
 199                 wsum += w0 - mptr[K-1].weight;
 200                 mptr[K-1].weight = w0;
 201                 mptr[K-1].mean = pix;
 202                 mptr[K-1].var = var0;
 203                 mptr[K-1].sortKey = sk0;
 204             }
 205             else
 206                 for( ; k < K; k++ )
 207                     wsum += mptr[k].weight;
 208
 209             dw = 1.f/wsum;
 210             for( k = 0; k < K; k++ )
 211             {
 212                 mptr[k].weight *= dw;
 213                 mptr[k].sortKey *= dw;
 214             }
 215         }
 216     }
 217 }
 218
 219 static void process8uC3( BackgroundSubtractorMOG& obj, const Mat& image, Mat& fgmask, double learningRate )
 220 {
 221     int x, y, k, k1, rows = image.rows, cols = image.cols;
 222     float alpha = (float)learningRate, T = (float)obj.backgroundRatio, vT = (float)obj.varThreshold;
 223     int K = obj.nmixtures;
 224
 225     const float w0 = (float)CV_BGFG_MOG_WEIGHT_INIT;
 226     const float sk0 = (float)(CV_BGFG_MOG_WEIGHT_INIT/CV_BGFG_MOG_SIGMA_INIT);
 227     const float var0 = (float)(CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT);
 228
 229     for( y = 0; y < rows; y++ )
 230     {
 231         const uchar* src = image.ptr<uchar>(y);
 232         uchar* dst = fgmask.ptr<uchar>(y);
 233         MixData<Vec3f>* mptr = (MixData<Vec3f>*)obj.bgmodel.ptr(y);
 234
 235         for( x = 0; x < cols; x++, mptr += K )
 236         {
 237             float wsum = 0, dw = 0;
 238             Vec3f pix(src[x*3], src[x*3+1], src[x*3+2]);
 239             for( k = 0; k < K; k++ )
 240             {
 241                 float w = mptr[k].weight;
 242                 Vec3f mu = mptr[k].mean[0];
 243                 Vec3f var = mptr[k].var[0];
 244                 Vec3f diff = pix - mu;
 245                 float d2 = diff.dot(diff);
 246                 if( d2 < vT*(var[0] + var[1] + var[2]) )
 247                 {
 248                     dw = alpha*(1.f - w);
 249                     mptr[k].weight = w + dw;
 250                     mptr[k].mean = mu + alpha*diff;
 251                     var = Vec3f(max(var[0] + alpha*(diff[0]*diff[0] - var[0]), FLT_EPSILON),
 252                                 max(var[1] + alpha*(diff[1]*diff[1] - var[1]), FLT_EPSILON),
 253                                 max(var[2] + alpha*(diff[2]*diff[2] - var[2]), FLT_EPSILON));
 254                     mptr[k].var = var;
 255                     mptr[k].sortKey = w/sqrt(var[0] + var[1] + var[2]);
 256
 257                     for( k1 = k-1; k1 >= 0; k1-- )
 258                     {
 259                         if( mptr[k1].sortKey > mptr[k1+1].sortKey )
 260                             break;
 261                         std::swap( mptr[k1], mptr[k1+1] );
 262                     }
 263                     break;
 264                 }
 265                 wsum += w;
 266             }
 267
 268             dst[x] = (uchar)(-(wsum >= T));
 269             wsum += dw;
 270
 271             if( k == K ) // no appropriate gaussian mixture found at all, remove the weakest mixture and create a new one
 272             {
 273                 wsum += w0 - mptr[K-1].weight;
 274                 mptr[K-1].weight = w0;
 275                 mptr[K-1].mean = pix;
 276                 mptr[K-1].var = Vec3f(var0, var0, var0);
 277                 mptr[K-1].sortKey = sk0;
 278             }
 279             else
 280                 for( ; k < K; k++ )
 281                     wsum += mptr[k].weight;
 282
 283             dw = 1.f/wsum;
 284             for( k = 0; k < K; k++ )
 285             {
 286                 mptr[k].weight *= dw;
 287                 mptr[k].sortKey *= dw;
 288             }
 289         }
 290     }
 291 }
 292
 293 void BackgroundSubtractorMOG::operator()(const Mat& image, Mat& fgmask, double learningRate)
 294 {
 295     bool needToInitialize = nframes == 0 || learningRate >= 1 || image.size() != frameSize || image.type() != frameType;
 296
 297     if( needToInitialize )
 298         initialize(frameSize, frameType);
 299
 300     CV_Assert( image.depth() == CV_8U );
 301     fgmask.create( image.size(), CV_8U );
 302
 303     ++nframes;
 304     learningRate = learningRate >= 0 && nframes > 1 ? learningRate : 1./min( nframes, history );
 305
 306     if( image.type() == CV_8UC1 )
 307         process8uC1( *this, image, fgmask, learningRate );
 308     else if( image.type() == CV_8UC3 )
 309         process8uC3( *this, image, fgmask, learningRate );
 310     else
 311         CV_Error( CV_StsUnsupportedFormat, "Only 1- and 3-channel 8-bit images are supported in BackgroundSubtractorMOG" );
 312 }
 313
 314 }
 315
 316
 317 static void CV_CDECL
 318 icvReleaseGaussianBGModel( CvGaussBGModel** bg_model )
 319 {
 320     if( !bg_model )
 321         CV_Error( CV_StsNullPtr, "" );
 322
 323     if( *bg_model )
 324     {
 325         delete (cv::Mat*)((*bg_model)->g_point);
 326         cvReleaseImage( &(*bg_model)->background );
 327         cvReleaseImage( &(*bg_model)->foreground );
 328         cvReleaseMemStorage(&(*bg_model)->storage);
 329         memset( *bg_model, 0, sizeof(**bg_model) );
 330         delete *bg_model;
 331         *bg_model = 0;
 332     }
 333 }
 334
 335
 336 static int CV_CDECL
 337 icvUpdateGaussianBGModel( IplImage* curr_frame, CvGaussBGModel*  bg_model, double learningRate )
 338 {
 339     int region_count = 0;
 340     CvSeq *first_seq = NULL, *prev_seq = NULL, *seq = NULL;
 341
 342     cv::Mat image = cv::cvarrToMat(curr_frame), mask = cv::cvarrToMat(bg_model->foreground);
 343
 344     cv::BackgroundSubtractorMOG mog;
 345     mog.bgmodel = *(cv::Mat*)bg_model->g_point;
 346     mog.frameSize = mog.bgmodel.data ? cv::Size(cvGetSize(curr_frame)) : cv::Size();
 347     mog.frameType = image.type();
 348
 349     mog.nframes = bg_model->countFrames;
 350     mog.history = bg_model->params.win_size;
 351     mog.nmixtures = bg_model->params.n_gauss;
 352     mog.varThreshold = bg_model->params.std_threshold;
 353     mog.backgroundRatio = bg_model->params.bg_threshold;
 354
 355     mog(image, mask, learningRate);
 356
 357     bg_model->countFrames = mog.nframes;
 358     if( ((cv::Mat*)bg_model->g_point)->data != mog.bgmodel.data )
 359         *((cv::Mat*)bg_model->g_point) = mog.bgmodel;
 360
 361     //foreground filtering
 362
 363     //filter small regions
 364     cvClearMemStorage(bg_model->storage);
 365
 366     //cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_OPEN, 1 );
 367     //cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_CLOSE, 1 );
 368
 369     cvFindContours( bg_model->foreground, bg_model->storage, &first_seq, sizeof(CvContour), CV_RETR_LIST );
 370     for( seq = first_seq; seq; seq = seq->h_next )
 371     {
 372         CvContour* cnt = (CvContour*)seq;
 373         if( cnt->rect.width * cnt->rect.height < bg_model->params.minArea )
 374         {
 375             //delete small contour
 376             prev_seq = seq->h_prev;
 377             if( prev_seq )
 378             {
 379                 prev_seq->h_next = seq->h_next;
 380                 if( seq->h_next ) seq->h_next->h_prev = prev_seq;
 381             }
 382             else
 383             {
 384                 first_seq = seq->h_next;
 385                 if( seq->h_next ) seq->h_next->h_prev = NULL;
 386             }
 387         }
 388         else
 389         {
 390             region_count++;
 391         }
 392     }
 393     bg_model->foreground_regions = first_seq;
 394     cvZero(bg_model->foreground);
 395     cvDrawContours(bg_model->foreground, first_seq, CV_RGB(0, 0, 255), CV_RGB(0, 0, 255), 10, -1);
 396
 397     return region_count;
 398 }
 399
 400 CV_IMPL CvBGStatModel*
 401 cvCreateGaussianBGModel( IplImage* first_frame, CvGaussBGStatModelParams* parameters )
 402 {
 403     CvGaussBGStatModelParams params;
 404
 405     CV_Assert( CV_IS_IMAGE(first_frame) );
 406
 407     //init parameters
 408     if( parameters == NULL )
 409     {                        /* These constants are defined in cvaux/include/cvaux.h: */
 410         params.win_size      = CV_BGFG_MOG_WINDOW_SIZE;
 411         params.bg_threshold  = CV_BGFG_MOG_BACKGROUND_THRESHOLD;
 412
 413         params.std_threshold = CV_BGFG_MOG_STD_THRESHOLD;
 414         params.weight_init   = CV_BGFG_MOG_WEIGHT_INIT;
 415
 416         params.variance_init = CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT;
 417         params.minArea       = CV_BGFG_MOG_MINAREA;
 418         params.n_gauss       = CV_BGFG_MOG_NGAUSSIANS;
 419     }
 420     else
 421         params = *parameters;
 422
 423     CvGaussBGModel* bg_model = new CvGaussBGModel;
 424     memset( bg_model, 0, sizeof(*bg_model) );
 425     bg_model->type = CV_BG_MODEL_MOG;
 426     bg_model->release = (CvReleaseBGStatModel)icvReleaseGaussianBGModel;
 427     bg_model->update = (CvUpdateBGStatModel)icvUpdateGaussianBGModel;
 428
 429     bg_model->params = params;
 430
 431     //prepare storages
 432     bg_model->g_point = (CvGaussBGPoint*)new cv::Mat();
 433
 434     bg_model->background = cvCreateImage(cvSize(first_frame->width,
 435         first_frame->height), IPL_DEPTH_8U, first_frame->nChannels);
 436     bg_model->foreground = cvCreateImage(cvSize(first_frame->width,
 437         first_frame->height), IPL_DEPTH_8U, 1);
 438
 439     bg_model->storage = cvCreateMemStorage();
 440
 441     bg_model->countFrames = 0;
 442
 443     icvUpdateGaussianBGModel( first_frame, bg_model, 1 );
 444
 445     return (CvBGStatModel*)bg_model;
 446 }
 447
 448 /* End of file. */