cvbgfg_gaussmix.cpp

混合高斯模型

            }
            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, p, 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;
}

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 decending order
        {
            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;
    }
}


static int icvMatchTest( double* src_pixel, int nChannels, int* match,
                         const CvGaussBGPoint* g_point,
                         const CvGaussBGStatModelParams *bg_model_params )
{
    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++) {
        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];
            sum_d2 += (d*d);
            var_threshold += g_point->g_values[k].variance[m];
        }  //difference < STD_LIMIT*STD_LIMIT or difference**2 < STD_LIMIT*STD_LIMIT*VAR
        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;
        }
    }
    
    return matchPosition;
}

/*
static int icvMatchTest2( double* src_pixel, int nChannels, int* match,
                          const CvGaussBGPoint* g_point,
                          const CvGaussBGStatModelParams *bg_model_params )
{
    int k, m;
    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++ )
    {
        double sum_d2 = 0.0, var_threshold;
        for( m = 0; m < nChannels; m++ )
        {
            double d = g_point->g_values[k].mean[m]- src_pixel[m];
            sum_d2 += (d*d) / (g_point->g_values[k].variance[m] * g_point->g_values[k].variance[m]);
        }  //difference < STD_LIMIT*STD_LIMIT or difference**2 < STD_LIMIT*STD_LIMIT*VAR
        
        var_threshold = bg_model_params->std_threshold*bg_model_params->std_threshold;
        if( sum_d2 < var_threshold )
        {
            match[k] = 1;
            matchPosition = k;
            break;
        }
    }
    
    return matchPosition;
}
*/

static void icvUpdateFullWindow( double* src_pixel, int nChannels, int* match,
                                 CvGaussBGPoint* g_point,
                                 const CvGaussBGStatModelParams *bg_model_params )
{
    const double learning_rate_weight = (1.0/(double)bg_model_params->win_size);
    for(int k = 0; k < bg_model_params->n_gauss; k++){
        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]){
            double learning_rate_gaussian = (double)match[k]/(g_point->g_values[k].weight*
                (double)bg_model_params->win_size);
            for(int 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]));
            }
        }
    }
}


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;
    
    for( k = 0; k < bg_model_params->n_gauss; k++ )
    {
        g_point->g_values[k].match_sum += match[k];
        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]));
            }
        }
    }
}

static void icvUpdateFullNoMatch( IplImage* gm_image, int p, int* match,
                                  CvGaussBGPoint* g_point,
                                  const CvGaussBGStatModelParams *bg_model_params)
{
    int k, m;
    double alpha;
    int match_sum_total = 0;

    //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
    
    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;
    for( m = 0; m < gm_image->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] = (unsigned char)gm_image->imageData[p + m];
    }
    
    alpha = 1.0 - (1.0/bg_model_params->win_size);
    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;
    }
}


static void
icvUpdatePartialNoMatch(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;
    }
}

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);
        }
        else
            sort_key[k]= 0.0;
    }
}


static void icvBackgroundTest( const int nChannels, int n, int p, int *match, CvGaussBGModel* bg_model )
{
    int m, b;
    uchar pixelValue = (uchar)255; // will switch to 0 if match found
    double weight_sum = 0.0;
    CvGaussBGPoint* g_point = bg_model->g_point;
    
    for( m = 0; m < nChannels; m++)
        bg_model->background->imageData[p+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( weight_sum > bg_model->params.bg_threshold )
            break;
    }
    
    bg_model->foreground->imageData[p/nChannels] = pixelValue;
}

/* End of file. */