1 /*M///////////////////////////////////////////////////////////////////////////////////////
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
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.
10 // Intel License Agreement
12 // Copyright (C) 2000, Intel Corporation, all rights reserved.
13 // Third party copyrights are property of their respective owners.
15 // Redistribution and use in source and binary forms, with or without modification,
16 // are permitted provided that the following conditions are met:
18 // * Redistribution's of source code must retain the above copyright notice,
19 // this list of conditions and the following disclaimer.
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.
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.
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.
42 //This is based on the "An Improved Adaptive Background Mixture Model for
43 //Real-time Tracking and Shadow Detection" by P. KaewTraKulPong and R. Bowden
44 //The windowing method is used, but not the shadow detection. I make some of my
45 //own modifications which make more sense. There are some errors in some of their
47 //IplImage values of image that are useful
48 //int nSize; /* sizeof(IplImage) */
49 //int depth; /* pixel depth in bits: IPL_DEPTH_8U ...*/
50 //int nChannels; /* OpenCV functions support 1,2,3 or 4 channels */
51 //int width; /* image width in pixels */
52 //int height; /* image height in pixels */
53 //int imageSize; /* image data size in bytes in case of interleaved data)*/
54 //char *imageData; /* pointer to aligned image data */
55 //char *imageDataOrigin; /* pointer to very origin of image -deallocation */
56 //Values useful for gaussian integral
57 //0.5 - 0.19146 - 0.38292
58 //1.0 - 0.34134 - 0.68268
59 //1.5 - 0.43319 - 0.86638
60 //2.0 - 0.47725 - 0.95450
61 //2.5 - 0.49379 - 0.98758
62 //3.0 - 0.49865 - 0.99730
63 //3.5 - 0.4997674 - 0.9995348
64 //4.0 - 0.4999683 - 0.9999366
69 //internal functions for gaussian background detection
70 static void icvInsertionSortGaussians( CvGaussBGPoint* g_point, double* sort_key, CvGaussBGStatModelParams *bg_model_params );
73 Test whether pixel can be explained by background model;
74 Return -1 if no match was found; otherwise the index in match[] is returned
76 icvMatchTest(...) assumes what all color channels component exhibit the same variance
77 icvMatchTest2(...) accounts for different variances per color channel
79 static int icvMatchTest( double* src_pixel, int nChannels, int* match,
80 const CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params );
81 /*static int icvMatchTest2( double* src_pixel, int nChannels, int* match,
82 const CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params );*/
86 The update procedure differs between
87 * the initialization phase (named *Partial* ) and
88 * the normal phase (named *Full* )
89 The initalization phase is defined as not having processed <win_size> frames yet
91 static void icvUpdateFullWindow( double* src_pixel, int nChannels,
93 CvGaussBGPoint* g_point,
94 const CvGaussBGStatModelParams *bg_model_params );
95 static void icvUpdateFullNoMatch( IplImage* gm_image, int p,
97 CvGaussBGPoint* g_point,
98 const CvGaussBGStatModelParams *bg_model_params);
99 static void icvUpdatePartialWindow( double* src_pixel, int nChannels, int* match,
100 CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params );
101 static void icvUpdatePartialNoMatch( double* src_pixel, int nChannels,
103 CvGaussBGPoint* g_point,
104 const CvGaussBGStatModelParams *bg_model_params);
107 static void icvGetSortKey( const int nChannels, double* sort_key, const CvGaussBGPoint* g_point,
108 const CvGaussBGStatModelParams *bg_model_params );
109 static void icvBackgroundTest( const int nChannels, int n, int p, int *match, CvGaussBGModel* bg_model );
111 static void CV_CDECL icvReleaseGaussianBGModel( CvGaussBGModel** bg_model );
112 static int CV_CDECL icvUpdateGaussianBGModel( IplImage* curr_frame, CvGaussBGModel* bg_model );
114 //#define for if(0);else for
116 //g = 1 for first gaussian in list that matches else g = 0
117 //Rw is the learning rate for weight and Rg is leaning rate for mean and variance
118 //Ms is the match_sum which is the sum of matches for a particular gaussian
119 //Ms values are incremented until the sum of Ms values in the list equals window size L
120 //SMs is the sum of match_sums for gaussians in the list
121 //Rw = 1/SMs note the smallest Rw gets is 1/L
122 //Rg = g/Ms for SMs < L and Rg = g/(w*L) for SMs = L
123 //The list is maintained in sorted order using w/sqrt(variance) as a key
124 //If there is no match the last gaussian in the list is replaced by the new gaussian
125 //This will result in changes to SMs which results in changes in Rw and Rg.
126 //If a gaussian is replaced and SMs previously equaled L values of Ms are computed from w
127 //w[n+1] = w[n] + Rw*(g - w[n]) weight
128 //u[n+1] = u[n] + Rg*(x[n+1] - u[n]) mean value Sg is sum n values of g
129 //v[n+1] = v[n] + Rg*((x[n+1] - u[n])*(x[n+1] - u[n])) - v[n]) variance
132 CV_IMPL CvBGStatModel*
133 cvCreateGaussianBGModel( IplImage* first_frame, CvGaussBGStatModelParams* parameters )
135 CvGaussBGModel* bg_model = 0;
137 CV_FUNCNAME( "cvCreateGaussianBGModel" );
142 CvGaussBGStatModelParams params;
143 int i, j, k, n, m, p;
146 if( parameters == NULL )
148 params.win_size = CV_BGFG_MOG_WINDOW_SIZE;
149 params.bg_threshold = CV_BGFG_MOG_BACKGROUND_THRESHOLD;
150 params.std_threshold = CV_BGFG_MOG_STD_THRESHOLD;
151 params.weight_init = CV_BGFG_MOG_WEIGHT_INIT;
152 params.variance_init = CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT;
153 params.minArea = CV_BGFG_MOG_MINAREA;
154 params.n_gauss = CV_BGFG_MOG_NGAUSSIANS;
158 params = *parameters;
161 if( !CV_IS_IMAGE(first_frame) )
162 CV_ERROR( CV_StsBadArg, "Invalid or NULL first_frame parameter" );
164 CV_CALL( bg_model = (CvGaussBGModel*)cvAlloc( sizeof(*bg_model) ));
165 memset( bg_model, 0, sizeof(*bg_model) );
166 bg_model->type = CV_BG_MODEL_MOG;
167 bg_model->release = (CvReleaseBGStatModel)icvReleaseGaussianBGModel;
168 bg_model->update = (CvUpdateBGStatModel)icvUpdateGaussianBGModel;
170 bg_model->params = params;
173 CV_CALL( bg_model->g_point = (CvGaussBGPoint*)cvAlloc(sizeof(CvGaussBGPoint)*
174 ((first_frame->width*first_frame->height) + 256)));
176 CV_CALL( bg_model->background = cvCreateImage(cvSize(first_frame->width,
177 first_frame->height), IPL_DEPTH_8U, first_frame->nChannels));
178 CV_CALL( bg_model->foreground = cvCreateImage(cvSize(first_frame->width,
179 first_frame->height), IPL_DEPTH_8U, 1));
181 CV_CALL( bg_model->storage = cvCreateMemStorage());
184 var_init = 2 * params.std_threshold * params.std_threshold;
185 CV_CALL( bg_model->g_point[0].g_values =
186 (CvGaussBGValues*)cvAlloc( sizeof(CvGaussBGValues)*params.n_gauss*
187 (first_frame->width*first_frame->height + 128)));
189 for( i = 0, p = 0, n = 0; i < first_frame->height; i++ )
191 for( j = 0; j < first_frame->width; j++, n++ )
193 bg_model->g_point[n].g_values =
194 bg_model->g_point[0].g_values + n*params.n_gauss;
195 bg_model->g_point[n].g_values[0].weight = 1; //the first value seen has weight one
196 bg_model->g_point[n].g_values[0].match_sum = 1;
197 for( m = 0; m < first_frame->nChannels; m++)
199 bg_model->g_point[n].g_values[0].variance[m] = var_init;
200 bg_model->g_point[n].g_values[0].mean[m] = (unsigned char)first_frame->imageData[p + m];
202 for( k = 1; k < params.n_gauss; k++)
204 bg_model->g_point[n].g_values[k].weight = 0;
205 bg_model->g_point[n].g_values[k].match_sum = 0;
206 for( m = 0; m < first_frame->nChannels; m++){
207 bg_model->g_point[n].g_values[k].variance[m] = var_init;
208 bg_model->g_point[n].g_values[k].mean[m] = 0;
211 p += first_frame->nChannels;
215 bg_model->countFrames = 0;
219 if( cvGetErrStatus() < 0 )
221 if( bg_model && bg_model->release )
222 bg_model->release( (CvBGStatModel**)&bg_model );
227 return (CvBGStatModel*)bg_model;
232 icvReleaseGaussianBGModel( CvGaussBGModel** _bg_model )
234 CV_FUNCNAME( "icvReleaseGaussianBGModel" );
239 CV_ERROR( CV_StsNullPtr, "" );
243 CvGaussBGModel* bg_model = *_bg_model;
244 if( bg_model->g_point )
246 cvFree( &bg_model->g_point[0].g_values );
247 cvFree( &bg_model->g_point );
250 cvReleaseImage( &bg_model->background );
251 cvReleaseImage( &bg_model->foreground );
252 cvReleaseMemStorage(&bg_model->storage);
253 memset( bg_model, 0, sizeof(*bg_model) );
262 icvUpdateGaussianBGModel( IplImage* curr_frame, CvGaussBGModel* bg_model )
265 int region_count = 0;
266 CvSeq *first_seq = NULL, *prev_seq = NULL, *seq = NULL;
268 bg_model->countFrames++;
270 for( i = 0; i < curr_frame->height; i++ )
272 for( j = 0; j < curr_frame->width; j++ )
274 int match[CV_BGFG_MOG_MAX_NGAUSSIANS];
275 double sort_key[CV_BGFG_MOG_MAX_NGAUSSIANS];
276 const int nChannels = curr_frame->nChannels;
277 const int n = i*curr_frame->width+j;
278 const int p = n*curr_frame->nChannels;
281 CvGaussBGPoint* g_point = &bg_model->g_point[n];
282 const CvGaussBGStatModelParams bg_model_params = bg_model->params;
286 for( k = 0; k < nChannels; k++ )
287 pixel[k] = (uchar)curr_frame->imageData[p+k];
289 no_match = icvMatchTest( pixel, nChannels, match, g_point, &bg_model_params );
290 if( bg_model->countFrames == bg_model->params.win_size )
292 icvUpdateFullWindow( pixel, nChannels, match, g_point, &bg_model->params );
294 icvUpdateFullNoMatch( curr_frame, p, match, g_point, &bg_model_params );
298 icvUpdatePartialWindow( pixel, nChannels, match, g_point, &bg_model_params );
300 icvUpdatePartialNoMatch( pixel, nChannels, match, g_point, &bg_model_params );
302 icvGetSortKey( nChannels, sort_key, g_point, &bg_model_params );
303 icvInsertionSortGaussians( g_point, sort_key, (CvGaussBGStatModelParams *)&bg_model_params );
304 icvBackgroundTest( nChannels, n, p, match, bg_model );
308 //foreground filtering
310 //filter small regions
311 cvClearMemStorage(bg_model->storage);
313 //cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_OPEN, 1 );
314 //cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_CLOSE, 1 );
316 cvFindContours( bg_model->foreground, bg_model->storage, &first_seq, sizeof(CvContour), CV_RETR_LIST );
317 for( seq = first_seq; seq; seq = seq->h_next )
319 CvContour* cnt = (CvContour*)seq;
320 if( cnt->rect.width * cnt->rect.height < bg_model->params.minArea )
322 //delete small contour
323 prev_seq = seq->h_prev;
326 prev_seq->h_next = seq->h_next;
327 if( seq->h_next ) seq->h_next->h_prev = prev_seq;
331 first_seq = seq->h_next;
332 if( seq->h_next ) seq->h_next->h_prev = NULL;
340 bg_model->foreground_regions = first_seq;
341 cvZero(bg_model->foreground);
342 cvDrawContours(bg_model->foreground, first_seq, CV_RGB(0, 0, 255), CV_RGB(0, 0, 255), 10, -1);
347 static void icvInsertionSortGaussians( CvGaussBGPoint* g_point, double* sort_key, CvGaussBGStatModelParams *bg_model_params )
350 for( i = 1; i < bg_model_params->n_gauss; i++ )
352 double index = sort_key[i];
353 for( j = i; j > 0 && sort_key[j-1] < index; j-- ) //sort decending order
355 double temp_sort_key = sort_key[j];
356 sort_key[j] = sort_key[j-1];
357 sort_key[j-1] = temp_sort_key;
359 CvGaussBGValues temp_gauss_values = g_point->g_values[j];
360 g_point->g_values[j] = g_point->g_values[j-1];
361 g_point->g_values[j-1] = temp_gauss_values;
363 // sort_key[j] = index;
368 static int icvMatchTest( double* src_pixel, int nChannels, int* match,
369 const CvGaussBGPoint* g_point,
370 const CvGaussBGStatModelParams *bg_model_params )
373 int matchPosition=-1;
374 for ( k = 0; k < bg_model_params->n_gauss; k++) match[k]=0;
376 for ( k = 0; k < bg_model_params->n_gauss; k++) {
378 double var_threshold = 0.0;
379 for(int m = 0; m < nChannels; m++){
380 double d = g_point->g_values[k].mean[m]- src_pixel[m];
382 var_threshold += g_point->g_values[k].variance[m];
383 } //difference < STD_LIMIT*STD_LIMIT or difference**2 < STD_LIMIT*STD_LIMIT*VAR
384 var_threshold = bg_model_params->std_threshold*bg_model_params->std_threshold*var_threshold;
385 if(sum_d2 < var_threshold){
392 return matchPosition;
396 static int icvMatchTest2( double* src_pixel, int nChannels, int* match,
397 const CvGaussBGPoint* g_point,
398 const CvGaussBGStatModelParams *bg_model_params )
401 int matchPosition=-1;
403 for( k = 0; k < bg_model_params->n_gauss; k++ )
406 for( k = 0; k < bg_model_params->n_gauss; k++ )
408 double sum_d2 = 0.0, var_threshold;
409 for( m = 0; m < nChannels; m++ )
411 double d = g_point->g_values[k].mean[m]- src_pixel[m];
412 sum_d2 += (d*d) / (g_point->g_values[k].variance[m] * g_point->g_values[k].variance[m]);
413 } //difference < STD_LIMIT*STD_LIMIT or difference**2 < STD_LIMIT*STD_LIMIT*VAR
415 var_threshold = bg_model_params->std_threshold*bg_model_params->std_threshold;
416 if( sum_d2 < var_threshold )
424 return matchPosition;
428 static void icvUpdateFullWindow( double* src_pixel, int nChannels, int* match,
429 CvGaussBGPoint* g_point,
430 const CvGaussBGStatModelParams *bg_model_params )
432 const double learning_rate_weight = (1.0/(double)bg_model_params->win_size);
433 for(int k = 0; k < bg_model_params->n_gauss; k++){
434 g_point->g_values[k].weight = g_point->g_values[k].weight +
435 (learning_rate_weight*((double)match[k] -
436 g_point->g_values[k].weight));
438 double learning_rate_gaussian = (double)match[k]/(g_point->g_values[k].weight*
439 (double)bg_model_params->win_size);
440 for(int m = 0; m < nChannels; m++){
441 const double tmpDiff = src_pixel[m] - g_point->g_values[k].mean[m];
442 g_point->g_values[k].mean[m] = g_point->g_values[k].mean[m] +
443 (learning_rate_gaussian * tmpDiff);
444 g_point->g_values[k].variance[m] = g_point->g_values[k].variance[m]+
445 (learning_rate_gaussian*((tmpDiff*tmpDiff) - g_point->g_values[k].variance[m]));
452 static void icvUpdatePartialWindow( double* src_pixel, int nChannels, int* match, CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params )
455 int window_current = 0;
457 for( k = 0; k < bg_model_params->n_gauss; k++ )
458 window_current += g_point->g_values[k].match_sum;
460 for( k = 0; k < bg_model_params->n_gauss; k++ )
462 g_point->g_values[k].match_sum += match[k];
463 double learning_rate_weight = (1.0/((double)window_current + 1.0)); //increased by one since sum
464 g_point->g_values[k].weight = g_point->g_values[k].weight +
465 (learning_rate_weight*((double)match[k] - g_point->g_values[k].weight));
467 if( g_point->g_values[k].match_sum > 0 && match[k] )
469 double learning_rate_gaussian = (double)match[k]/((double)g_point->g_values[k].match_sum);
470 for( m = 0; m < nChannels; m++ )
472 const double tmpDiff = src_pixel[m] - g_point->g_values[k].mean[m];
473 g_point->g_values[k].mean[m] = g_point->g_values[k].mean[m] +
474 (learning_rate_gaussian*tmpDiff);
475 g_point->g_values[k].variance[m] = g_point->g_values[k].variance[m]+
476 (learning_rate_gaussian*((tmpDiff*tmpDiff) - g_point->g_values[k].variance[m]));
482 static void icvUpdateFullNoMatch( IplImage* gm_image, int p, int* match,
483 CvGaussBGPoint* g_point,
484 const CvGaussBGStatModelParams *bg_model_params)
488 int match_sum_total = 0;
490 //new value of last one
491 g_point->g_values[bg_model_params->n_gauss - 1].match_sum = 1;
493 //get sum of all but last value of match_sum
495 for( k = 0; k < bg_model_params->n_gauss ; k++ )
496 match_sum_total += g_point->g_values[k].match_sum;
498 g_point->g_values[bg_model_params->n_gauss - 1].weight = 1./(double)match_sum_total;
499 for( m = 0; m < gm_image->nChannels ; m++ )
501 // first pass mean is image value
502 g_point->g_values[bg_model_params->n_gauss - 1].variance[m] = bg_model_params->variance_init;
503 g_point->g_values[bg_model_params->n_gauss - 1].mean[m] = (unsigned char)gm_image->imageData[p + m];
506 alpha = 1.0 - (1.0/bg_model_params->win_size);
507 for( k = 0; k < bg_model_params->n_gauss - 1; k++ )
509 g_point->g_values[k].weight *= alpha;
511 g_point->g_values[k].weight += alpha;
517 icvUpdatePartialNoMatch(double *pixel,
520 CvGaussBGPoint* g_point,
521 const CvGaussBGStatModelParams *bg_model_params)
524 //new value of last one
525 g_point->g_values[bg_model_params->n_gauss - 1].match_sum = 1;
527 //get sum of all but last value of match_sum
528 int match_sum_total = 0;
529 for(k = 0; k < bg_model_params->n_gauss ; k++)
530 match_sum_total += g_point->g_values[k].match_sum;
532 for(m = 0; m < nChannels; m++)
534 //first pass mean is image value
535 g_point->g_values[bg_model_params->n_gauss - 1].variance[m] = bg_model_params->variance_init;
536 g_point->g_values[bg_model_params->n_gauss - 1].mean[m] = pixel[m];
538 for(k = 0; k < bg_model_params->n_gauss; k++)
540 g_point->g_values[k].weight = (double)g_point->g_values[k].match_sum /
541 (double)match_sum_total;
545 static void icvGetSortKey( const int nChannels, double* sort_key, const CvGaussBGPoint* g_point,
546 const CvGaussBGStatModelParams *bg_model_params )
549 for( k = 0; k < bg_model_params->n_gauss; k++ )
551 // Avoid division by zero
552 if( g_point->g_values[k].match_sum > 0 )
554 // Independence assumption between components
555 double variance_sum = 0.0;
556 for( m = 0; m < nChannels; m++ )
557 variance_sum += g_point->g_values[k].variance[m];
559 sort_key[k] = g_point->g_values[k].weight/sqrt(variance_sum);
567 static void icvBackgroundTest( const int nChannels, int n, int p, int *match, CvGaussBGModel* bg_model )
570 uchar pixelValue = (uchar)255; // will switch to 0 if match found
571 double weight_sum = 0.0;
572 CvGaussBGPoint* g_point = bg_model->g_point;
574 for( m = 0; m < nChannels; m++)
575 bg_model->background->imageData[p+m] = (unsigned char)(g_point[n].g_values[0].mean[m]+0.5);
577 for( b = 0; b < bg_model->params.n_gauss; b++)
579 weight_sum += g_point[n].g_values[b].weight;
582 if( weight_sum > bg_model->params.bg_threshold )
586 bg_model->foreground->imageData[p/nChannels] = pixelValue;