OpenCv中混合高斯模型的实现

来源:互联网 发布:头戴耳机推荐 知乎 编辑:程序博客网 时间:2024/06/10 15:29

////////////////////////cvCreateGaussianBGModel///////////////////////////////////////////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, m, n;        // init parameters    if( parameters == NULL )          {                                params.win_size      = CV_BGFG_MOG_WINDOW_SIZE;   // 初始化阶段的帧数;用户自定义模型学 习率a=1/win_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;    //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, n = 0; i < first_frame->height; i++ ) //行    {        for( j = 0; j < first_frame->width; j++, n++ ) //列        {             const int p = i*first_frame->widthStep+j*first_frame->nChannels;           //以下几步是对第一个高斯函数做初始化            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;    //权值赋为1            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;                }            }        }    } //g_point[]:像素,g_values[]:高斯分布函数,mean[]:通道        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;} cvUpdateBGStatModel(videoFrame,bgModel);typedef int (CV_CDECL * CvUpdateBGStatModel)( IplImage* curr_frame, struct CvBGStatModel* bg_model );/////////////////////////cvUpdateBGStatModel////////////////////////////////////函数功能:背景模型的更新,不仅要更新高斯分布函数的参数,还要更新各高斯函数的权重static int CV_CDECL icvUpdateGaussianBGModel( IplImage* curr_frame, CvGaussBGModel* bg_model ){    int i, j, k, n;    int region_count = 0;    CvSeq *first_seq = NULL, *prev_seq = NULL, *seq = NULL;        bg_model->countFrames++;        for( i = 0, n = 0; i < curr_frame->height; i++ )    {        for( j = 0; j < curr_frame->width; j++, n++ )        {            int match[CV_BGFG_MOG_MAX_NGAUSSIANS];   //对高斯函数做标记,match[m]=1表示函数m为匹配的高斯分布函数            double sort_key[CV_BGFG_MOG_MAX_NGAUSSIANS]; //此数组存贮每个高斯函数的均值与方差比值            const int nChannels = curr_frame->nChannels;            const int p = curr_frame->widthStep*i+j*nChannels;                       CvGaussBGPoint* g_point = &bg_model->g_point[n];            const CvGaussBGStatModelParams bg_model_params = bg_model->params;            double pixel[4];   // pixel[]存贮当前像素的各通道RGB值            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 );            }            else            {                icvUpdatePartialWindow( pixel, nChannels, match, g_point, &bg_model_params );                if( no_match == -1)                    icvUpdatePartialNoMatch( pixel, nChannels, match, g_point, &bg_model_params );            }            icvGetSortKey( nChannels, sort_key, g_point, &bg_model_params );            icvInsertionSortGaussians( g_point, sort_key, (CvGaussBGStatModelParams *)&bg_model_params );            icvBackgroundTest( nChannels, n, i, j, match, bg_model );        }    }        //foreground filtering        //filter small regions    cvClearMemStorage(bg_model->storage);        //cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_OPEN, 1 );    //cvMorphologyEx( bg_model->foreground, bg_model->foreground, 0, 0, CV_MOP_CLOSE, 1 );        cvFindContours( bg_model->foreground, bg_model->storage, &first_seq, sizeof(CvContour), CV_RETR_LIST );    for( seq = first_seq; seq; seq = seq->h_next )    {        CvContour* cnt = (CvContour*)seq;        if( cnt->rect.width * cnt->rect.height < bg_model->params.minArea )        {            //delete small contour            prev_seq = seq->h_prev;            if( prev_seq )            {                prev_seq->h_next = seq->h_next;                if( seq->h_next ) seq->h_next->h_prev = prev_seq;            }            else            {                first_seq = seq->h_next;                if( seq->h_next ) seq->h_next->h_prev = NULL;            }        }        else        {            region_count++;        }    }    bg_model->foreground_regions = first_seq;    cvZero(bg_model->foreground);    cvDrawContours(bg_model->foreground, first_seq, CV_RGB(0, 0, 255), CV_RGB(0, 0, 255), 10, -1);        return region_count;} /////////////////////////////////////icvMatchTest////////////////////////////////////////////函数功能:拿当前像素的值与已存在的高斯分布函数比较,查找是否存在匹配的的高斯分布函数,如果有则返回 k值(高斯分布函数的序号)static int icvMatchTest( double* src_pixel, int nChannels, int* match,                         const CvGaussBGPoint* g_point,                         const CvGaussBGStatModelParams *bg_model_params ){    //参数的传递:src_pixel为piexl[]:即当前像素的各通道值    int k;    int matchPosition=-1;    for ( k = 0; k < bg_model_params->n_gauss; k++)    match[k]=0;        for ( k = 0; k < bg_model_params->n_gauss; k++)    if (g_point->g_values[k].match_sum > 0) {        double sum_d2 = 0.0;        double var_threshold = 0.0;        for(int m = 0; m < nChannels; m++)   {             double d = g_point->g_values[k].mean[m]- src_pixel[m]; //通道m的原始模型值与当前像素的值之差            sum_d2 += (d*d);            var_threshold += g_point->g_values[k].variance[m];        }               //当前sum_d2为d0,d1,d2的平方和,var_threshold的值为像素各通道方差之和        var_threshold = bg_model_params->std_threshold*                                 bg_model_params- >std_threshold*var_threshold;        if(sum_d2 < var_threshold) //查看是否可以与某高斯分布匹配 ????????????????   {            match[k] = 1;            matchPosition = k;            break; //如果和第k个高斯函数匹配,则终止与后续函数的匹配        }    }        return matchPosition;} ///////////////////////icvUpdateFullWindow//////////////////////////////////////////函数功能:更新每个高斯分布的权值(对匹配的高斯函数k加大权值,其余的则减小权值),如果前面的结果中存在匹配的高斯分布函数k,则需要再对第k个高斯分布函数的均值mean和方差variance做修正static void icvUpdateFullWindow( double* src_pixel, int nChannels, int* match,                                 CvGaussBGPoint* g_point,                                 const CvGaussBGStatModelParams *bg_model_params ){ //参数的传递:src_pixel为piexl[]:即当前帧中该像素的RGB值    const double learning_rate_weight = (1.0/(double)bg_model_params->win_size); //用户自定义模型学习率a    for(int k = 0; k < bg_model_params->n_gauss; k++){      //对每个高斯分布的权值做修正:w=(1-a)w+a*m (a:模型学习率,m是匹配,匹配就是1,不匹配就是0)        g_point->g_values[k].weight = g_point->g_values[k].weight +            (learning_rate_weight*((double)match[k] -g_point->g_values[k].weight));        if(match[k]) //如果存在匹配的高斯分布函数k(当前像素为背景像素),则需要再对第k个高斯分布函数的均值mean和方差variance更新   {            double learning_rate_gaussian = (double)match[k]/(g_point->g_values[k].weight*                (double)bg_model_params->win_size); //参数学习率p(p=a/w)            for(int m = 0; m < nChannels; m++)    {   //参数更新公式:u=(1-p)*u0+p*x; o*o=(1-p)*o*o+p*tmpDiff*tmpDiff                const double tmpDiff = src_pixel[m] - g_point->g_values[k].mean[m]; //当前像素的通道m的值与原始模型值之差                g_point->g_values[k].mean[m] = g_point->g_values[k].mean[m] + (learning_rate_gaussian * tmpDiff);                g_point->g_values[k].variance[m] = g_point->g_values[k].variance[m]+                    (learning_rate_gaussian*((tmpDiff*tmpDiff) - g_point->g_values[k].variance[m]));            }        }    }}   /////////////////////////icvUpdatePartialWindow///////////////////////////////函数功能:对所有的高斯分布函数做更新.至少每个高斯分布的权值必须修正,如果前面的结果中存在匹配的高斯分布函数k,则需要再对第k个高斯分布函数的match_sum修改,最终对那些匹配的高斯分布函数k的参数match_sum>0的做均值mean和方差variance修正static void icvUpdatePartialWindow( double* src_pixel, int nChannels, int* match, CvGaussBGPoint* g_point, const CvGaussBGStatModelParams *bg_model_params ){    int k, m;    int window_current = 0;        for( k = 0; k < bg_model_params->n_gauss; k++ )        window_current += g_point->g_values[k].match_sum;   //window_current为k个高斯分布函数的match_sum值之和        for( k = 0; k < bg_model_params->n_gauss; k++ )    {        g_point->g_values[k].match_sum += match[k]; //修正匹配的高斯分布函数k的match_sum值        double learning_rate_weight = (1.0/((double)window_current + 1.0)); //increased by one since sum         //修正每个高斯分布的权值        g_point->g_values[k].weight = g_point->g_values[k].weight +            (learning_rate_weight*((double)match[k] - g_point->g_values[k].weight));                if( g_point->g_values[k].match_sum > 0 && match[k] )        {            double learning_rate_gaussian = (double)match[k]/((double)g_point->g_values[k].match_sum);            for( m = 0; m < nChannels; m++ )            {                const double tmpDiff = src_pixel[m] - g_point->g_values[k].mean[m];                g_point->g_values[k].mean[m] = g_point->g_values[k].mean[m] +                    (learning_rate_gaussian*tmpDiff);                g_point->g_values[k].variance[m] = g_point->g_values[k].variance[m]+                    (learning_rate_gaussian*((tmpDiff*tmpDiff) - g_point->g_values[k].variance[m]));            }        }    }} //////////////////////////icvUpdateFullNoMatch////////////////////////////函数功能:当所有的高斯函数均不匹配时,说明有新的分布出现,需要将原高斯函数中sort_key最小的替换为新的高斯函数(权值小,方差大),其余的高斯函数对应的只需更新权值static void icvUpdateFullNoMatch( IplImage* gm_image, int p, int* match,                                  CvGaussBGPoint* g_point,                                  const CvGaussBGStatModelParams *bg_model_params){ //参数的传递:gm_image为当前帧curr_frame    int k, m;    double alpha;    int match_sum_total = 0;    g_point->g_values[bg_model_params->n_gauss - 1].match_sum = 1;    //将新的高斯分布函数的match_sum置为1            for( k = 0; k < bg_model_params->n_gauss ; k++ )        match_sum_total += g_point->g_values[k].match_sum;    g_point->g_values[bg_model_params->n_gauss - 1].weight = 1./(double)match_sum_total; //要给新的高斯分布函数赋一个较小的权值    //将新的高斯分布函数的variance[m]全部置为variance_init;mean[m]的值置为当前像素各通道的值    for( m = 0; m < gm_image->nChannels ; m++ )    {        g_point->g_values[bg_model_params->n_gauss - 1].variance[m] = bg_model_params->variance_init;        g_point->g_values[bg_model_params->n_gauss - 1].mean[m] = (unsigned char)gm_image->imageData[p + m];    }    //对其他的高斯分布函数做权值更新:w=(1-a)*w+a*m (a:模型学习率,m是匹配,匹配就是1,不匹配就是0)    alpha = 1.0 - (1.0/bg_model_params->win_size);   //alpha=1-a;    for( k = 0; k < bg_model_params->n_gauss - 1; k++ )    {           g_point->g_values[k].weight *= alpha;        if( match[k] )             g_point->g_values[k].weight += alpha;    }}  ////////////////////////////icvUpdatePartialNoMatch////////////////////////////////static voidicvUpdatePartialNoMatch(double *pixel,                        int nChannels,                        int* /*match*/,                        CvGaussBGPoint* g_point,                        const CvGaussBGStatModelParams *bg_model_params){    int k, m;    //new value of last one    g_point->g_values[bg_model_params->n_gauss - 1].match_sum = 1;        //get sum of all but last value of match_sum    int match_sum_total = 0;    for(k = 0; k < bg_model_params->n_gauss ; k++)        match_sum_total += g_point->g_values[k].match_sum;    for(m = 0; m < nChannels; m++)    {        //first pass mean is image value        g_point->g_values[bg_model_params->n_gauss - 1].variance[m] = bg_model_params->variance_init;        g_point->g_values[bg_model_params->n_gauss - 1].mean[m] = pixel[m];    }    for(k = 0; k < bg_model_params->n_gauss; k++)    { //更新所有高斯分布的权值        g_point->g_values[k].weight = (double)g_point->g_values[k].match_sum /            (double)match_sum_total;    }} /////////////////////////////////icvGetSortKey/////////////////////////////////////函数功能:计算各个高斯分布weight/sqrt(variance_sum)的值,后面将对该值进行排序(该值越大则表示背景的可能性就越大)static void icvGetSortKey( const int nChannels, double* sort_key, const CvGaussBGPoint* g_point,                           const CvGaussBGStatModelParams *bg_model_params ){    int k, m;    for( k = 0; k < bg_model_params->n_gauss; k++ )    {        // Avoid division by zero        if( g_point->g_values[k].match_sum > 0 )        {            // Independence assumption between components            double variance_sum = 0.0;            for( m = 0; m < nChannels; m++ )                variance_sum += g_point->g_values[k].variance[m];                        sort_key[k] = g_point->g_values[k].weight/sqrt(variance_sum);     //sort_key=w/(o*o)        }        else            sort_key[k]= 0.0;    }} //////////////////////////////icvInsertionSortGaussians////////////////////////////static void icvInsertionSortGaussians( CvGaussBGPoint* g_point, double* sort_key, CvGaussBGStatModelParams *bg_model_params ){    int i, j;    for( i = 1; i < bg_model_params->n_gauss; i++ )    {        double index = sort_key[i];        for( j = i; j > 0 && sort_key[j-1] < index; j-- )   //对sort_key[]按降序排序        {            double temp_sort_key = sort_key[j];            sort_key[j] = sort_key[j-1];            sort_key[j-1] = temp_sort_key;                        CvGaussBGValues temp_gauss_values = g_point->g_values[j];            g_point->g_values[j] = g_point->g_values[j-1];            g_point->g_values[j-1] = temp_gauss_values;        }//        sort_key[j] = index;    }} ///////////////////////////////////icvBackgroundTest/////////////////////////static void icvBackgroundTest( const int nChannels, int n, int i, int j, int *match, CvGaussBGModel* bg_model ){    int m, b;    uchar pixelValue = (uchar)255;   // 像素默认都为前景    double weight_sum = 0.0;    CvGaussBGPoint* g_point = bg_model->g_point;        for( m = 0; m < nChannels; m++)?????????????        bg_model->background->imageData[ bg_model->background->widthStep*i + j*nChannels + m] = (unsigned char)(g_point[n].g_values[0].mean[m]+0.5);         for( b = 0; b < bg_model->params.n_gauss; b++)    {        weight_sum += g_point[n].g_values[b].weight;         if( match[b] )            pixelValue = 0;   //if为真,说明该像素已与某高斯函数匹配,该像素为背景        if( weight_sum > bg_model->params.bg_threshold )               break; //如果if语句为真,则前b个高斯分布被选为描述背景的函数    }        bg_model->foreground->imageData[ bg_model->foreground->widthStep*i + j] = pixelValue;}

基于论文:An Improved Adaptive Background Mixture Model for Real-time Tracking with Shadow Detection" by P. KaewTraKulPong and R. Bowden

0 0
原创粉丝点击