vcs.maemo.org Git - opencv/blob - cvaux/src/cvbgfg_acmmm2003.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
  43 #include <math.h>
  44 #include <stdio.h>
  45 #include <stdlib.h>
  46 //#include <algorithm>
  47
  48 static double* _cv_max_element( double* start, double* end )
  49 {
  50     double* p = start++;
  51     for( ; start != end; start++ )
  52         if( *p < *start )
  53             p = start;
  54     return p;
  55 }
  56
  57 static void CV_CDECL icvReleaseFGDStatModel( CvFGDStatModel** model );
  58 static int CV_CDECL icvUpdateFGDStatModel( IplImage* curr_frame,
  59                                            CvFGDStatModel*  model );
  60
  61 //  Function cvCreateFGDStatModel initializes foreground detection process
  62 // parameters:
  63 //      first_frame - frame from video sequence
  64 //      parameters  - (optional) if NULL default parameters of the algorithm will be used
  65 //      p_model     - pointer to CvFGDStatModel structure
  66 CV_IMPL CvBGStatModel*
  67 cvCreateFGDStatModel( IplImage* first_frame, CvFGDStatModelParams* parameters )
  68 {
  69     CvFGDStatModel* p_model = 0;
  70
  71     CV_FUNCNAME( "cvCreateFGDStatModel" );
  72
  73     __BEGIN__;
  74
  75     int i, j, k, pixel_count, buf_size;
  76     CvFGDStatModelParams params;
  77
  78     if( !CV_IS_IMAGE(first_frame) )
  79         CV_ERROR( CV_StsBadArg, "Invalid or NULL first_frame parameter" );
  80
  81     if (first_frame->nChannels != 3)
  82         CV_ERROR( CV_StsBadArg, "first_frame must have 3 color channels" );
  83
  84     //init parameters
  85     if( parameters == NULL )
  86     {
  87         params.Lc = CV_BGFG_FGD_LC;
  88         params.N1c = CV_BGFG_FGD_N1C;
  89         params.N2c = CV_BGFG_FGD_N2C;
  90         params.Lcc = CV_BGFG_FGD_LCC;
  91         params.N1cc = CV_BGFG_FGD_N1CC;
  92         params.N2cc = CV_BGFG_FGD_N2CC;
  93         params.delta = CV_BGFG_FGD_DELTA;
  94         params.alpha1 = CV_BGFG_FGD_ALPHA_1;
  95         params.alpha2 = CV_BGFG_FGD_ALPHA_2;
  96         params.alpha3 = CV_BGFG_FGD_ALPHA_3;
  97         params.T = CV_BGFG_FGD_T;
  98         params.minArea = CV_BGFG_FGD_MINAREA;
  99         params.is_obj_without_holes = 1;
 100         params.perform_morphing = 1;
 101     }
 102     else
 103     {
 104         params = *parameters;
 105     }
 106
 107     CV_CALL( p_model = (CvFGDStatModel*)cvAlloc( sizeof(*p_model) ));
 108     memset( p_model, 0, sizeof(*p_model) );
 109     p_model->type = CV_BG_MODEL_FGD;
 110     p_model->release = (CvReleaseBGStatModel)icvReleaseFGDStatModel;
 111     p_model->update = (CvUpdateBGStatModel)icvUpdateFGDStatModel;;
 112     p_model->params = params;
 113
 114     //init storages
 115     pixel_count = first_frame->width*first_frame->height;
 116
 117     buf_size = pixel_count*sizeof(p_model->pixel_stat[0]);
 118     CV_CALL( p_model->pixel_stat = (CvBGPixelStat*)cvAlloc(buf_size) );
 119     memset( p_model->pixel_stat, 0, buf_size );
 120
 121     buf_size = pixel_count*params.N2c*sizeof(p_model->pixel_stat[0].ctable[0]);
 122     CV_CALL( p_model->pixel_stat[0].ctable = (CvBGPixelCStatTable*)cvAlloc(buf_size) );
 123     memset( p_model->pixel_stat[0].ctable, 0, buf_size );
 124
 125     buf_size = pixel_count*params.N2cc*sizeof(p_model->pixel_stat[0].cctable[0]);
 126     CV_CALL( p_model->pixel_stat[0].cctable = (CvBGPixelCCStatTable*)cvAlloc(buf_size) );
 127     memset( p_model->pixel_stat[0].cctable, 0, buf_size );
 128
 129     for( i = 0, k = 0; i < first_frame->height; i++ )
 130         for( j = 0; j < first_frame->width; j++, k++ )
 131         {
 132             p_model->pixel_stat[k].ctable = p_model->pixel_stat[0].ctable + k*params.N2c;
 133             p_model->pixel_stat[k].cctable = p_model->pixel_stat[0].cctable + k*params.N2cc;
 134         }
 135
 136     //init temporary images
 137     CV_CALL( p_model->Ftd = cvCreateImage(cvSize(first_frame->width, first_frame->height), IPL_DEPTH_8U, 1));
 138     CV_CALL( p_model->Fbd = cvCreateImage(cvSize(first_frame->width, first_frame->height), IPL_DEPTH_8U, 1));
 139     CV_CALL( p_model->foreground = cvCreateImage(cvSize(first_frame->width, first_frame->height), IPL_DEPTH_8U, 1));
 140
 141         CV_CALL( p_model->background = cvCloneImage(first_frame));
 142     CV_CALL( p_model->prev_frame = cvCloneImage(first_frame));
 143     CV_CALL( p_model->storage = cvCreateMemStorage());
 144
 145     __END__;
 146
 147     if( cvGetErrStatus() < 0 )
 148     {
 149         CvBGStatModel* base_ptr = (CvBGStatModel*)p_model;
 150
 151         if( p_model && p_model->release )
 152             p_model->release( &base_ptr );
 153         else
 154             cvFree( &p_model );
 155         p_model = 0;
 156     }
 157
 158     return (CvBGStatModel*)p_model;
 159 }
 160
 161
 162 static void CV_CDECL
 163 icvReleaseFGDStatModel( CvFGDStatModel** _model )
 164 {
 165     CV_FUNCNAME( "icvReleaseFGDStatModel" );
 166
 167     __BEGIN__;
 168
 169     if( !_model )
 170         CV_ERROR( CV_StsNullPtr, "" );
 171
 172     if( *_model )
 173     {
 174         CvFGDStatModel* model = *_model;
 175         if( model->pixel_stat )
 176         {
 177             cvFree( &model->pixel_stat[0].ctable );
 178             cvFree( &model->pixel_stat[0].cctable );
 179             cvFree( &model->pixel_stat );
 180         }
 181
 182         cvReleaseImage( &model->Ftd );
 183         cvReleaseImage( &model->Fbd );
 184         cvReleaseImage( &model->foreground );
 185         cvReleaseImage( &model->background );
 186         cvReleaseImage( &model->prev_frame );
 187         cvReleaseMemStorage(&model->storage);
 188
 189         cvFree( _model );
 190     }
 191
 192     __END__;
 193 }
 194
 195 //  Function cvChangeDetection performs change detection for Foreground detection algorithm
 196 // parameters:
 197 //      prev_frame -
 198 //      curr_frame -
 199 //      change_mask -
 200 CV_IMPL int
 201 cvChangeDetection( IplImage*  prev_frame,
 202                    IplImage*  curr_frame,
 203                    IplImage*  change_mask )
 204 {
 205     int i, j, b, x, y, thres;
 206     const int PIXELRANGE=256;
 207
 208     if( !prev_frame || !curr_frame || !change_mask ||
 209         prev_frame->nChannels != 3 || curr_frame->nChannels != 3 || change_mask->nChannels != 1 ||
 210         prev_frame->depth != IPL_DEPTH_8U || curr_frame->depth != IPL_DEPTH_8U || change_mask->depth != IPL_DEPTH_8U ||
 211         !CV_ARE_SIZES_EQ( prev_frame, curr_frame ) || !CV_ARE_SIZES_EQ( prev_frame, change_mask ) ) return 0;
 212
 213     cvZero ( change_mask );
 214
 215     // All operations per colour
 216     for (b=0 ; b<prev_frame->nChannels ; b++) {
 217         // create histogram
 218
 219         long HISTOGRAM[PIXELRANGE];
 220         for (i=0 ; i<PIXELRANGE; i++) HISTOGRAM[i]=0;
 221
 222         for (y=0 ; y<curr_frame->height ; y++)
 223         {
 224             uchar* rowStart1 = (uchar*)curr_frame->imageData + y * curr_frame->widthStep + b;
 225             uchar* rowStart2 = (uchar*)prev_frame->imageData + y * prev_frame->widthStep + b;
 226             for (x=0 ; x<curr_frame->width ; x++, rowStart1+=curr_frame->nChannels, rowStart2+=prev_frame->nChannels) {
 227                 int diff = abs( int(*rowStart1) - int(*rowStart2) );
 228                 HISTOGRAM[diff]++;
 229             }
 230         }
 231
 232         double relativeVariance[PIXELRANGE];
 233         for (i=0 ; i<PIXELRANGE; i++) relativeVariance[i]=0;
 234
 235         for (thres=PIXELRANGE-2; thres>=0 ; thres--)
 236         {
 237             //            fprintf(stderr, "Iter %d\n", thres);
 238             double sum=0;
 239             double sqsum=0;
 240             int count=0;
 241             //            fprintf(stderr, "Iter %d entering loop\n", thres);
 242             for (j=thres ; j<PIXELRANGE ; j++) {
 243                 sum   += double(j)*double(HISTOGRAM[j]);
 244                 sqsum += double(j*j)*double(HISTOGRAM[j]);
 245                 count += HISTOGRAM[j];
 246             }
 247             count = count == 0 ? 1 : count;
 248             //            fprintf(stderr, "Iter %d finishing loop\n", thres);
 249             double my = sum / count;
 250             double sigma = sqrt( sqsum/count - my*my);
 251             //            fprintf(stderr, "Iter %d sum=%g sqsum=%g count=%d sigma = %g\n", thres, sum, sqsum, count, sigma);
 252             //            fprintf(stderr, "Writing to %x\n", &(relativeVariance[thres]));
 253             relativeVariance[thres] = sigma;
 254             //            fprintf(stderr, "Iter %d finished\n", thres);
 255         }
 256         // find maximum
 257         uchar bestThres = 0;
 258
 259         double* pBestThres = _cv_max_element(relativeVariance, relativeVariance+PIXELRANGE);
 260         bestThres = (uchar)(*pBestThres); if (bestThres <10) bestThres=10;
 261
 262         for (y=0 ; y<prev_frame->height ; y++)
 263         {
 264             uchar* rowStart1 = (uchar*)(curr_frame->imageData) + y * curr_frame->widthStep + b;
 265             uchar* rowStart2 = (uchar*)(prev_frame->imageData) + y * prev_frame->widthStep + b;
 266             uchar* rowStart3 = (uchar*)(change_mask->imageData) + y * change_mask->widthStep;
 267             for (x = 0; x < curr_frame->width; x++, rowStart1+=curr_frame->nChannels,
 268                 rowStart2+=prev_frame->nChannels, rowStart3+=change_mask->nChannels) {
 269                 // OR between different color channels
 270                 int diff = abs( int(*rowStart1) - int(*rowStart2) );
 271                 if ( diff > bestThres)
 272                     *rowStart3 |=255;
 273             }
 274         }
 275     }
 276
 277     return 1;
 278 }
 279
 280
 281 #define MIN_PV 1E-10
 282
 283
 284 #define V_C(k,l) ctable[k].v[l]
 285 #define PV_C(k) ctable[k].Pv
 286 #define PVB_C(k) ctable[k].Pvb
 287 #define V_CC(k,l) cctable[k].v[l]
 288 #define PV_CC(k) cctable[k].Pv
 289 #define PVB_CC(k) cctable[k].Pvb
 290
 291
 292 //  Function cvUpdateFGDStatModel updates statistical model and returns number of foreground regions
 293 // parameters:
 294 //      curr_frame  - current frame from video sequence
 295 //      p_model     - pointer to CvFGDStatModel structure
 296 static int CV_CDECL
 297 icvUpdateFGDStatModel( IplImage* curr_frame, CvFGDStatModel*  model )
 298 {
 299     int            mask_step = model->Ftd->widthStep;
 300     CvSeq         *first_seq = NULL, *prev_seq = NULL, *seq = NULL;
 301     IplImage*      prev_frame = model->prev_frame;
 302     int            region_count = 0;
 303     int            FG_pixels_count = 0;
 304     int            deltaC  = cvRound(model->params.delta * 256 / model->params.Lc);
 305     int            deltaCC = cvRound(model->params.delta * 256 / model->params.Lcc);
 306     int            i, j, k, l;
 307
 308     //clear storages
 309     cvClearMemStorage(model->storage);
 310     cvZero(model->foreground);
 311
 312     //form FG pixels candidates using image differencing with adaptive threshold [P.Rosin, Thresholding for change detection, ICCV, 1998 ]
 313     cvChangeDetection( prev_frame, curr_frame, model->Ftd );
 314     cvChangeDetection( model->background, curr_frame, model->Fbd );
 315
 316     for( i = 0; i < model->Ftd->height; i++ )
 317     {
 318         for( j = 0; j < model->Ftd->width; j++ )
 319         {
 320             if( ((uchar*)model->Fbd->imageData)[i*mask_step+j] || ((uchar*)model->Ftd->imageData)[i*mask_step+j] )
 321             {
 322                 float Pb=0, Pv=0, Pvb=0;
 323                 CvBGPixelStat* stat = model->pixel_stat + i * model->Ftd->width + j;
 324                 CvBGPixelCStatTable* ctable = stat->ctable;
 325                 CvBGPixelCCStatTable* cctable = stat->cctable;
 326
 327                 uchar* curr_data = (uchar*)(curr_frame->imageData)+i*curr_frame->widthStep+j*3;
 328                 uchar* prev_data = (uchar*)(prev_frame->imageData)+i*prev_frame->widthStep+j*3;
 329
 330                 int val = 0;
 331                 // is it a motion pixel?
 332                 if( ((uchar*)model->Ftd->imageData)[i*mask_step+j] )
 333                 {
 334                     if( !stat->is_trained_dyn_model ) val = 1;
 335                     else
 336                     {
 337                         //compare with stored CCt vectors
 338                         for( k = 0; PV_CC(k) > model->params.alpha2 && k < model->params.N1cc; k++ )
 339                         {
 340                             if ( abs( V_CC(k,0) - prev_data[0]) <=  deltaCC &&
 341                                  abs( V_CC(k,1) - prev_data[1]) <=  deltaCC &&
 342                                  abs( V_CC(k,2) - prev_data[2]) <=  deltaCC &&
 343                                  abs( V_CC(k,3) - curr_data[0]) <=  deltaCC &&
 344                                  abs( V_CC(k,4) - curr_data[1]) <=  deltaCC &&
 345                                  abs( V_CC(k,5) - curr_data[2]) <=  deltaCC)
 346                             {
 347                                 Pv += PV_CC(k);
 348                                 Pvb += PVB_CC(k);
 349                             }
 350                         }
 351                         Pb = stat->Pbcc;
 352                         if( 2 * Pvb * Pb <= Pv ) val = 1;
 353                     }
 354                 }
 355                 else if( stat->is_trained_st_model )
 356                 {
 357                     //compare with stored Ct vectors
 358                     for( k = 0; PV_C(k) > model->params.alpha2 && k < model->params.N1c; k++ )
 359                     {
 360                         if ( abs( V_C(k,0) - curr_data[0]) <=  deltaC &&
 361                              abs( V_C(k,1) - curr_data[1]) <=  deltaC &&
 362                              abs( V_C(k,2) - curr_data[2]) <=  deltaC )
 363                         {
 364                             Pv += PV_C(k);
 365                             Pvb += PVB_C(k);
 366                         }
 367                     }
 368                     Pb = stat->Pbc;
 369                     if( 2 * Pvb * Pb <= Pv ) val = 1;
 370                 }
 371                 //update FG
 372                 ((uchar*)model->foreground->imageData)[i*mask_step+j] = (uchar)(val*255);
 373                 FG_pixels_count += val;
 374             }// end if( change detection...
 375         }//for j...
 376     } //for i...
 377     //end BG/FG classification
 378
 379     //foreground segmentation
 380         //smooth FG map
 381     if( model->params.perform_morphing ){
 382         cvMorphologyEx( model->foreground, model->foreground, 0, 0, CV_MOP_OPEN, 1 );
 383         cvMorphologyEx( model->foreground, model->foreground, 0, 0, CV_MOP_CLOSE, 1 );
 384     }
 385
 386
 387     if( model->params.minArea > 0 || model->params.is_obj_without_holes ){
 388         //filter small regions
 389         cvFindContours( model->foreground, model->storage, &first_seq, sizeof(CvContour), CV_RETR_LIST );
 390         for( seq = first_seq; seq; seq = seq->h_next )
 391         {
 392             CvContour* cnt = (CvContour*)seq;
 393             if( cnt->rect.width * cnt->rect.height < model->params.minArea ||
 394                 (model->params.is_obj_without_holes && CV_IS_SEQ_HOLE(seq)) )
 395             {
 396                 //delete small contour
 397                 prev_seq = seq->h_prev;
 398                 if( prev_seq )
 399                 {
 400                     prev_seq->h_next = seq->h_next;
 401                     if( seq->h_next ) seq->h_next->h_prev = prev_seq;
 402                 }
 403                 else
 404                 {
 405                     first_seq = seq->h_next;
 406                     if( seq->h_next ) seq->h_next->h_prev = NULL;
 407                 }
 408             }
 409             else
 410             {
 411                 region_count++;
 412             }
 413         }
 414         model->foreground_regions = first_seq;
 415         cvZero(model->foreground);
 416         cvDrawContours(model->foreground, first_seq, CV_RGB(0, 0, 255), CV_RGB(0, 0, 255), 10, -1);
 417     }
 418     else{
 419         model->foreground_regions = NULL;
 420     }
 421
 422     //check ALL BG update condition
 423     if( ((float)FG_pixels_count/(model->Ftd->width*model->Ftd->height)) > CV_BGFG_FGD_BG_UPDATE_TRESH )
 424     {
 425          for( i = 0; i < model->Ftd->height; i++ )
 426              for( j = 0; j < model->Ftd->width; j++ )
 427              {
 428                  CvBGPixelStat* stat = model->pixel_stat + i * model->Ftd->width + j;
 429                  stat->is_trained_st_model = stat->is_trained_dyn_model = 1;
 430              }
 431     }
 432
 433
 434     //update BG model
 435     for( i = 0; i < model->Ftd->height; i++ )
 436     {
 437         for( j = 0; j < model->Ftd->width; j++ )
 438         {
 439             CvBGPixelStat* stat = model->pixel_stat + i * model->Ftd->width + j;
 440             CvBGPixelCStatTable* ctable = stat->ctable;
 441             CvBGPixelCCStatTable* cctable = stat->cctable;
 442
 443             uchar *curr_data = (uchar*)(curr_frame->imageData)+i*curr_frame->widthStep+j*3;
 444             uchar *prev_data = (uchar*)(prev_frame->imageData)+i*prev_frame->widthStep+j*3;
 445
 446             if( ((uchar*)model->Ftd->imageData)[i*mask_step+j] || !stat->is_trained_dyn_model )
 447             {
 448                 float alpha = stat->is_trained_dyn_model ? model->params.alpha2 : model->params.alpha3;
 449                 float diff = 0;
 450                 int dist, min_dist = 2147483647, indx = -1;
 451
 452                 //update Pb
 453                 stat->Pbcc *= (1.f-alpha);
 454                 if( !((uchar*)model->foreground->imageData)[i*mask_step+j] )
 455                 {
 456                     stat->Pbcc += alpha;
 457                 }
 458
 459                 // find best Vi match
 460                 for(k = 0; PV_CC(k) && k < model->params.N2cc; k++ )
 461                 {
 462                     // Exponential decay of memory
 463                     PV_CC(k)  *= (1-alpha);
 464                     PVB_CC(k) *= (1-alpha);
 465                     if( PV_CC(k) < MIN_PV )
 466                     {
 467                         PV_CC(k) = 0;
 468                         PVB_CC(k) = 0;
 469                         continue;
 470                     }
 471
 472                     dist = 0;
 473                     for( l = 0; l < 3; l++ )
 474                     {
 475                         int val = abs( V_CC(k,l) - prev_data[l] );
 476                         if( val > deltaCC ) break;
 477                         dist += val;
 478                         val = abs( V_CC(k,l+3) - curr_data[l] );
 479                         if( val > deltaCC) break;
 480                         dist += val;
 481                     }
 482                     if( l == 3 && dist < min_dist )
 483                     {
 484                         min_dist = dist;
 485                         indx = k;
 486                     }
 487                 }
 488
 489
 490                 if( indx < 0 )
 491                 {//N2th elem in the table is replaced by a new features
 492                     indx = model->params.N2cc - 1;
 493                     PV_CC(indx) = alpha;
 494                     PVB_CC(indx) = alpha;
 495                     //udate Vt
 496                     for( l = 0; l < 3; l++ )
 497                     {
 498                         V_CC(indx,l) = prev_data[l];
 499                         V_CC(indx,l+3) = curr_data[l];
 500                     }
 501                 }
 502                 else
 503                 {//update
 504                     PV_CC(indx) += alpha;
 505                     if( !((uchar*)model->foreground->imageData)[i*mask_step+j] )
 506                     {
 507                         PVB_CC(indx) += alpha;
 508                     }
 509                 }
 510
 511                 //re-sort CCt table by Pv
 512                 for( k = 0; k < indx; k++ )
 513                 {
 514                     if( PV_CC(k) <= PV_CC(indx) )
 515                     {
 516                         //shift elements
 517                         CvBGPixelCCStatTable tmp1, tmp2 = cctable[indx];
 518                         for( l = k; l <= indx; l++ )
 519                         {
 520                             tmp1 = cctable[l];
 521                             cctable[l] = tmp2;
 522                             tmp2 = tmp1;
 523                         }
 524                         break;
 525                     }
 526                 }
 527
 528
 529                 float sum1=0, sum2=0;
 530                 //check "once-off" changes
 531                 for(k = 0; PV_CC(k) && k < model->params.N1cc; k++ )
 532                 {
 533                     sum1 += PV_CC(k);
 534                     sum2 += PVB_CC(k);
 535                 }
 536                 if( sum1 > model->params.T ) stat->is_trained_dyn_model = 1;
 537
 538                 diff = sum1 - stat->Pbcc * sum2;
 539                 //update stat table
 540                 if( diff >  model->params.T )
 541                 {
 542                     //printf("once off change at motion mode\n");
 543                     //new BG features are discovered
 544                     for( k = 0; PV_CC(k) && k < model->params.N1cc; k++ )
 545                     {
 546                         PVB_CC(k) =
 547                             (PV_CC(k)-stat->Pbcc*PVB_CC(k))/(1-stat->Pbcc);
 548                     }
 549                     assert(stat->Pbcc<=1 && stat->Pbcc>=0);
 550                 }
 551             }
 552
 553             //case of stational pixel
 554             if( !((uchar*)model->Ftd->imageData)[i*mask_step+j] )
 555             {
 556                 float alpha = stat->is_trained_st_model ? model->params.alpha2 : model->params.alpha3;
 557                 float diff = 0;
 558                 int dist, min_dist = 2147483647, indx = -1;
 559
 560                 //update Pb
 561                 stat->Pbc *= (1.f-alpha);
 562                 if( !((uchar*)model->foreground->imageData)[i*mask_step+j] )
 563                 {
 564                     stat->Pbc += alpha;
 565                 }
 566
 567                 //find best Vi match
 568                 for( k = 0; k < model->params.N2c; k++ )
 569                 {
 570                     // Exponential decay of memory
 571                     PV_C(k) *= (1-alpha);
 572                     PVB_C(k) *= (1-alpha);
 573                     if( PV_C(k) < MIN_PV )
 574                     {
 575                         PV_C(k) = 0;
 576                         PVB_C(k) = 0;
 577                         continue;
 578                     }
 579
 580                     dist = 0;
 581                     for( l = 0; l < 3; l++ )
 582                     {
 583                         int val = abs( V_C(k,l) - curr_data[l] );
 584                         if( val > deltaC ) break;
 585                         dist += val;
 586                     }
 587                     if( l == 3 && dist < min_dist )
 588                     {
 589                         min_dist = dist;
 590                         indx = k;
 591                     }
 592                 }
 593
 594                 if( indx < 0 )
 595                 {//N2th elem in the table is replaced by a new features
 596                     indx = model->params.N2c - 1;
 597                     PV_C(indx) = alpha;
 598                     PVB_C(indx) = alpha;
 599                     //udate Vt
 600                     for( l = 0; l < 3; l++ )
 601                     {
 602                         V_C(indx,l) = curr_data[l];
 603                     }
 604                 } else
 605                 {//update
 606                     PV_C(indx) += alpha;
 607                     if( !((uchar*)model->foreground->imageData)[i*mask_step+j] )
 608                     {
 609                         PVB_C(indx) += alpha;
 610                     }
 611                 }
 612
 613                 //re-sort Ct table by Pv
 614                 for( k = 0; k < indx; k++ )
 615                 {
 616                     if( PV_C(k) <= PV_C(indx) )
 617                     {
 618                         //shift elements
 619                         CvBGPixelCStatTable tmp1, tmp2 = ctable[indx];
 620                         for( l = k; l <= indx; l++ )
 621                         {
 622                             tmp1 = ctable[l];
 623                             ctable[l] = tmp2;
 624                             tmp2 = tmp1;
 625                         }
 626                         break;
 627                     }
 628                 }
 629
 630                 //check "once-off" changes
 631                 float sum1=0, sum2=0;
 632                 for( k = 0; PV_C(k) && k < model->params.N1c; k++ )
 633                 {
 634                     sum1 += PV_C(k);
 635                     sum2 += PVB_C(k);
 636                 }
 637                 diff = sum1 - stat->Pbc * sum2;
 638                 if( sum1 > model->params.T ) stat->is_trained_st_model = 1;
 639
 640                 //update stat table
 641                 if( diff >  model->params.T )
 642                 {
 643                     //printf("once off change at stat mode\n");
 644                     //new BG features are discovered
 645                     for( k = 0; PV_C(k) && k < model->params.N1c; k++ )
 646                     {
 647                         PVB_C(k) = (PV_C(k)-stat->Pbc*PVB_C(k))/(1-stat->Pbc);
 648                     }
 649                     stat->Pbc = 1 - stat->Pbc;
 650                 }
 651             }//if !(change detection) at pixel (i,j)
 652
 653             //update the reference BG image
 654             if( !((uchar*)model->foreground->imageData)[i*mask_step+j])
 655             {
 656                 uchar* ptr = ((uchar*)model->background->imageData) + i*model->background->widthStep+j*3;
 657
 658                 if( !((uchar*)model->Ftd->imageData)[i*mask_step+j] &&
 659                     !((uchar*)model->Fbd->imageData)[i*mask_step+j] )
 660                 {
 661                     //apply IIR filter
 662                     for( l = 0; l < 3; l++ )
 663                     {
 664                         int a = cvRound(ptr[l]*(1 - model->params.alpha1) + model->params.alpha1*curr_data[l]);
 665                         ptr[l] = (uchar)a;
 666                         //((uchar*)model->background->imageData)[i*model->background->widthStep+j*3+l]*=(1 - model->params.alpha1);
 667                         //((uchar*)model->background->imageData)[i*model->background->widthStep+j*3+l] += model->params.alpha1*curr_data[l];
 668                     }
 669                 }
 670                 else
 671                 {
 672                     //background change is detected
 673                     for( l = 0; l < 3; l++ )
 674                     {
 675                         //((uchar*)model->background->imageData)[i*model->background->widthStep+j*3+l] = curr_data[l];
 676                         ptr[l] = curr_data[l];
 677                     }
 678                 }
 679             }
 680         }//j
 681     }//i
 682
 683     //keep prev frame
 684     cvCopy( curr_frame, model->prev_frame );
 685
 686     return region_count;
 687 }
 688
 689 /* End of file. */