| /*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 |
| // For Open Source Computer Vision Library |
| // |
| // 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*/ |
| |
| |
| #include "_cvaux.h" |
| |
| #if 0 |
| |
| #define LN2PI 1.837877f |
| #define BIG_FLT 1.e+10f |
| |
| |
| #define _CV_ERGODIC 1 |
| #define _CV_CAUSAL 2 |
| |
| #define _CV_LAST_STATE 1 |
| #define _CV_BEST_STATE 2 |
| |
| //*F/////////////////////////////////////////////////////////////////////////////////////// |
| // Name: icvForward1DHMM |
| // Purpose: The function performs baum-welsh algorithm |
| // Context: |
| // Parameters: obs_info - addres of pointer to CvImgObsInfo structure |
| // num_hor_obs - number of horizontal observation vectors |
| // num_ver_obs - number of horizontal observation vectors |
| // obs_size - length of observation vector |
| // |
| // Returns: error status |
| // |
| // Notes: |
| //F*/ |
| #if 0 |
| CvStatus icvForward1DHMM( int num_states, int num_obs, CvMatr64d A, |
| CvMatr64d B, |
| double* scales) |
| { |
| // assume that observation and transition |
| // probabilities already computed |
| int m_HMMType = _CV_CAUSAL; |
| double* m_pi = icvAlloc( num_states* sizeof( double) ); |
| |
| /* alpha is matrix |
| rows throuhg states |
| columns through time |
| */ |
| double* alpha = icvAlloc( num_states*num_obs * sizeof( double ) ); |
| |
| /* All calculations will be in non-logarithmic domain */ |
| |
| /* Initialization */ |
| /* set initial state probabilities */ |
| m_pi[0] = 1; |
| for (i = 1; i < num_states; i++) |
| { |
| m_pi[i] = 0.0; |
| } |
| |
| for (i = 0; i < num_states; i++) |
| { |
| alpha[i] = m_pi[i] * m_b[ i]; |
| } |
| |
| /******************************************************************/ |
| /* Induction */ |
| |
| if ( m_HMMType == _CV_ERGODIC ) |
| { |
| int t; |
| for (t = 1 ; t < num_obs; t++) |
| { |
| for (j = 0; j < num_states; j++) |
| { |
| double sum = 0.0; |
| int i; |
| |
| for (i = 0; i < num_states; i++) |
| { |
| sum += alpha[(t - 1) * num_states + i] * A[i * num_states + j]; |
| } |
| |
| alpha[(t - 1) * num_states + j] = sum * B[t * num_states + j]; |
| |
| /* add computed alpha to scale factor */ |
| sum_alpha += alpha[(t - 1) * num_states + j]; |
| } |
| |
| double scale = 1/sum_alpha; |
| |
| /* scale alpha */ |
| for (j = 0; j < num_states; j++) |
| { |
| alpha[(t - 1) * num_states + j] *= scale; |
| } |
| |
| scales[t] = scale; |
| |
| } |
| } |
| |
| #endif |
| |
| |
| |
| //*F/////////////////////////////////////////////////////////////////////////////////////// |
| // Name: icvCreateObsInfo |
| // Purpose: The function allocates memory for CvImgObsInfo structure |
| // and its inner stuff |
| // Context: |
| // Parameters: obs_info - addres of pointer to CvImgObsInfo structure |
| // num_hor_obs - number of horizontal observation vectors |
| // num_ver_obs - number of horizontal observation vectors |
| // obs_size - length of observation vector |
| // |
| // Returns: error status |
| // |
| // Notes: |
| //F*/ |
| /*CvStatus icvCreateObsInfo( CvImgObsInfo** obs_info, |
| CvSize num_obs, int obs_size ) |
| { |
| int total = num_obs.height * num_obs.width; |
| |
| CvImgObsInfo* obs = (CvImgObsInfo*)icvAlloc( sizeof( CvImgObsInfo) ); |
| |
| obs->obs_x = num_obs.width; |
| obs->obs_y = num_obs.height; |
| |
| obs->obs = (float*)icvAlloc( total * obs_size * sizeof(float) ); |
| |
| obs->state = (int*)icvAlloc( 2 * total * sizeof(int) ); |
| obs->mix = (int*)icvAlloc( total * sizeof(int) ); |
| |
| obs->obs_size = obs_size; |
| |
| obs_info[0] = obs; |
| |
| return CV_NO_ERR; |
| }*/ |
| |
| /*CvStatus icvReleaseObsInfo( CvImgObsInfo** p_obs_info ) |
| { |
| CvImgObsInfo* obs_info = p_obs_info[0]; |
| |
| icvFree( &(obs_info->obs) ); |
| icvFree( &(obs_info->mix) ); |
| icvFree( &(obs_info->state) ); |
| icvFree( &(obs_info) ); |
| |
| p_obs_info[0] = NULL; |
| |
| return CV_NO_ERR; |
| } */ |
| |
| |
| //*F/////////////////////////////////////////////////////////////////////////////////////// |
| // Name: icvCreate1DHMM |
| // Purpose: The function allocates memory for 1-dimensional HMM |
| // and its inner stuff |
| // Context: |
| // Parameters: hmm - addres of pointer to CvEHMM structure |
| // state_number - number of states in HMM |
| // num_mix - number of gaussian mixtures in HMM states |
| // size of array is defined by previous parameter |
| // obs_size - length of observation vectors |
| // |
| // Returns: error status |
| // Notes: |
| //F*/ |
| CvStatus icvCreate1DHMM( CvEHMM** this_hmm, |
| int state_number, int* num_mix, int obs_size ) |
| { |
| int i; |
| int real_states = state_number; |
| |
| CvEHMMState* all_states; |
| CvEHMM* hmm; |
| int total_mix = 0; |
| float* pointers; |
| |
| /* allocate memory for hmm */ |
| hmm = (CvEHMM*)icvAlloc( sizeof(CvEHMM) ); |
| |
| /* set number of superstates */ |
| hmm->num_states = state_number; |
| hmm->level = 0; |
| |
| /* allocate memory for all states */ |
| all_states = (CvEHMMState *)icvAlloc( real_states * sizeof( CvEHMMState ) ); |
| |
| /* assign number of mixtures */ |
| for( i = 0; i < real_states; i++ ) |
| { |
| all_states[i].num_mix = num_mix[i]; |
| } |
| |
| /* compute size of inner of all real states */ |
| for( i = 0; i < real_states; i++ ) |
| { |
| total_mix += num_mix[i]; |
| } |
| /* allocate memory for states stuff */ |
| pointers = (float*)icvAlloc( total_mix * (2/*for mu invvar */ * obs_size + |
| 2/*for weight and log_var_val*/ ) * sizeof( float) ); |
| |
| /* organize memory */ |
| for( i = 0; i < real_states; i++ ) |
| { |
| all_states[i].mu = pointers; pointers += num_mix[i] * obs_size; |
| all_states[i].inv_var = pointers; pointers += num_mix[i] * obs_size; |
| |
| all_states[i].log_var_val = pointers; pointers += num_mix[i]; |
| all_states[i].weight = pointers; pointers += num_mix[i]; |
| } |
| hmm->u.state = all_states; |
| |
| hmm->transP = icvCreateMatrix_32f( hmm->num_states, hmm->num_states ); |
| hmm->obsProb = NULL; |
| |
| /* if all ok - return pointer */ |
| *this_hmm = hmm; |
| return CV_NO_ERR; |
| } |
| |
| CvStatus icvRelease1DHMM( CvEHMM** phmm ) |
| { |
| CvEHMM* hmm = phmm[0]; |
| icvDeleteMatrix( hmm->transP ); |
| |
| if (hmm->obsProb != NULL) |
| { |
| int* tmp = ((int*)(hmm->obsProb)) - 3; |
| icvFree( &(tmp) ); |
| } |
| |
| icvFree( &(hmm->u.state->mu) ); |
| icvFree( &(hmm->u.state) ); |
| |
| phmm[0] = NULL; |
| |
| return CV_NO_ERR; |
| } |
| |
| /*can be used in CHMM & DHMM */ |
| CvStatus icvUniform1DSegm( Cv1DObsInfo* obs_info, CvEHMM* hmm ) |
| { |
| /* implementation is very bad */ |
| int i; |
| CvEHMMState* first_state; |
| |
| /* check arguments */ |
| if ( !obs_info || !hmm ) return CV_NULLPTR_ERR; |
| |
| first_state = hmm->u.state; |
| |
| for (i = 0; i < obs_info->obs_x; i++) |
| { |
| //bad line (division ) |
| int state = (i * hmm->num_states)/obs_info->obs_x; |
| obs_info->state[i] = state; |
| } |
| return CV_NO_ERR; |
| } |
| |
| |
| |
| /*F/////////////////////////////////////////////////////////////////////////////////////// |
| // Name: InitMixSegm |
| // Purpose: The function implements the mixture segmentation of the states of the embedded HMM |
| // Context: used with the Viterbi training of the embedded HMM |
| // Function uses K-Means algorithm for clustering |
| // |
| // Parameters: obs_info_array - array of pointers to image observations |
| // num_img - length of above array |
| // hmm - pointer to HMM structure |
| // |
| // Returns: error status |
| // |
| // Notes: |
| //F*/ |
| CvStatus icvInit1DMixSegm(Cv1DObsInfo** obs_info_array, int num_img, CvEHMM* hmm) |
| { |
| int k, i, j; |
| int* num_samples; /* number of observations in every state */ |
| int* counter; /* array of counters for every state */ |
| |
| int** a_class; /* for every state - characteristic array */ |
| |
| CvVect32f** samples; /* for every state - pointer to observation vectors */ |
| int*** samples_mix; /* for every state - array of pointers to vectors mixtures */ |
| |
| CvTermCriteria criteria = cvTermCriteria( CV_TERMCRIT_EPS|CV_TERMCRIT_ITER, |
| 1000, /* iter */ |
| 0.01f ); /* eps */ |
| |
| int total = hmm->num_states; |
| CvEHMMState* first_state = hmm->u.state; |
| |
| /* for every state integer is allocated - number of vectors in state */ |
| num_samples = (int*)icvAlloc( total * sizeof(int) ); |
| |
| /* integer counter is allocated for every state */ |
| counter = (int*)icvAlloc( total * sizeof(int) ); |
| |
| samples = (CvVect32f**)icvAlloc( total * sizeof(CvVect32f*) ); |
| samples_mix = (int***)icvAlloc( total * sizeof(int**) ); |
| |
| /* clear */ |
| memset( num_samples, 0 , total*sizeof(int) ); |
| memset( counter, 0 , total*sizeof(int) ); |
| |
| |
| /* for every state the number of vectors which belong to it is computed (smth. like histogram) */ |
| for (k = 0; k < num_img; k++) |
| { |
| CvImgObsInfo* obs = obs_info_array[k]; |
| |
| for (i = 0; i < obs->obs_x; i++) |
| { |
| int state = obs->state[ i ]; |
| num_samples[state] += 1; |
| } |
| } |
| |
| /* for every state int* is allocated */ |
| a_class = (int**)icvAlloc( total*sizeof(int*) ); |
| |
| for (i = 0; i < total; i++) |
| { |
| a_class[i] = (int*)icvAlloc( num_samples[i] * sizeof(int) ); |
| samples[i] = (CvVect32f*)icvAlloc( num_samples[i] * sizeof(CvVect32f) ); |
| samples_mix[i] = (int**)icvAlloc( num_samples[i] * sizeof(int*) ); |
| } |
| |
| /* for every state vectors which belong to state are gathered */ |
| for (k = 0; k < num_img; k++) |
| { |
| CvImgObsInfo* obs = obs_info_array[k]; |
| int num_obs = obs->obs_x; |
| float* vector = obs->obs; |
| |
| for (i = 0; i < num_obs; i++, vector+=obs->obs_size ) |
| { |
| int state = obs->state[i]; |
| |
| samples[state][counter[state]] = vector; |
| samples_mix[state][counter[state]] = &(obs->mix[i]); |
| counter[state]++; |
| } |
| } |
| |
| /* clear counters */ |
| memset( counter, 0, total*sizeof(int) ); |
| |
| /* do the actual clustering using the K Means algorithm */ |
| for (i = 0; i < total; i++) |
| { |
| if ( first_state[i].num_mix == 1) |
| { |
| for (k = 0; k < num_samples[i]; k++) |
| { |
| /* all vectors belong to one mixture */ |
| a_class[i][k] = 0; |
| } |
| } |
| else if( num_samples[i] ) |
| { |
| /* clusterize vectors */ |
| icvKMeans( first_state[i].num_mix, samples[i], num_samples[i], |
| obs_info_array[0]->obs_size, criteria, a_class[i] ); |
| } |
| } |
| |
| /* for every vector number of mixture is assigned */ |
| for( i = 0; i < total; i++ ) |
| { |
| for (j = 0; j < num_samples[i]; j++) |
| { |
| samples_mix[i][j][0] = a_class[i][j]; |
| } |
| } |
| |
| for (i = 0; i < total; i++) |
| { |
| icvFree( &(a_class[i]) ); |
| icvFree( &(samples[i]) ); |
| icvFree( &(samples_mix[i]) ); |
| } |
| |
| icvFree( &a_class ); |
| icvFree( &samples ); |
| icvFree( &samples_mix ); |
| icvFree( &counter ); |
| icvFree( &num_samples ); |
| |
| |
| return CV_NO_ERR; |
| } |
| |
| /*F/////////////////////////////////////////////////////////////////////////////////////// |
| // Name: ComputeUniModeGauss |
| // Purpose: The function computes the Gaussian pdf for a sample vector |
| // Context: |
| // Parameters: obsVeq - pointer to the sample vector |
| // mu - pointer to the mean vector of the Gaussian pdf |
| // var - pointer to the variance vector of the Gaussian pdf |
| // VecSize - the size of sample vector |
| // |
| // Returns: the pdf of the sample vector given the specified Gaussian |
| // |
| // Notes: |
| //F*/ |
| /*float icvComputeUniModeGauss(CvVect32f vect, CvVect32f mu, |
| CvVect32f inv_var, float log_var_val, int vect_size) |
| { |
| int n; |
| double tmp; |
| double prob; |
| |
| prob = -log_var_val; |
| |
| for (n = 0; n < vect_size; n++) |
| { |
| tmp = (vect[n] - mu[n]) * inv_var[n]; |
| prob = prob - tmp * tmp; |
| } |
| //prob *= 0.5f; |
| |
| return (float)prob; |
| }*/ |
| |
| /*F/////////////////////////////////////////////////////////////////////////////////////// |
| // Name: ComputeGaussMixture |
| // Purpose: The function computes the mixture Gaussian pdf of a sample vector. |
| // Context: |
| // Parameters: obsVeq - pointer to the sample vector |
| // mu - two-dimensional pointer to the mean vector of the Gaussian pdf; |
| // the first dimension is indexed over the number of mixtures and |
| // the second dimension is indexed along the size of the mean vector |
| // var - two-dimensional pointer to the variance vector of the Gaussian pdf; |
| // the first dimension is indexed over the number of mixtures and |
| // the second dimension is indexed along the size of the variance vector |
| // VecSize - the size of sample vector |
| // weight - pointer to the wights of the Gaussian mixture |
| // NumMix - the number of Gaussian mixtures |
| // |
| // Returns: the pdf of the sample vector given the specified Gaussian mixture. |
| // |
| // Notes: |
| //F*/ |
| /* Calculate probability of observation at state in logarithmic scale*/ |
| /*float icvComputeGaussMixture( CvVect32f vect, float* mu, |
| float* inv_var, float* log_var_val, |
| int vect_size, float* weight, int num_mix ) |
| { |
| double prob, l_prob; |
| |
| prob = 0.0f; |
| |
| if (num_mix == 1) |
| { |
| return icvComputeUniModeGauss( vect, mu, inv_var, log_var_val[0], vect_size); |
| } |
| else |
| { |
| int m; |
| for (m = 0; m < num_mix; m++) |
| { |
| if ( weight[m] > 0.0) |
| { |
| l_prob = icvComputeUniModeGauss(vect, mu + m*vect_size, |
| inv_var + m * vect_size, |
| log_var_val[m], |
| vect_size); |
| |
| prob = prob + weight[m]*exp((double)l_prob); |
| } |
| } |
| prob = log(prob); |
| } |
| return (float)prob; |
| } |
| */ |
| |
| /*F/////////////////////////////////////////////////////////////////////////////////////// |
| // Name: EstimateObsProb |
| // Purpose: The function computes the probability of every observation in every state |
| // Context: |
| // Parameters: obs_info - observations |
| // hmm - hmm |
| // Returns: error status |
| // |
| // Notes: |
| //F*/ |
| CvStatus icvEstimate1DObsProb(CvImgObsInfo* obs_info, CvEHMM* hmm ) |
| { |
| int j; |
| int total_states = 0; |
| |
| /* check if matrix exist and check current size |
| if not sufficient - realloc */ |
| int status = 0; /* 1 - not allocated, 2 - allocated but small size, |
| 3 - size is enough, but distribution is bad, 0 - all ok */ |
| |
| /*for( j = 0; j < hmm->num_states; j++ ) |
| { |
| total_states += hmm->u.ehmm[j].num_states; |
| }*/ |
| total_states = hmm->num_states; |
| |
| if ( hmm->obsProb == NULL ) |
| { |
| /* allocare memory */ |
| int need_size = ( obs_info->obs_x /* * obs_info->obs_y*/ * total_states * sizeof(float) /* + |
| obs_info->obs_y * hmm->num_states * sizeof( CvMatr32f) */); |
| |
| int* buffer = (int*)icvAlloc( need_size + 3 * sizeof(int) ); |
| buffer[0] = need_size; |
| buffer[1] = obs_info->obs_y; |
| buffer[2] = obs_info->obs_x; |
| hmm->obsProb = (float**) (buffer + 3); |
| status = 3; |
| |
| } |
| else |
| { |
| /* check current size */ |
| int* total= (int*)(((int*)(hmm->obsProb)) - 3); |
| int need_size = ( obs_info->obs_x /* * obs_info->obs_y*/ * total_states * sizeof(float) /* + |
| obs_info->obs_y * hmm->num_states * sizeof( CvMatr32f(float*) )*/ ); |
| |
| assert( sizeof(float*) == sizeof(int) ); |
| |
| if ( need_size > (*total) ) |
| { |
| int* buffer = ((int*)(hmm->obsProb)) - 3; |
| icvFree( &buffer); |
| buffer = (int*)icvAlloc( need_size + 3); |
| buffer[0] = need_size; |
| buffer[1] = obs_info->obs_y; |
| buffer[2] = obs_info->obs_x; |
| |
| hmm->obsProb = (float**)(buffer + 3); |
| |
| status = 3; |
| } |
| } |
| if (!status) |
| { |
| int* obsx = ((int*)(hmm->obsProb)) - 1; |
| //int* obsy = ((int*)(hmm->obsProb)) - 2; |
| |
| assert( /*(*obsy > 0) &&*/ (*obsx > 0) ); |
| |
| /* is good distribution? */ |
| if ( (obs_info->obs_x > (*obsx) ) /* || (obs_info->obs_y > (*obsy) ) */ ) |
| status = 3; |
| } |
| |
| assert( (status == 0) || (status == 3) ); |
| /* if bad status - do reallocation actions */ |
| if ( status ) |
| { |
| float** tmp = hmm->obsProb; |
| //float* tmpf; |
| |
| /* distribute pointers of ehmm->obsProb */ |
| /* for( i = 0; i < hmm->num_states; i++ ) |
| { |
| hmm->u.ehmm[i].obsProb = tmp; |
| tmp += obs_info->obs_y; |
| } |
| */ |
| //tmpf = (float*)tmp; |
| |
| /* distribute pointers of ehmm->obsProb[j] */ |
| /* for( i = 0; i < hmm->num_states; i++ ) |
| { |
| CvEHMM* ehmm = &( hmm->u.ehmm[i] ); |
| |
| for( j = 0; j < obs_info->obs_y; j++ ) |
| { |
| ehmm->obsProb[j] = tmpf; |
| tmpf += ehmm->num_states * obs_info->obs_x; |
| } |
| } |
| */ |
| hmm->obsProb = tmp; |
| |
| }/* end of pointer distribution */ |
| |
| #if 1 |
| { |
| #define MAX_BUF_SIZE 1200 |
| float local_log_mix_prob[MAX_BUF_SIZE]; |
| double local_mix_prob[MAX_BUF_SIZE]; |
| int vect_size = obs_info->obs_size; |
| CvStatus res = CV_NO_ERR; |
| |
| float* log_mix_prob = local_log_mix_prob; |
| double* mix_prob = local_mix_prob; |
| |
| int max_size = 0; |
| int obs_x = obs_info->obs_x; |
| |
| /* calculate temporary buffer size */ |
| //for( i = 0; i < hmm->num_states; i++ ) |
| //{ |
| // CvEHMM* ehmm = &(hmm->u.ehmm[i]); |
| CvEHMMState* state = hmm->u.state; |
| |
| int max_mix = 0; |
| for( j = 0; j < hmm->num_states; j++ ) |
| { |
| int t = state[j].num_mix; |
| if( max_mix < t ) max_mix = t; |
| } |
| max_mix *= hmm->num_states; |
| /*if( max_size < max_mix )*/ max_size = max_mix; |
| //} |
| |
| max_size *= obs_x * vect_size; |
| |
| /* allocate buffer */ |
| if( max_size > MAX_BUF_SIZE ) |
| { |
| log_mix_prob = (float*)icvAlloc( max_size*(sizeof(float) + sizeof(double))); |
| if( !log_mix_prob ) return CV_OUTOFMEM_ERR; |
| mix_prob = (double*)(log_mix_prob + max_size); |
| } |
| |
| memset( log_mix_prob, 0, max_size*sizeof(float)); |
| |
| /*****************computing probabilities***********************/ |
| |
| /* loop through external states */ |
| //for( i = 0; i < hmm->num_states; i++ ) |
| { |
| // CvEHMM* ehmm = &(hmm->u.ehmm[i]); |
| CvEHMMState* state = hmm->u.state; |
| |
| int max_mix = 0; |
| int n_states = hmm->num_states; |
| |
| /* determine maximal number of mixtures (again) */ |
| for( j = 0; j < hmm->num_states; j++ ) |
| { |
| int t = state[j].num_mix; |
| if( max_mix < t ) max_mix = t; |
| } |
| |
| /* loop through rows of the observation matrix */ |
| //for( j = 0; j < obs_info->obs_y; j++ ) |
| { |
| int m, n; |
| |
| float* obs = obs_info->obs;/* + j * obs_x * vect_size; */ |
| float* log_mp = max_mix > 1 ? log_mix_prob : (float*)(hmm->obsProb); |
| double* mp = mix_prob; |
| |
| /* several passes are done below */ |
| |
| /* 1. calculate logarithms of probabilities for each mixture */ |
| |
| /* loop through mixtures */ |
| /* !!!! */ for( m = 0; m < max_mix; m++ ) |
| { |
| /* set pointer to first observation in the line */ |
| float* vect = obs; |
| |
| /* cycles through obseravtions in the line */ |
| for( n = 0; n < obs_x; n++, vect += vect_size, log_mp += n_states ) |
| { |
| int k, l; |
| for( l = 0; l < n_states; l++ ) |
| { |
| if( state[l].num_mix > m ) |
| { |
| float* mu = state[l].mu + m*vect_size; |
| float* inv_var = state[l].inv_var + m*vect_size; |
| double prob = -state[l].log_var_val[m]; |
| for( k = 0; k < vect_size; k++ ) |
| { |
| double t = (vect[k] - mu[k])*inv_var[k]; |
| prob -= t*t; |
| } |
| log_mp[l] = MAX( (float)prob, -500 ); |
| } |
| } |
| } |
| } |
| |
| /* skip the rest if there is a single mixture */ |
| if( max_mix != 1 ) |
| { |
| /* 2. calculate exponent of log_mix_prob |
| (i.e. probability for each mixture) */ |
| res = icvbExp_32f64f( log_mix_prob, mix_prob, |
| max_mix * obs_x * n_states ); |
| if( res < 0 ) goto processing_exit; |
| |
| /* 3. sum all mixtures with weights */ |
| /* 3a. first mixture - simply scale by weight */ |
| for( n = 0; n < obs_x; n++, mp += n_states ) |
| { |
| int l; |
| for( l = 0; l < n_states; l++ ) |
| { |
| mp[l] *= state[l].weight[0]; |
| } |
| } |
| |
| /* 3b. add other mixtures */ |
| for( m = 1; m < max_mix; m++ ) |
| { |
| int ofs = -m*obs_x*n_states; |
| for( n = 0; n < obs_x; n++, mp += n_states ) |
| { |
| int l; |
| for( l = 0; l < n_states; l++ ) |
| { |
| if( m < state[l].num_mix ) |
| { |
| mp[l + ofs] += mp[l] * state[l].weight[m]; |
| } |
| } |
| } |
| } |
| |
| /* 4. Put logarithms of summary probabilities to the destination matrix */ |
| 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 |
| |
