cvhmm.cpp.svn-base

来自「非结构化路识别」· SVN-BASE 代码 · 共 1,633 行 · 第 1/4 页

    for( last_obs = m_minNumObs-1, i = 0; last_obs < m_maxNumObs; last_obs++, i++ )
    {
        int t;

        /******************************************************************/
        /*    Viterbi termination                                         */
    
        if ( m_ProbType == _CV_LAST_STATE )
        {
            m_MaxGamma[i] = m_Gamma[last_obs * num_states + num_states - 1];
            q[i][last_obs] = num_states - 1;
        }
        else if( m_ProbType == _CV_BEST_STATE )
        {
            int k;
            q[i][last_obs] = 0;  
            m_MaxGamma[i] = m_Gamma[last_obs * num_states + 0]; 
        
            for(k = 1; k < num_states; k++)
            { 
                if ( m_Gamma[last_obs * num_states + k] > m_MaxGamma[i] )
                {
                    m_MaxGamma[i] = m_Gamma[last_obs * num_states + k];
                    q[i][last_obs] = k;
                }    
            }
        } 
    
        /******************************************************************/
        /*    Viterbi backtracking                                        */
        for  (t = last_obs-1; t >= 0; t--)
        {
            q[i][t] = m_csi[(t+1) * num_states + q[i][t+1] ];   
        } 
    }            
    
    /* memory free */
    icvFree( &m_pi );
    icvFree( &m_csi );
    icvDeleteMatrix( m_Gamma );   
       
    return CV_NO_ERR;
} 

/*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: icvEViterbi
//    Purpose: The function calculates the embedded Viterbi algorithm
//             for 1 image 
//    Context:
//    Parameters:  
//             obs_info - observations
//             hmm      - HMM
//                
//    Returns: the Embedded Viterbi probability (float) 
//             and do state segmentation of observations
//
//    Notes: 
//F*/
IPCVAPI_IMPL(float, icvEViterbi, (CvImgObsInfo* obs_info, CvEHMM* hmm))
{
    int    i, j, counter;
    float  log_likelihood;

    float inv_obs_x = 1.f / obs_info->obs_x;

    CvEHMMState* first_state = hmm->u.ehmm->u.state;
    
    /* memory allocation for superB */
    CvMatr32f superB = icvCreateMatrix_32f(hmm->num_states, obs_info->obs_y );
    
    /* memory allocation for q */
    int*** q = (int***)icvAlloc( hmm->num_states * sizeof(int**) );
    int* super_q = (int*)icvAlloc( obs_info->obs_y * sizeof(int) );
    
    for (i = 0; i < hmm->num_states; i++)
    {
        q[i] = (int**)icvAlloc( obs_info->obs_y * sizeof(int*) );
        
        for (j = 0; j < obs_info->obs_y ; j++)
        {
            q[i][j] = (int*)icvAlloc( obs_info->obs_x * sizeof(int) );
        }
    }             
    
    /* start Viterbi segmentation */
    for (i = 0; i < hmm->num_states; i++)
    {
        CvEHMM* ehmm = &(hmm->u.ehmm[i]);
        
        for (j = 0; j < obs_info->obs_y; j++)
        {
            float max_gamma;
            
            /* 1D HMM Viterbi segmentation */
            icvViterbiSegmentation( ehmm->num_states, obs_info->obs_x, 
                ehmm->transP, ehmm->obsProb[j], 0, 
                _CV_LAST_STATE, &q[i][j], obs_info->obs_x,
                obs_info->obs_x, &max_gamma);
            
            superB[j * hmm->num_states + i] = max_gamma * inv_obs_x;
        }
    }
    
    /* perform global Viterbi segmentation (i.e. process higher-level HMM) */
    
    icvViterbiSegmentation( hmm->num_states, obs_info->obs_y, 
                             hmm->transP, superB, 0, 
                             _CV_LAST_STATE, &super_q, obs_info->obs_y,
                             obs_info->obs_y, &log_likelihood );
    
    log_likelihood /= obs_info->obs_y ;   
    
    
    counter = 0;
    /* assign new state to observation vectors */
    for (i = 0; i < obs_info->obs_y; i++)
    {   
        for (j = 0; j < obs_info->obs_x; j++, counter++)
        {
            int superstate = super_q[i];
            int state = hmm->u.ehmm[superstate].u.state - first_state;
            
            obs_info->state[2 * counter] = superstate;
            obs_info->state[2 * counter + 1] = state + q[superstate][i][j];
        }
    }
    
    /* memory deallocation for superB */
    icvDeleteMatrix( superB );
    
    /*memory deallocation for q */
    for (i = 0; i < hmm->num_states; i++)
    {
        for (j = 0; j < obs_info->obs_y ; j++)
        {
            icvFree( &q[i][j] );
        }
        icvFree( &q[i] );
    }
    
    icvFree( &q );
    icvFree( &super_q );
    
    return log_likelihood;
}  

IPCVAPI_IMPL(CvStatus, icvEstimateHMMStateParams, (CvImgObsInfo** obs_info_array,
                                                   int num_img, CvEHMM* hmm))

{
    /* compute gamma, weights, means, vars */
    int k, i, j, m;
    int total = 0;
    int vect_len = obs_info_array[0]->obs_size;

    float start_log_var_val = LN2PI * vect_len;

    CvVect32f tmp_vect = icvCreateVector_32f( vect_len );
    
    CvEHMMState* first_state = hmm->u.ehmm[0].u.state;

    assert( sizeof(float) == sizeof(int) );

    for(i = 0; i < hmm->num_states; i++ )
    {
        total+= hmm->u.ehmm[i].num_states;
    }

    /***************Gamma***********************/
    /* initialize gamma */
    for( i = 0; i < total; i++ )
    {
        for (m = 0; m < first_state[i].num_mix; m++)
        {
            ((int*)(first_state[i].weight))[m] = 0;
        }     
    }
    
    /* maybe gamma must be computed in mixsegm process ?? */

    /* compute gamma */
    for (k = 0; k < num_img; k++)
    {
        CvImgObsInfo* info = obs_info_array[k];
        int num_obs = info->obs_y * info->obs_x;

        for (i = 0; i < num_obs; i++)
        {
            int state, mixture;
            state = info->state[2*i + 1];
            mixture = info->mix[i];
            /* computes gamma - number of observations corresponding 
               to every mixture of every state */ 
            ((int*)(first_state[state].weight))[mixture] += 1;              
        }
    }     
    /***************Mean and Var***********************/
    /* compute means and variances of every item */
    /* initially variance placed to inv_var */
    /* zero mean and variance */
    for (i = 0; i < total; i++)
    {
        memset( (void*)first_state[i].mu, 0, first_state[i].num_mix * vect_len * 
                                                                         sizeof(float) );
        memset( (void*)first_state[i].inv_var, 0, first_state[i].num_mix * vect_len * 
                                                                         sizeof(float) );
    }
    
    /* compute sums */
    for (i = 0; i < num_img; i++)
    {
        CvImgObsInfo* info = obs_info_array[i];
        int total_obs = info->obs_x * info->obs_y;

        float* vector = info->obs;

        for (j = 0; j < total_obs; j++, vector+=vect_len )
        {   
            int state = info->state[2 * j + 1];
            int mixture = info->mix[j]; 
            
            CvVect32f mean  = first_state[state].mu + mixture * vect_len;
            CvVect32f mean2 = first_state[state].inv_var + mixture * vect_len;
            
            icvAddVector_32f( mean, vector, mean, vect_len );
            icvAddSquare_32f_C1IR( vector, vect_len * sizeof(float),
                                    mean2, vect_len * sizeof(float), cvSize(vect_len, 1) ); 
        }   
    }
    
    /*compute the means and variances */
    /* assume gamma already computed */
    for (i = 0; i < total; i++)
    {           
        CvEHMMState* state = &(first_state[i]);

        for (m = 0; m < state->num_mix; m++)
        {
            int k;
            CvVect32f mu  = state->mu + m * vect_len;
            CvVect32f invar = state->inv_var + m * vect_len;             
            
            if ( ((int*)state->weight)[m] > 1)
            {
                float inv_gamma = 1.f/((int*)(state->weight))[m];
            
                icvScaleVector_32f( mu, mu, vect_len, inv_gamma);
                icvScaleVector_32f( invar, invar, vect_len, inv_gamma);
            }

            icvMulVectors_32f(mu, mu, tmp_vect, vect_len);
            icvSubVector_32f( invar, tmp_vect, invar, vect_len);     
            
            /* low bound of variance - 100 (Ara's experimental result) */
            for( k = 0; k < vect_len; k++ )
            {
                invar[k] = (invar[k] > 100.f) ? invar[k] : 100.f;
            }

            /* compute log_var */
            state->log_var_val[m] = start_log_var_val;
            for( k = 0; k < vect_len; k++ )
            {
                state->log_var_val[m] += (float)log( invar[k] );
            }           

            /* SMOLI 27.10.2000 */
            state->log_var_val[m] *= 0.5;


            /* compute inv_var = 1/sqrt(2*variance) */
            icvScaleVector_32f(invar, invar, vect_len, 2.f );
            icvbInvSqrt_32f(invar, invar, vect_len );
        }
    }
  
    /***************Weights***********************/
    /* normilize gammas - i.e. compute mixture weights */
    
    //compute weights
    for (i = 0; i < total; i++)
    {           
        int gamma_total = 0;
        float norm;

        for (m = 0; m < first_state[i].num_mix; m++)
        {
            gamma_total += ((int*)(first_state[i].weight))[m];  
        }

        norm = gamma_total ? (1.f/(float)gamma_total) : 0.f;
            
        for (m = 0; m < first_state[i].num_mix; m++)
        {
            first_state[i].weight[m] = ((int*)(first_state[i].weight))[m] * norm;
        } 
    }                                               

    icvDeleteVector( tmp_vect);
    return CV_NO_ERR; 
}   

/*
CvStatus icvLightingCorrection8uC1R( uchar* img, CvSize roi, int src_step )
{
    int i, j;
    int width = roi.width;
    int height = roi.height;
    
    float x1, x2, y1, y2;
    int f[3] = {0, 0, 0};
    float a[3] = {0, 0, 0};
    
    float h1;
    float h2;
    
    float c1,c2;
    
    float min = FLT_MAX;
    float max = -FLT_MAX;
    float correction;
    
    float* float_img = icvAlloc( width * height * sizeof(float) );
    
    x1 = width * (width + 1) / 2.0f; // Sum (1, ... , width)
    x2 = width * (width + 1 ) * (2 * width + 1) / 6.0f; // Sum (1^2, ... , width^2)
    y1 = height * (height + 1)/2.0f; // Sum (1, ... , width)
    y2 = height * (height + 1 ) * (2 * height + 1) / 6.0f; // Sum (1^2, ... , width^2)
    
    
    // extract grayvalues
    for (i = 0; i < height; i++)
    {
        for (j = 0; j < width; j++)
        {
            f[2] = f[2] + j * img[i*src_step + j];
            f[1] = f[1] + i * img[i*src_step + j];
            f[0] = f[0] +     img[i*src_step + j];
        }
    }
    
    h1 = (float)f[0] * (float)x1 / (float)width;
    h2 = (float)f[0] * (float)y1 / (float)height;
    
    a[2] = ((float)f[2] - h1) / (float)(x2*height - x1*x1*height/(float)width);
    a[1] = ((float)f[1] - h2) / (float)(y2*width - y1*y1*width/(float)height);
    a[0] = (float)f[0]/(float)(width*height) - (float)y1*a[1]/(float)height - 
        (float)x1*a[2]/(float)width;
    
    for (i = 0; i < height; i++) 
    {    
        for (j = 0; j < width; j++)
        {
            
            correction = a[0] + a[1]*(float)i + a[2]*(float)j;
            
            float_img[i*width + j] = img[i*src_step + j] - correction;
            
            if (float_img[i*width + j] < min) min = float_img[i*width+j];
            if (float_img[i*width + j] > max) max = float_img[i*width+j];
        }
    }
    
    //rescaling to the range 0:255
    c2 = 0;
    if (max == min)
        c2 = 255.0f;
    else
        c2 = 255.0f/(float)(max - min);
    
    c1 = (-(float)min)*c2;
    
    for (i = 0; i < height; i++)
    {
        for (j = 0; j < width; j++)
        {
            int value = (int)floor(c2*float_img[i*width + j] + c1);
            if (value < 0) value = 0;
            if (value > 255) value = 255;
            img[i*src_step + j] = (uchar)value;
        }
    }

    icvFree( &float_img );
    return CV_NO_ERR;
}
                              

CvStatus icvLightingCorrection( icvImage* img ) 
{
    CvSize roi;
    if ( img->type != IPL_DEPTH_8U || img->channels != 1 )
    return CV_BADFACTOR_ERR;

    roi = _cvSize( img->roi.width, img->roi.height );
    
    return _cvLightingCorrection8uC1R( img->data + img->roi.y * img->step + img->roi.x, 
                                        roi, img->step );

}

*/

/* End of file */