cvhmm1d.cpp.svn-base

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

                    res = icvbLog_64f32f( mix_prob, (float*)(hmm->obsProb),//[j],
                                            obs_x * n_states );
                    if( res < 0 ) goto processing_exit;
                }
            }
        }

processing_exit:

        if( log_mix_prob != local_log_mix_prob ) icvFree( &log_mix_prob );
        return res;
#undef MAX_BUF_SIZE
    }
#else
/*    for( i = 0; i < hmm->num_states; i++ )
    {
        CvEHMM* ehmm = &(hmm->u.ehmm[i]);
        CvEHMMState* state = ehmm->u.state;

        for( j = 0; j < obs_info->obs_y; j++ )
        {
            int k,m;
                       
            int obs_index = j * obs_info->obs_x;

            float* B = ehmm->obsProb[j];
            
            // cycles through obs and states
            for( k = 0; k < obs_info->obs_x; k++ )
            {
                CvVect32f vect = (obs_info->obs) + (obs_index + k) * vect_size;
                
                float* matr_line = B + k * ehmm->num_states;

                for( m = 0; m < ehmm->num_states; m++ )
                {
                    matr_line[m] = icvComputeGaussMixture( vect, state[m].mu, state[m].inv_var, 
                                                             state[m].log_var_val, vect_size, state[m].weight,
                                                             state[m].num_mix );
                }
            }
        }
    }
*/
#endif
}


/*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: EstimateTransProb
//    Purpose: The function calculates the state and super state transition probabilities
//             of the model given the images,
//             the state segmentation and the input parameters
//    Context:
//    Parameters: obs_info_array - array of pointers to image observations
//                num_img - length of above array
//                hmm - pointer to HMM structure
//    Returns: void
//
//    Notes:
//F*/
CvStatus icvEstimate1DTransProb( Cv1DObsInfo** obs_info_array,
                                 int num_seq,
                                 CvEHMM* hmm )
{
    int    i, j, k;

    /* as a counter we will use transP matrix */

    /* initialization */

    /* clear transP */
    icvSetZero_32f( hmm->transP, hmm->num_states, hmm->num_states );


    /* compute the counters */
    for (i = 0; i < num_seq; i++)
    {
        int counter = 0;
        Cv1DObsInfo* info = obs_info_array[i];

        for (k = 0; k < info->obs_x; k++, counter++)
        {
            /* compute how many transitions from state to state
               occured */ 
            int state;
            int nextstate;
            
            state = info->state[counter];

            if (k < info->obs_x - 1)
            {   
                int transP_size = hmm->num_states;

                nextstate = info->state[counter+1];
                hmm->transP[ state * transP_size + nextstate] += 1;
            }            
        }
    }

    /* estimate superstate matrix */
    for( i = 0; i < hmm->num_states; i++)
    {
        float total = 0;
        float inv_total;
        for( j = 0; j < hmm->num_states; j++)
        {
            total += hmm->transP[i * hmm->num_states + j];
        }
        //assert( total );

        inv_total = total ? 1.f/total : 0;
        
        for( j = 0; j < hmm->num_states; j++)
        {                   
            hmm->transP[i * hmm->num_states + j] = 
                hmm->transP[i * hmm->num_states + j] ? 
                (float)log( hmm->transP[i * hmm->num_states + j] * inv_total ) : -BIG_FLT;
        }
    }
    
    return CV_NO_ERR;
} 
                      

/*F///////////////////////////////////////////////////////////////////////////////////////
//    Name: MixSegmL2
//    Purpose: The function implements the mixture segmentation of the states of the embedded HMM
//    Context: used with the Viterbi training of the embedded HMM
//
//    Parameters:  
//             obs_info_array
//             num_img
//             hmm
//    Returns: void
//
//    Notes: 
//F*/
CvStatus icv1DMixSegmL2(CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm )
{
    int     k, i, m;
       
    CvEHMMState* state = hmm->u.state;
    
    for (k = 0; k < num_img; k++)
    {   
        //int counter = 0;
        CvImgObsInfo* info = obs_info_array[k];

        for (i = 0; i < info->obs_x; i++)
        {
            int e_state = info->state[i];
            float min_dist;
                                                
            min_dist = icvSquareDistance((info->obs) + (i * info->obs_size), 
                                               state[e_state].mu, info->obs_size);
            info->mix[i] = 0;  
                
            for (m = 1; m < state[e_state].num_mix; m++)
            {                 
                float dist=icvSquareDistance( (info->obs) + (i * info->obs_size),
                                               state[e_state].mu + m * info->obs_size,
                                               info->obs_size);
                if (dist < min_dist)
                {
                    min_dist = dist;
                    /* assign mixture with smallest distance */ 
                    info->mix[i] = m;
                }
            }
        }
    }
    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*/
float icvViterbi(Cv1DObsInfo* obs_info, CvEHMM* hmm)
{
    int    i, counter;
    float  log_likelihood;

    //CvEHMMState* first_state = hmm->u.state;
    
    /* memory allocation for superB */
    /*CvMatr32f superB = picvCreateMatrix_32f(hmm->num_states, obs_info->obs_x );*/
    
    /* memory allocation for q */
    int* super_q = (int*)icvAlloc( obs_info->obs_x * sizeof(int) );
    
    /* perform Viterbi segmentation (process 1D HMM) */
    icvViterbiSegmentation( hmm->num_states, obs_info->obs_x, 
                            hmm->transP, (float*)(hmm->obsProb), 0, 
                            _CV_LAST_STATE, &super_q, obs_info->obs_x,
                             obs_info->obs_x, &log_likelihood );
    
    log_likelihood /= obs_info->obs_x ;   

    counter = 0;
    /* assign new state to observation vectors */
    for (i = 0; i < obs_info->obs_x; i++)
    {   
         int state = super_q[i];
         obs_info->state[i] = state;
    }
    
    /* memory deallocation for superB */
    /*picvDeleteMatrix( superB );*/
    icvFree( &super_q );
    
    return log_likelihood;
}  

CvStatus icvEstimate1DHMMStateParams(CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm)

{
    /* compute gamma, weights, means, vars */
    int k, i, j, m;
    int counter = 0;
    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.state;

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

    total+= hmm->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 */
    counter = 0;
    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[i];
            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[j];
            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 */
    counter = 0;
    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 - 0.01 (Ara's experimental result) */
            for( k = 0; k < vect_len; k++ )
            {
                invar[k] = (invar[k] > 0.01f) ? invar[k] : 0.01f;
            }

            /* 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] );
            }    
                        
            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; 
}





#endif