git.maemo.org Git - opencv/blob - cvaux/src/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
  42 // This is based on the "An Improved Adaptive Background Mixture Model for
  43 // Real-time Tracking with Shadow Detection" by P. KaewTraKulPong and R. Bowden
  44 // http://personal.ee.surrey.ac.uk/Personal/R.Bowden/publications/avbs01/avbs01.pdf
  45 //
  46 // The windowing method is used, but not the shadow detection. I make some of my
  47 // own modifications which make more sense. There are some errors in some of their
  48 // equations.
  49 //
  50 //IplImage values of image that are useful
  51 //int  nSize;         /* sizeof(IplImage) */
  52 //int  depth;         /* pixel depth in bits: IPL_DEPTH_8U ...*/
  53 //int  nChannels;     /* OpenCV functions support 1,2,3 or 4 channels */
  54 //int  width;         /* image width in pixels */
  55 //int  height;        /* image height in pixels */
  56 //int  imageSize;     /* image data size in bytes in case of interleaved data)*/
  57 //char *imageData;    /* pointer to aligned image data */
  58 //char *imageDataOrigin; /* pointer to very origin of image -deallocation */
  59 //Values useful for gaussian integral
  60 //0.5 - 0.19146 - 0.38292
  61 //1.0 - 0.34134 - 0.68268
  62 //1.5 - 0.43319 - 0.86638
  63 //2.0 - 0.47725 - 0.95450
  64 //2.5 - 0.49379 - 0.98758
  65 //3.0 - 0.49865 - 0.99730
  66 //3.5 - 0.4997674 - 0.9995348
  67 //4.0 - 0.4999683 - 0.9999366
  68
  69 #include "_cvaux.h"
  70
  71
  72 //internal functions for gaussian background detection
  73 static void icvInsertionSortGaussians( CvGaussBGPoint* g_point, double* sort_key, CvGaussBGStatModelParams *bg_model_params );
  74
  75 /*
  76    Test whether pixel can be explained by background model;
  77    Return -1 if no match was found; otherwise the index in match[] is returned
  78
  79    icvMatchTest(...) assumes what all color channels component exhibit the same variance
  80    icvMatchTest2(...) accounts for different variances per color channel
  81  */
  82 static int icvMatchTest( double* src_pixel, int nChannels, int* match,
  83                  const CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params );
  84 /*static int icvMatchTest2( double* src_pixel, int nChannels, int* match,
  85                  const CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params );*/
  86
  87
  88 /*
  89    The update procedure differs between
  90       * the initialization phase (named *Partial* ) and
  91       * the normal phase (named *Full* )
  92    The initalization phase is defined as not having processed <win_size> frames yet
  93  */
  94 static void icvUpdateFullWindow( double* src_pixel, int nChannels,
  95                          int* match,
  96                          CvGaussBGPoint* g_point,
  97                          const CvGaussBGStatModelParams *bg_model_params );
  98 static void icvUpdateFullNoMatch( IplImage* gm_image, int p,
  99                           int* match,
 100                           CvGaussBGPoint* g_point,
 101                           const CvGaussBGStatModelParams *bg_model_params);
 102 static void icvUpdatePartialWindow( double* src_pixel, int nChannels, int* match,
 103                             CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params );
 104 static void icvUpdatePartialNoMatch( double* src_pixel, int nChannels,
 105                              int* match,
 106                              CvGaussBGPoint* g_point,
 107                              const CvGaussBGStatModelParams *bg_model_params);
 108
 109
 110 static void icvGetSortKey( const int nChannels, double* sort_key, const CvGaussBGPoint* g_point,
 111                     const CvGaussBGStatModelParams *bg_model_params );
 112 static void icvBackgroundTest( const int nChannels, int n, int i, int j, int *match, CvGaussBGModel* bg_model );
 113
 114 static void CV_CDECL icvReleaseGaussianBGModel( CvGaussBGModel** bg_model );
 115 static int CV_CDECL icvUpdateGaussianBGModel( IplImage* curr_frame, CvGaussBGModel*  bg_model );
 116
 117 //#define for if(0);else for
 118
 119 //g = 1 for first gaussian in list that matches else g = 0
 120 //Rw is the learning rate for weight and Rg is leaning rate for mean and variance
 121 //Ms is the match_sum which is the sum of matches for a particular gaussian
 122 //Ms values are incremented until the sum of Ms values in the list equals window size L
 123 //SMs is the sum of match_sums for gaussians in the list
 124 //Rw = 1/SMs note the smallest Rw gets is 1/L
 125 //Rg = g/Ms for SMs < L and Rg = g/(w*L) for SMs = L
 126 //The list is maintained in sorted order using w/sqrt(variance) as a key
 127 //If there is no match the last gaussian in the list is replaced by the new gaussian
 128 //This will result in changes to SMs which results in changes in Rw and Rg.
 129 //If a gaussian is replaced and SMs previously equaled L values of Ms are computed from w
 130 //w[n+1] = w[n] + Rw*(g - w[n])   weight
 131 //u[n+1] = u[n] + Rg*(x[n+1] - u[n]) mean value Sg is sum n values of g
 132 //v[n+1] = v[n] + Rg*((x[n+1] - u[n])*(x[n+1] - u[n])) - v[n]) variance
 133 //
 134
 135 CV_IMPL CvBGStatModel*
 136 cvCreateGaussianBGModel( IplImage* first_frame, CvGaussBGStatModelParams* parameters )
 137 {
 138     CvGaussBGModel* bg_model = 0;
 139
 140     CV_FUNCNAME( "cvCreateGaussianBGModel" );
 141
 142     __BEGIN__;
 143
 144     double var_init;
 145     CvGaussBGStatModelParams params;
 146     int i, j, k, m, n;
 147
 148     //init parameters
 149     if( parameters == NULL )
 150       {                        /* These constants are defined in cvaux/include/cvaux.h: */
 151         params.win_size      = CV_BGFG_MOG_WINDOW_SIZE;
 152         params.bg_threshold  = CV_BGFG_MOG_BACKGROUND_THRESHOLD;
 153
 154         params.std_threshold = CV_BGFG_MOG_STD_THRESHOLD;
 155         params.weight_init   = CV_BGFG_MOG_WEIGHT_INIT;
 156
 157         params.variance_init = CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT;
 158         params.minArea       = CV_BGFG_MOG_MINAREA;
 159         params.n_gauss       = CV_BGFG_MOG_NGAUSSIANS;
 160     }
 161     else
 162     {
 163         params = *parameters;
 164     }
 165
 166     if( !CV_IS_IMAGE(first_frame) )
 167         CV_ERROR( CV_StsBadArg, "Invalid or NULL first_frame parameter" );
 168
 169     CV_CALL( bg_model = (CvGaussBGModel*)cvAlloc( sizeof(*bg_model) ));
 170     memset( bg_model, 0, sizeof(*bg_model) );
 171     bg_model->type = CV_BG_MODEL_MOG;
 172     bg_model->release = (CvReleaseBGStatModel)icvReleaseGaussianBGModel;
 173     bg_model->update = (CvUpdateBGStatModel)icvUpdateGaussianBGModel;
 174
 175     bg_model->params = params;
 176
 177     //prepare storages
 178     CV_CALL( bg_model->g_point = (CvGaussBGPoint*)cvAlloc(sizeof(CvGaussBGPoint)*
 179         ((first_frame->width*first_frame->height) + 256)));
 180
 181     CV_CALL( bg_model->background = cvCreateImage(cvSize(first_frame->width,
 182         first_frame->height), IPL_DEPTH_8U, first_frame->nChannels));
 183     CV_CALL( bg_model->foreground = cvCreateImage(cvSize(first_frame->width,
 184         first_frame->height), IPL_DEPTH_8U, 1));
 185
 186     CV_CALL( bg_model->storage = cvCreateMemStorage());
 187
 188     //initializing
 189     var_init = 2 * params.std_threshold * params.std_threshold;
 190     CV_CALL( bg_model->g_point[0].g_values =
 191         (CvGaussBGValues*)cvAlloc( sizeof(CvGaussBGValues)*params.n_gauss*
 192         (first_frame->width*first_frame->height + 128)));
 193
 194     for( i = 0, n = 0; i < first_frame->height; i++ )
 195     {
 196         for( j = 0; j < first_frame->width; j++, n++ )
 197         {
 198             const int p = i*first_frame->widthStep+j*first_frame->nChannels;
 199
 200             bg_model->g_point[n].g_values =
 201                 bg_model->g_point[0].g_values + n*params.n_gauss;
 202             bg_model->g_point[n].g_values[0].weight = 1;    //the first value seen has weight one
 203             bg_model->g_point[n].g_values[0].match_sum = 1;
 204             for( m = 0; m < first_frame->nChannels; m++)
 205             {
 206                 bg_model->g_point[n].g_values[0].variance[m] = var_init;
 207                 bg_model->g_point[n].g_values[0].mean[m] = (unsigned char)first_frame->imageData[p + m];
 208             }
 209             for( k = 1; k < params.n_gauss; k++)
 210             {
 211                 bg_model->g_point[n].g_values[k].weight = 0;
 212                 bg_model->g_point[n].g_values[k].match_sum = 0;
 213                 for( m = 0; m < first_frame->nChannels; m++){
 214                     bg_model->g_point[n].g_values[k].variance[m] = var_init;
 215                     bg_model->g_point[n].g_values[k].mean[m] = 0;
 216                 }
 217             }
 218         }
 219     }
 220
 221     bg_model->countFrames = 0;
 222
 223     __END__;
 224
 225     if( cvGetErrStatus() < 0 )
 226     {
 227         CvBGStatModel* base_ptr = (CvBGStatModel*)bg_model;
 228
 229         if( bg_model && bg_model->release )
 230             bg_model->release( &base_ptr );
 231         else
 232             cvFree( &bg_model );
 233         bg_model = 0;
 234     }
 235
 236     return (CvBGStatModel*)bg_model;
 237 }
 238
 239
 240 static void CV_CDECL
 241 icvReleaseGaussianBGModel( CvGaussBGModel** _bg_model )
 242 {
 243     CV_FUNCNAME( "icvReleaseGaussianBGModel" );
 244
 245     __BEGIN__;
 246
 247     if( !_bg_model )
 248         CV_ERROR( CV_StsNullPtr, "" );
 249
 250     if( *_bg_model )
 251     {
 252         CvGaussBGModel* bg_model = *_bg_model;
 253         if( bg_model->g_point )
 254         {
 255             cvFree( &bg_model->g_point[0].g_values );
 256             cvFree( &bg_model->g_point );
 257         }
 258
 259         cvReleaseImage( &bg_model->background );
 260         cvReleaseImage( &bg_model->foreground );
 261         cvReleaseMemStorage(&bg_model->storage);
 262         memset( bg_model, 0, sizeof(*bg_model) );
 263         cvFree( _bg_model );
 264     }
 265
 266     __END__;
 267 }
 268
 269
 270 static int CV_CDECL
 271 icvUpdateGaussianBGModel( IplImage* curr_frame, CvGaussBGModel*  bg_model )
 272 {
 273     int i, j, k, n;
 274     int region_count = 0;
 275     CvSeq *first_seq = NULL, *prev_seq = NULL, *seq = NULL;
 276
 277     bg_model->countFrames++;
 278
 279     for( i = 0, n = 0; i < curr_frame->height; i++ )
 280     {
 281         for( j = 0; j < curr_frame->width; j++, n++ )
 282         {
 283             int match[CV_BGFG_MOG_MAX_NGAUSSIANS];
 284             double sort_key[CV_BGFG_MOG_MAX_NGAUSSIANS];
 285             const int nChannels = curr_frame->nChannels;
 286             const int p = curr_frame->widthStep*i+j*nChannels;
 287
 288             // A few short cuts
 289             CvGaussBGPoint* g_point = &bg_model->g_point[n];
 290             const CvGaussBGStatModelParams bg_model_params = bg_model->params;
 291             double pixel[4];
 292             int no_match;
 293
 294             for( k = 0; k < nChannels; k++ )
 295                 pixel[k] = (uchar)curr_frame->imageData[p+k];
 296
 297             no_match = icvMatchTest( pixel, nChannels, match, g_point, &bg_model_params );
 298             if( bg_model->countFrames >= bg_model->params.win_size )
 299             {
 300                 icvUpdateFullWindow( pixel, nChannels, match, g_point, &bg_model->params );
 301                 if( no_match == -1)
 302                     icvUpdateFullNoMatch( curr_frame, p, match, g_point, &bg_model_params );
 303             }
 304             else
 305             {
 306                 icvUpdatePartialWindow( pixel, nChannels, match, g_point, &bg_model_params );
 307                 if( no_match == -1)
 308                     icvUpdatePartialNoMatch( pixel, nChannels, match, g_point, &bg_model_params );
 309             }
 310             icvGetSortKey( nChannels, sort_key, g_point, &bg_model_params );
 311             icvInsertionSortGaussians( g_point, sort_key, (CvGaussBGStatModelParams *)&bg_model_params );
 312             icvBackgroundTest( nChannels, n, i, j, match, bg_model );
 313         }
 314     }
 315
 316     //foreground filtering
 317
 318     //filter small regions
 319     cvClearMemStorage(bg_model->storage);
 320
 321     //cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_OPEN, 1 );
 322     //cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_CLOSE, 1 );
 323
 324     cvFindContours( bg_model->foreground, bg_model->storage, &first_seq, sizeof(CvContour), CV_RETR_LIST );
 325     for( seq = first_seq; seq; seq = seq->h_next )
 326     {
 327         CvContour* cnt = (CvContour*)seq;
 328         if( cnt->rect.width * cnt->rect.height < bg_model->params.minArea )
 329         {
 330             //delete small contour
 331             prev_seq = seq->h_prev;
 332             if( prev_seq )
 333             {
 334                 prev_seq->h_next = seq->h_next;
 335                 if( seq->h_next ) seq->h_next->h_prev = prev_seq;
 336             }
 337             else
 338             {
 339                 first_seq = seq->h_next;
 340                 if( seq->h_next ) seq->h_next->h_prev = NULL;
 341             }
 342         }
 343         else
 344         {
 345             region_count++;
 346         }
 347     }
 348     bg_model->foreground_regions = first_seq;
 349     cvZero(bg_model->foreground);
 350     cvDrawContours(bg_model->foreground, first_seq, CV_RGB(0, 0, 255), CV_RGB(0, 0, 255), 10, -1);
 351
 352     return region_count;
 353 }
 354
 355 static void icvInsertionSortGaussians( CvGaussBGPoint* g_point, double* sort_key, CvGaussBGStatModelParams *bg_model_params )
 356 {
 357     int i, j;
 358     for( i = 1; i < bg_model_params->n_gauss; i++ )
 359     {
 360         double index = sort_key[i];
 361         for( j = i; j > 0 && sort_key[j-1] < index; j-- ) //sort decending order
 362         {
 363             double temp_sort_key = sort_key[j];
 364             sort_key[j] = sort_key[j-1];
 365             sort_key[j-1] = temp_sort_key;
 366
 367             CvGaussBGValues temp_gauss_values = g_point->g_values[j];
 368             g_point->g_values[j] = g_point->g_values[j-1];
 369             g_point->g_values[j-1] = temp_gauss_values;
 370         }
 371 //        sort_key[j] = index;
 372     }
 373 }
 374
 375
 376 static int icvMatchTest( double* src_pixel, int nChannels, int* match,
 377                          const CvGaussBGPoint* g_point,
 378                          const CvGaussBGStatModelParams *bg_model_params )
 379 {
 380     int k;
 381     int matchPosition=-1;
 382     for ( k = 0; k < bg_model_params->n_gauss; k++) match[k]=0;
 383
 384     for ( k = 0; k < bg_model_params->n_gauss; k++) {
 385         double sum_d2 = 0.0;
 386         double var_threshold = 0.0;
 387         for(int m = 0; m < nChannels; m++){
 388             double d = g_point->g_values[k].mean[m]- src_pixel[m];
 389             sum_d2 += (d*d);
 390             var_threshold += g_point->g_values[k].variance[m];
 391         }  //difference < STD_LIMIT*STD_LIMIT or difference**2 < STD_LIMIT*STD_LIMIT*VAR
 392         var_threshold = bg_model_params->std_threshold*bg_model_params->std_threshold*var_threshold;
 393         if(sum_d2 < var_threshold){
 394             match[k] = 1;
 395             matchPosition = k;
 396             break;
 397         }
 398     }
 399
 400     return matchPosition;
 401 }
 402
 403 /*
 404 static int icvMatchTest2( double* src_pixel, int nChannels, int* match,
 405                           const CvGaussBGPoint* g_point,
 406                           const CvGaussBGStatModelParams *bg_model_params )
 407 {
 408     int k, m;
 409     int matchPosition=-1;
 410
 411     for( k = 0; k < bg_model_params->n_gauss; k++ )
 412         match[k] = 0;
 413
 414     for( k = 0; k < bg_model_params->n_gauss; k++ )
 415     {
 416         double sum_d2 = 0.0, var_threshold;
 417         for( m = 0; m < nChannels; m++ )
 418         {
 419             double d = g_point->g_values[k].mean[m]- src_pixel[m];
 420             sum_d2 += (d*d) / (g_point->g_values[k].variance[m] * g_point->g_values[k].variance[m]);
 421         }  //difference < STD_LIMIT*STD_LIMIT or difference**2 < STD_LIMIT*STD_LIMIT*VAR
 422
 423         var_threshold = bg_model_params->std_threshold*bg_model_params->std_threshold;
 424         if( sum_d2 < var_threshold )
 425         {
 426             match[k] = 1;
 427             matchPosition = k;
 428             break;
 429         }
 430     }
 431
 432     return matchPosition;
 433 }
 434 */
 435
 436 static void icvUpdateFullWindow( double* src_pixel, int nChannels, int* match,
 437                                  CvGaussBGPoint* g_point,
 438                                  const CvGaussBGStatModelParams *bg_model_params )
 439 {
 440     const double learning_rate_weight = (1.0/(double)bg_model_params->win_size);
 441     for(int k = 0; k < bg_model_params->n_gauss; k++){
 442         g_point->g_values[k].weight = g_point->g_values[k].weight +
 443             (learning_rate_weight*((double)match[k] -
 444             g_point->g_values[k].weight));
 445         if(match[k]){
 446             double learning_rate_gaussian = (double)match[k]/(g_point->g_values[k].weight*
 447                 (double)bg_model_params->win_size);
 448             for(int m = 0; m < nChannels; m++){
 449                 const double tmpDiff = src_pixel[m] - g_point->g_values[k].mean[m];
 450                 g_point->g_values[k].mean[m] = g_point->g_values[k].mean[m] +
 451                     (learning_rate_gaussian * tmpDiff);
 452                 g_point->g_values[k].variance[m] = g_point->g_values[k].variance[m]+
 453                     (learning_rate_gaussian*((tmpDiff*tmpDiff) - g_point->g_values[k].variance[m]));
 454             }
 455         }
 456     }
 457 }
 458
 459
 460 static void icvUpdatePartialWindow( double* src_pixel, int nChannels, int* match, CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params )
 461 {
 462     int k, m;
 463     int window_current = 0;
 464
 465     for( k = 0; k < bg_model_params->n_gauss; k++ )
 466         window_current += g_point->g_values[k].match_sum;
 467
 468     for( k = 0; k < bg_model_params->n_gauss; k++ )
 469     {
 470         g_point->g_values[k].match_sum += match[k];
 471         double learning_rate_weight = (1.0/((double)window_current + 1.0)); //increased by one since sum
 472         g_point->g_values[k].weight = g_point->g_values[k].weight +
 473             (learning_rate_weight*((double)match[k] - g_point->g_values[k].weight));
 474
 475         if( g_point->g_values[k].match_sum > 0 && match[k] )
 476         {
 477             double learning_rate_gaussian = (double)match[k]/((double)g_point->g_values[k].match_sum);
 478             for( m = 0; m < nChannels; m++ )
 479             {
 480                 const double tmpDiff = src_pixel[m] - g_point->g_values[k].mean[m];
 481                 g_point->g_values[k].mean[m] = g_point->g_values[k].mean[m] +
 482                     (learning_rate_gaussian*tmpDiff);
 483                 g_point->g_values[k].variance[m] = g_point->g_values[k].variance[m]+
 484                     (learning_rate_gaussian*((tmpDiff*tmpDiff) - g_point->g_values[k].variance[m]));
 485             }
 486         }
 487     }
 488 }
 489
 490 static void icvUpdateFullNoMatch( IplImage* gm_image, int p, int* match,
 491                                   CvGaussBGPoint* g_point,
 492                                   const CvGaussBGStatModelParams *bg_model_params)
 493 {
 494     int k, m;
 495     double alpha;
 496     int match_sum_total = 0;
 497
 498     //new value of last one
 499     g_point->g_values[bg_model_params->n_gauss - 1].match_sum = 1;
 500
 501     //get sum of all but last value of match_sum
 502
 503     for( k = 0; k < bg_model_params->n_gauss ; k++ )
 504         match_sum_total += g_point->g_values[k].match_sum;
 505
 506     g_point->g_values[bg_model_params->n_gauss - 1].weight = 1./(double)match_sum_total;
 507     for( m = 0; m < gm_image->nChannels ; m++ )
 508     {
 509         // first pass mean is image value
 510         g_point->g_values[bg_model_params->n_gauss - 1].variance[m] = bg_model_params->variance_init;
 511         g_point->g_values[bg_model_params->n_gauss - 1].mean[m] = (unsigned char)gm_image->imageData[p + m];
 512     }
 513
 514     alpha = 1.0 - (1.0/bg_model_params->win_size);
 515     for( k = 0; k < bg_model_params->n_gauss - 1; k++ )
 516     {
 517         g_point->g_values[k].weight *= alpha;
 518         if( match[k] )
 519             g_point->g_values[k].weight += alpha;
 520     }
 521 }
 522
 523
 524 static void
 525 icvUpdatePartialNoMatch(double *pixel,
 526                         int nChannels,
 527                         int* /*match*/,
 528                         CvGaussBGPoint* g_point,
 529                         const CvGaussBGStatModelParams *bg_model_params)
 530 {
 531     int k, m;
 532     //new value of last one
 533     g_point->g_values[bg_model_params->n_gauss - 1].match_sum = 1;
 534
 535     //get sum of all but last value of match_sum
 536     int match_sum_total = 0;
 537     for(k = 0; k < bg_model_params->n_gauss ; k++)
 538         match_sum_total += g_point->g_values[k].match_sum;
 539
 540     for(m = 0; m < nChannels; m++)
 541     {
 542         //first pass mean is image value
 543         g_point->g_values[bg_model_params->n_gauss - 1].variance[m] = bg_model_params->variance_init;
 544         g_point->g_values[bg_model_params->n_gauss - 1].mean[m] = pixel[m];
 545     }
 546     for(k = 0; k < bg_model_params->n_gauss; k++)
 547     {
 548         g_point->g_values[k].weight = (double)g_point->g_values[k].match_sum /
 549             (double)match_sum_total;
 550     }
 551 }
 552
 553 static void icvGetSortKey( const int nChannels, double* sort_key, const CvGaussBGPoint* g_point,
 554                            const CvGaussBGStatModelParams *bg_model_params )
 555 {
 556     int k, m;
 557     for( k = 0; k < bg_model_params->n_gauss; k++ )
 558     {
 559         // Avoid division by zero
 560         if( g_point->g_values[k].match_sum > 0 )
 561         {
 562             // Independence assumption between components
 563             double variance_sum = 0.0;
 564             for( m = 0; m < nChannels; m++ )
 565                 variance_sum += g_point->g_values[k].variance[m];
 566
 567             sort_key[k] = g_point->g_values[k].weight/sqrt(variance_sum);
 568         }
 569         else
 570             sort_key[k]= 0.0;
 571     }
 572 }
 573
 574
 575 static void icvBackgroundTest( const int nChannels, int n, int i, int j, int *match, CvGaussBGModel* bg_model )
 576 {
 577     int m, b;
 578     uchar pixelValue = (uchar)255; // will switch to 0 if match found
 579     double weight_sum = 0.0;
 580     CvGaussBGPoint* g_point = bg_model->g_point;
 581
 582     for( m = 0; m < nChannels; m++)
 583         bg_model->background->imageData[ bg_model->background->widthStep*i + j*nChannels + m]  = (unsigned char)(g_point[n].g_values[0].mean[m]+0.5);
 584
 585     for( b = 0; b < bg_model->params.n_gauss; b++)
 586     {
 587         weight_sum += g_point[n].g_values[b].weight;
 588         if( match[b] )
 589             pixelValue = 0;
 590         if( weight_sum > bg_model->params.bg_threshold )
 591             break;
 592     }
 593
 594     bg_model->foreground->imageData[ bg_model->foreground->widthStep*i + j] = pixelValue;
 595 }
 596
 597 /* End of file. */