cvbgfg_gaussmix.cpp

混合高斯模型

/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                        Intel License Agreement
//
// Copyright (C) 2000, Intel Corporation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of Intel Corporation may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/


//This is based on the "An Improved Adaptive Background Mixture Model for
//Real-time Tracking and Shadow Detection" by P. KaewTraKulPong and R. Bowden
//The windowing method is used, but not the shadow detection. I make some of my
//own modifications which make more sense. There are some errors in some of their
//equations.
//IplImage values of image that are useful
//int  nSize;         /* sizeof(IplImage) */
//int  depth;         /* pixel depth in bits: IPL_DEPTH_8U ...*/
//int  nChannels;     /* OpenCV functions support 1,2,3 or 4 channels */
//int  width;         /* image width in pixels */
//int  height;        /* image height in pixels */
//int  imageSize;     /* image data size in bytes in case of interleaved data)*/
//char *imageData;    /* pointer to aligned image data */
//char *imageDataOrigin; /* pointer to very origin of image -deallocation */
//Values useful for gaussian integral
//0.5 - 0.19146 - 0.38292
//1.0 - 0.34134 - 0.68268
//1.5 - 0.43319 - 0.86638
//2.0 - 0.47725 - 0.95450
//2.5 - 0.49379 - 0.98758
//3.0 - 0.49865 - 0.99730
//3.5 - 0.4997674 - 0.9995348
//4.0 - 0.4999683 - 0.9999366

#include "_cvaux.h"


//internal functions for gaussian background detection
static void icvInsertionSortGaussians( CvGaussBGPoint* g_point, double* sort_key, CvGaussBGStatModelParams *bg_model_params );

/* 
   Test whether pixel can be explained by background model; 
   Return -1 if no match was found; otherwise the index in match[] is returned

   icvMatchTest(...) assumes what all color channels component exhibit the same variance
   icvMatchTest2(...) accounts for different variances per color channel
 */
static int icvMatchTest( double* src_pixel, int nChannels, int* match, 
                 const CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params );
/*static int icvMatchTest2( double* src_pixel, int nChannels, int* match, 
                 const CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params );*/


/* 
   The update procedure differs between  
      * the initialization phase (named *Partial* ) and
      * the normal phase (named *Full* )
   The initalization phase is defined as not having processed <win_size> frames yet
 */
static void icvUpdateFullWindow( double* src_pixel, int nChannels, 
                         int* match,
                         CvGaussBGPoint* g_point, 
                         const CvGaussBGStatModelParams *bg_model_params );
static void icvUpdateFullNoMatch( IplImage* gm_image, int p, 
                          int* match, 
                          CvGaussBGPoint* g_point, 
                          const CvGaussBGStatModelParams *bg_model_params);
static void icvUpdatePartialWindow( double* src_pixel, int nChannels, int* match, 
                            CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params );
static void icvUpdatePartialNoMatch( double* src_pixel, int nChannels, 
                             int* match, 
                             CvGaussBGPoint* g_point, 
                             const CvGaussBGStatModelParams *bg_model_params);


static void icvGetSortKey( const int nChannels, double* sort_key, const CvGaussBGPoint* g_point, 
                    const CvGaussBGStatModelParams *bg_model_params );
static void icvBackgroundTest( const int nChannels, int n, int p, int *match, CvGaussBGModel* bg_model );

static void CV_CDECL icvReleaseGaussianBGModel( CvGaussBGModel** bg_model );
static int CV_CDECL icvUpdateGaussianBGModel( IplImage* curr_frame, CvGaussBGModel*  bg_model );

//#define for if(0);else for

//g = 1 for first gaussian in list that matches else g = 0
//Rw is the learning rate for weight and Rg is leaning rate for mean and variance
//Ms is the match_sum which is the sum of matches for a particular gaussian
//Ms values are incremented until the sum of Ms values in the list equals window size L
//SMs is the sum of match_sums for gaussians in the list
//Rw = 1/SMs note the smallest Rw gets is 1/L
//Rg = g/Ms for SMs < L and Rg = g/(w*L) for SMs = L
//The list is maintained in sorted order using w/sqrt(variance) as a key
//If there is no match the last gaussian in the list is replaced by the new gaussian
//This will result in changes to SMs which results in changes in Rw and Rg.
//If a gaussian is replaced and SMs previously equaled L values of Ms are computed from w
//w[n+1] = w[n] + Rw*(g - w[n])   weight
//u[n+1] = u[n] + Rg*(x[n+1] - u[n]) mean value Sg is sum n values of g
//v[n+1] = v[n] + Rg*((x[n+1] - u[n])*(x[n+1] - u[n])) - v[n]) variance
//

CV_IMPL CvBGStatModel*
cvCreateGaussianBGModel( IplImage* first_frame, CvGaussBGStatModelParams* parameters )
{
    CvGaussBGModel* bg_model = 0;
    
    CV_FUNCNAME( "cvCreateGaussianBGModel" );
    
    __BEGIN__;
    
    double var_init;
    CvGaussBGStatModelParams params;
    int i, j, k, n, m, p;
    
    //init parameters
    if( parameters == NULL )
    {
        params.win_size = CV_BGFG_MOG_WINDOW_SIZE;
        params.bg_threshold = CV_BGFG_MOG_BACKGROUND_THRESHOLD;
        params.std_threshold = CV_BGFG_MOG_STD_THRESHOLD;
        params.weight_init = CV_BGFG_MOG_WEIGHT_INIT;
        params.variance_init = CV_BGFG_MOG_SIGMA_INIT*CV_BGFG_MOG_SIGMA_INIT;
        params.minArea = CV_BGFG_MOG_MINAREA;
        params.n_gauss = CV_BGFG_MOG_NGAUSSIANS;
    }
    else
    {
        params = *parameters;
    }
    
    if( !CV_IS_IMAGE(first_frame) )
        CV_ERROR( CV_StsBadArg, "Invalid or NULL first_frame parameter" );
    
    CV_CALL( bg_model = (CvGaussBGModel*)cvAlloc( sizeof(*bg_model) ));
    memset( bg_model, 0, sizeof(*bg_model) );
    bg_model->type = CV_BG_MODEL_MOG;
    bg_model->release = (CvReleaseBGStatModel)icvReleaseGaussianBGModel;
    bg_model->update = (CvUpdateBGStatModel)icvUpdateGaussianBGModel;
    
    bg_model->params = params;
    
    //prepare storages
    CV_CALL( bg_model->g_point = (CvGaussBGPoint*)cvAlloc(sizeof(CvGaussBGPoint)*
        ((first_frame->width*first_frame->height) + 256)));
    
    CV_CALL( bg_model->background = cvCreateImage(cvSize(first_frame->width,
        first_frame->height), IPL_DEPTH_8U, first_frame->nChannels));
    CV_CALL( bg_model->foreground = cvCreateImage(cvSize(first_frame->width,
        first_frame->height), IPL_DEPTH_8U, 1));
    
    CV_CALL( bg_model->storage = cvCreateMemStorage());
    
    //initializing
    var_init = 2 * params.std_threshold * params.std_threshold;
    CV_CALL( bg_model->g_point[0].g_values =
        (CvGaussBGValues*)cvAlloc( sizeof(CvGaussBGValues)*params.n_gauss*
        (first_frame->width*first_frame->height + 128)));
    
    for( i = 0, p = 0, n = 0; i < first_frame->height; i++ )
    {
        for( j = 0; j < first_frame->width; j++, n++ )
        {
            bg_model->g_point[n].g_values =
                bg_model->g_point[0].g_values + n*params.n_gauss;
            bg_model->g_point[n].g_values[0].weight = 1;    //the first value seen has weight one
            bg_model->g_point[n].g_values[0].match_sum = 1;
            for( m = 0; m < first_frame->nChannels; m++)
            {
                bg_model->g_point[n].g_values[0].variance[m] = var_init;
                bg_model->g_point[n].g_values[0].mean[m] = (unsigned char)first_frame->imageData[p + m];
            }
            for( k = 1; k < params.n_gauss; k++)
            {
                bg_model->g_point[n].g_values[k].weight = 0;
                bg_model->g_point[n].g_values[k].match_sum = 0;
                for( m = 0; m < first_frame->nChannels; m++){
                    bg_model->g_point[n].g_values[k].variance[m] = var_init;
                    bg_model->g_point[n].g_values[k].mean[m] = 0;
                }
            }
            p += first_frame->nChannels;
        }
    }
    
    bg_model->countFrames = 0;
    
    __END__;
    
    if( cvGetErrStatus() < 0 )
    {
        CvBGStatModel* base_ptr = (CvBGStatModel*)bg_model;
        
        if( bg_model && bg_model->release )
            bg_model->release( &base_ptr );
        else
            cvFree( &bg_model );
        bg_model = 0;
    }
    
    return (CvBGStatModel*)bg_model;
}


static void CV_CDECL
icvReleaseGaussianBGModel( CvGaussBGModel** _bg_model )
{
    CV_FUNCNAME( "icvReleaseGaussianBGModel" );

    __BEGIN__;
    
    if( !_bg_model )
        CV_ERROR( CV_StsNullPtr, "" );

    if( *_bg_model )
    {
        CvGaussBGModel* bg_model = *_bg_model;
        if( bg_model->g_point )
        {
            cvFree( &bg_model->g_point[0].g_values );
            cvFree( &bg_model->g_point );
        }
        
        cvReleaseImage( &bg_model->background );
        cvReleaseImage( &bg_model->foreground );
        cvReleaseMemStorage(&bg_model->storage);
        memset( bg_model, 0, sizeof(*bg_model) );
        cvFree( _bg_model );
    }

    __END__;
}


static int CV_CDECL
icvUpdateGaussianBGModel( IplImage* curr_frame, CvGaussBGModel*  bg_model )
{
    int i, j, k;
    int region_count = 0;
    CvSeq *first_seq = NULL, *prev_seq = NULL, *seq = NULL;
    
    bg_model->countFrames++;
    
    for( i = 0; i < curr_frame->height; i++ )
    {
        for( j = 0; j < curr_frame->width; j++ )
        {
            int match[CV_BGFG_MOG_MAX_NGAUSSIANS];
            double sort_key[CV_BGFG_MOG_MAX_NGAUSSIANS];
            const int nChannels = curr_frame->nChannels;
            const int n = i*curr_frame->width+j;
            const int p = n*curr_frame->nChannels;
            
            // A few short cuts
            CvGaussBGPoint* g_point = &bg_model->g_point[n];
            const CvGaussBGStatModelParams bg_model_params = bg_model->params;
            double pixel[4];
            int no_match;
            
            for( k = 0; k < nChannels; k++ )
                pixel[k] = (uchar)curr_frame->imageData[p+k];
            
            no_match = icvMatchTest( pixel, nChannels, match, g_point, &bg_model_params );
            if( bg_model->countFrames == bg_model->params.win_size )
            {
                icvUpdateFullWindow( pixel, nChannels, match, g_point, &bg_model->params );
                if( no_match == -1)
                    icvUpdateFullNoMatch( curr_frame, p, match, g_point, &bg_model_params );