1 /*M///////////////////////////////////////////////////////////////////////////////////////
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install,
7 // copy or use the software.
10 // Intel License Agreement
11 // For Open Source Computer Vision Library
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Third party copyrights are property of their respective owners.
16 // Redistribution and use in source and binary forms, with or without modification,
17 // are permitted provided that the following conditions are met:
19 // * Redistribution's of source code must retain the above copyright notice,
20 // this list of conditions and the following disclaimer.
22 // * Redistribution's in binary form must reproduce the above copyright notice,
23 // this list of conditions and the following disclaimer in the documentation
24 // and/or other materials provided with the distribution.
26 // * The name of Intel Corporation may not be used to endorse or promote products
27 // derived from this software without specific prior written permission.
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the Intel Corporation or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
44 #define LN2PI 1.837877f
45 #define BIG_FLT 1.e+10f
51 #define _CV_LAST_STATE 1
52 #define _CV_BEST_STATE 2
55 //*F///////////////////////////////////////////////////////////////////////////////////////
56 // Name: _cvCreateObsInfo
57 // Purpose: The function allocates memory for CvImgObsInfo structure
58 // and its inner stuff
60 // Parameters: obs_info - addres of pointer to CvImgObsInfo structure
61 // num_hor_obs - number of horizontal observation vectors
62 // num_ver_obs - number of horizontal observation vectors
63 // obs_size - length of observation vector
65 // Returns: error status
69 static CvStatus CV_STDCALL icvCreateObsInfo( CvImgObsInfo** obs_info,
70 CvSize num_obs, int obs_size )
72 int total = num_obs.height * num_obs.width;
74 CvImgObsInfo* obs = (CvImgObsInfo*)cvAlloc( sizeof( CvImgObsInfo) );
76 obs->obs_x = num_obs.width;
77 obs->obs_y = num_obs.height;
79 obs->obs = (float*)cvAlloc( total * obs_size * sizeof(float) );
81 obs->state = (int*)cvAlloc( 2 * total * sizeof(int) );
82 obs->mix = (int*)cvAlloc( total * sizeof(int) );
84 obs->obs_size = obs_size;
91 static CvStatus CV_STDCALL icvReleaseObsInfo( CvImgObsInfo** p_obs_info )
93 CvImgObsInfo* obs_info = p_obs_info[0];
95 cvFree( &(obs_info->obs) );
96 cvFree( &(obs_info->mix) );
97 cvFree( &(obs_info->state) );
98 cvFree( &(obs_info) );
100 p_obs_info[0] = NULL;
106 //*F///////////////////////////////////////////////////////////////////////////////////////
107 // Name: icvCreate2DHMM
108 // Purpose: The function allocates memory for 2-dimensional embedded HMM model
109 // and its inner stuff
111 // Parameters: hmm - addres of pointer to CvEHMM structure
112 // state_number - array of hmm sizes (size of array == state_number[0]+1 )
113 // num_mix - number of gaussian mixtures in low-level HMM states
114 // size of array is defined by previous array values
115 // obs_size - length of observation vectors
117 // Returns: error status
119 // Notes: state_number[0] - number of states in external HMM.
120 // state_number[i] - number of states in embedded HMM
122 // example for face recognition: state_number = { 5 3 6 6 6 3 },
123 // length of num_mix array = 3+6+6+6+3 = 24//
126 static CvStatus CV_STDCALL icvCreate2DHMM( CvEHMM** this_hmm,
127 int* state_number, int* num_mix, int obs_size )
132 CvEHMMState* all_states;
137 //compute total number of states of all level in 2d EHMM
138 for( i = 1; i <= state_number[0]; i++ )
140 real_states += state_number[i];
143 /* allocate memory for all hmms (from all levels) */
144 hmm = (CvEHMM*)cvAlloc( (state_number[0] + 1) * sizeof(CvEHMM) );
146 /* set number of superstates */
147 hmm[0].num_states = state_number[0];
150 /* allocate memory for all states */
151 all_states = (CvEHMMState *)cvAlloc( real_states * sizeof( CvEHMMState ) );
153 /* assign number of mixtures */
154 for( i = 0; i < real_states; i++ )
156 all_states[i].num_mix = num_mix[i];
159 /* compute size of inner of all real states */
160 for( i = 0; i < real_states; i++ )
162 total_mix += num_mix[i];
164 /* allocate memory for states stuff */
165 pointers = (float*)cvAlloc( total_mix * (2/*for mu invvar */ * obs_size +
166 2/*for weight and log_var_val*/ ) * sizeof( float) );
168 /* organize memory */
169 for( i = 0; i < real_states; i++ )
171 all_states[i].mu = pointers; pointers += num_mix[i] * obs_size;
172 all_states[i].inv_var = pointers; pointers += num_mix[i] * obs_size;
174 all_states[i].log_var_val = pointers; pointers += num_mix[i];
175 all_states[i].weight = pointers; pointers += num_mix[i];
178 /* set pointer to embedded hmm array */
179 hmm->u.ehmm = hmm + 1;
181 for( i = 0; i < hmm[0].num_states; i++ )
183 hmm[i+1].u.state = all_states;
184 all_states += state_number[i+1];
185 hmm[i+1].num_states = state_number[i+1];
188 for( i = 0; i <= state_number[0]; i++ )
190 hmm[i].transP = icvCreateMatrix_32f( hmm[i].num_states, hmm[i].num_states );
191 hmm[i].obsProb = NULL;
192 hmm[i].level = i ? 0 : 1;
195 /* if all ok - return pointer */
200 static CvStatus CV_STDCALL icvRelease2DHMM( CvEHMM** phmm )
202 CvEHMM* hmm = phmm[0];
204 for( i = 0; i < hmm[0].num_states + 1; i++ )
206 icvDeleteMatrix( hmm[i].transP );
209 if (hmm->obsProb != NULL)
211 int* tmp = ((int*)(hmm->obsProb)) - 3;
215 cvFree( &(hmm->u.ehmm->u.state->mu) );
216 cvFree( &(hmm->u.ehmm->u.state) );
219 /* free hmm structures */
227 /* distance between 2 vectors */
228 static float icvSquareDistance( CvVect32f v1, CvVect32f v2, int len )
234 for( i = 0; i <= len - 4; i += 4 )
236 double t0 = v1[i] - v2[i];
237 double t1 = v1[i+1] - v2[i+1];
241 t0 = v1[i+2] - v2[i+2];
242 t1 = v1[i+3] - v2[i+3];
247 for( ; i < len; i++ )
249 double t0 = v1[i] - v2[i];
253 return (float)(dist0 + dist1);
256 /*can be used in CHMM & DHMM */
257 static CvStatus CV_STDCALL
258 icvUniformImgSegm( CvImgObsInfo* obs_info, CvEHMM* hmm )
261 /* implementation is very bad */
262 int i, j, counter = 0;
263 CvEHMMState* first_state;
264 float inv_x = 1.f/obs_info->obs_x;
265 float inv_y = 1.f/obs_info->obs_y;
267 /* check arguments */
268 if ( !obs_info || !hmm ) return CV_NULLPTR_ERR;
270 first_state = hmm->u.ehmm->u.state;
272 for (i = 0; i < obs_info->obs_y; i++)
274 //bad line (division )
275 int superstate = (int)((i * hmm->num_states)*inv_y);/* /obs_info->obs_y; */
277 int index = (int)(hmm->u.ehmm[superstate].u.state - first_state);
279 for (j = 0; j < obs_info->obs_x; j++, counter++)
281 int state = (int)((j * hmm->u.ehmm[superstate].num_states)* inv_x); /* / obs_info->obs_x; */
283 obs_info->state[2 * counter] = superstate;
284 obs_info->state[2 * counter + 1] = state + index;
288 //this is not ready yet
291 CvEHMMState* first_state = hmm->u.ehmm->u.state;
293 /* check bad arguments */
294 if ( hmm->num_states > obs_info->obs_y ) return CV_BADSIZE_ERR;
296 //compute vertical subdivision
297 float row_per_state = (float)obs_info->obs_y / hmm->num_states;
298 float col_per_state[1024]; /* maximum 1024 superstates */
300 //for every horizontal band compute subdivision
301 for( i = 0; i < hmm->num_states; i++ )
303 CvEHMM* ehmm = &(hmm->u.ehmm[i]);
304 col_per_state[i] = (float)obs_info->obs_x / ehmm->num_states;
307 //compute state bounds
309 for( i = 0; i < hmm->num_states - 1; i++ )
311 ss_bound[i] = floor( row_per_state * ( i+1 ) );
313 ss_bound[hmm->num_states - 1] = obs_info->obs_y;
315 //work inside every superstate
319 for( i = 0; i < hmm->num_states; i++ )
321 CvEHMM* ehmm = &(hmm->u.ehmm[i]);
322 int index = ehmm->u.state - first_state;
324 //calc distribution in superstate
326 for( j = 0; j < ehmm->num_states - 1; j++ )
328 es_bound[j] = floor( col_per_state[i] * ( j+1 ) );
330 es_bound[ehmm->num_states - 1] = obs_info->obs_x;
332 //assign states to first row of superstate
334 for( j = 0; j < ehmm->num_states; j++ )
336 for( k = col; k < es_bound[j]; k++, col++ )
338 obs_info->state[row * obs_info->obs_x + 2 * k] = i;
339 obs_info->state[row * obs_info->obs_x + 2 * k + 1] = j + index;
344 //copy the same to other rows of superstate
345 for( m = row; m < ss_bound[i]; m++ )
347 memcpy( &(obs_info->state[m * obs_info->obs_x * 2]),
348 &(obs_info->state[row * obs_info->obs_x * 2]), obs_info->obs_x * 2 * sizeof(int) );
360 /*F///////////////////////////////////////////////////////////////////////////////////////
362 // Purpose: The function implements the mixture segmentation of the states of the
364 // Context: used with the Viterbi training of the embedded HMM
365 // Function uses K-Means algorithm for clustering
367 // Parameters: obs_info_array - array of pointers to image observations
368 // num_img - length of above array
369 // hmm - pointer to HMM structure
371 // Returns: error status
375 static CvStatus CV_STDCALL
376 icvInitMixSegm( CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm )
379 int* num_samples; /* number of observations in every state */
380 int* counter; /* array of counters for every state */
382 int** a_class; /* for every state - characteristic array */
384 CvVect32f** samples; /* for every state - pointer to observation vectors */
385 int*** samples_mix; /* for every state - array of pointers to vectors mixtures */
387 CvTermCriteria criteria = cvTermCriteria( CV_TERMCRIT_EPS|CV_TERMCRIT_ITER,
393 CvEHMMState* first_state = hmm->u.ehmm->u.state;
395 for( i = 0 ; i < hmm->num_states; i++ )
397 total += hmm->u.ehmm[i].num_states;
400 /* for every state integer is allocated - number of vectors in state */
401 num_samples = (int*)cvAlloc( total * sizeof(int) );
403 /* integer counter is allocated for every state */
404 counter = (int*)cvAlloc( total * sizeof(int) );
406 samples = (CvVect32f**)cvAlloc( total * sizeof(CvVect32f*) );
407 samples_mix = (int***)cvAlloc( total * sizeof(int**) );
410 memset( num_samples, 0 , total*sizeof(int) );
411 memset( counter, 0 , total*sizeof(int) );
414 /* for every state the number of vectors which belong to it is computed (smth. like histogram) */
415 for (k = 0; k < num_img; k++)
417 CvImgObsInfo* obs = obs_info_array[k];
420 for (i = 0; i < obs->obs_y; i++)
422 for (j = 0; j < obs->obs_x; j++, count++)
424 int state = obs->state[ 2 * count + 1];
425 num_samples[state] += 1;
430 /* for every state int* is allocated */
431 a_class = (int**)cvAlloc( total*sizeof(int*) );
433 for (i = 0; i < total; i++)
435 a_class[i] = (int*)cvAlloc( num_samples[i] * sizeof(int) );
436 samples[i] = (CvVect32f*)cvAlloc( num_samples[i] * sizeof(CvVect32f) );
437 samples_mix[i] = (int**)cvAlloc( num_samples[i] * sizeof(int*) );
440 /* for every state vectors which belong to state are gathered */
441 for (k = 0; k < num_img; k++)
443 CvImgObsInfo* obs = obs_info_array[k];
444 int num_obs = ( obs->obs_x ) * ( obs->obs_y );
445 float* vector = obs->obs;
447 for (i = 0; i < num_obs; i++, vector+=obs->obs_size )
449 int state = obs->state[2*i+1];
451 samples[state][counter[state]] = vector;
452 samples_mix[state][counter[state]] = &(obs->mix[i]);
458 memset( counter, 0, total*sizeof(int) );
460 /* do the actual clustering using the K Means algorithm */
461 for (i = 0; i < total; i++)
463 if ( first_state[i].num_mix == 1)
465 for (k = 0; k < num_samples[i]; k++)
467 /* all vectors belong to one mixture */
471 else if( num_samples[i] )
473 /* clusterize vectors */
474 cvKMeans( first_state[i].num_mix, samples[i], num_samples[i],
475 obs_info_array[0]->obs_size, criteria, a_class[i] );
479 /* for every vector number of mixture is assigned */
480 for( i = 0; i < total; i++ )
482 for (j = 0; j < num_samples[i]; j++)
484 samples_mix[i][j][0] = a_class[i][j];
488 for (i = 0; i < total; i++)
490 cvFree( &(a_class[i]) );
491 cvFree( &(samples[i]) );
492 cvFree( &(samples_mix[i]) );
497 cvFree( &samples_mix );
499 cvFree( &num_samples );
504 /*F///////////////////////////////////////////////////////////////////////////////////////
505 // Name: ComputeUniModeGauss
506 // Purpose: The function computes the Gaussian pdf for a sample vector
508 // Parameters: obsVeq - pointer to the sample vector
509 // mu - pointer to the mean vector of the Gaussian pdf
510 // var - pointer to the variance vector of the Gaussian pdf
511 // VecSize - the size of sample vector
513 // Returns: the pdf of the sample vector given the specified Gaussian
517 /*static float icvComputeUniModeGauss(CvVect32f vect, CvVect32f mu,
518 CvVect32f inv_var, float log_var_val, int vect_size)
526 for (n = 0; n < vect_size; n++)
528 tmp = (vect[n] - mu[n]) * inv_var[n];
529 prob = prob - tmp * tmp;
536 /*F///////////////////////////////////////////////////////////////////////////////////////
537 // Name: ComputeGaussMixture
538 // Purpose: The function computes the mixture Gaussian pdf of a sample vector.
540 // Parameters: obsVeq - pointer to the sample vector
541 // mu - two-dimensional pointer to the mean vector of the Gaussian pdf;
542 // the first dimension is indexed over the number of mixtures and
543 // the second dimension is indexed along the size of the mean vector
544 // var - two-dimensional pointer to the variance vector of the Gaussian pdf;
545 // the first dimension is indexed over the number of mixtures and
546 // the second dimension is indexed along the size of the variance vector
547 // VecSize - the size of sample vector
548 // weight - pointer to the wights of the Gaussian mixture
549 // NumMix - the number of Gaussian mixtures
551 // Returns: the pdf of the sample vector given the specified Gaussian mixture.
555 /* Calculate probability of observation at state in logarithmic scale*/
557 icvComputeGaussMixture( CvVect32f vect, float* mu,
558 float* inv_var, float* log_var_val,
559 int vect_size, float* weight, int num_mix )
567 return icvComputeUniModeGauss( vect, mu, inv_var, log_var_val[0], vect_size);
572 for (m = 0; m < num_mix; m++)
574 if ( weight[m] > 0.0)
576 l_prob = icvComputeUniModeGauss(vect, mu + m*vect_size,
577 inv_var + m * vect_size,
581 prob = prob + weight[m]*exp((double)l_prob);
590 /*F///////////////////////////////////////////////////////////////////////////////////////
591 // Name: EstimateObsProb
592 // Purpose: The function computes the probability of every observation in every state
594 // Parameters: obs_info - observations
596 // Returns: error status
600 static CvStatus CV_STDCALL icvEstimateObsProb( CvImgObsInfo* obs_info, CvEHMM* hmm )
603 int total_states = 0;
605 /* check if matrix exist and check current size
606 if not sufficient - realloc */
607 int status = 0; /* 1 - not allocated, 2 - allocated but small size,
608 3 - size is enough, but distribution is bad, 0 - all ok */
610 for( j = 0; j < hmm->num_states; j++ )
612 total_states += hmm->u.ehmm[j].num_states;
615 if ( hmm->obsProb == NULL )
617 /* allocare memory */
618 int need_size = ( obs_info->obs_x * obs_info->obs_y * total_states * sizeof(float) +
619 obs_info->obs_y * hmm->num_states * sizeof( CvMatr32f) );
621 int* buffer = (int*)cvAlloc( need_size + 3 * sizeof(int) );
622 buffer[0] = need_size;
623 buffer[1] = obs_info->obs_y;
624 buffer[2] = obs_info->obs_x;
625 hmm->obsProb = (float**) (buffer + 3);
631 /* check current size */
632 int* total= (int*)(((int*)(hmm->obsProb)) - 3);
633 int need_size = ( obs_info->obs_x * obs_info->obs_y * total_states * sizeof(float) +
634 obs_info->obs_y * hmm->num_states * sizeof( CvMatr32f/*(float*)*/ ) );
636 assert( sizeof(float*) == sizeof(int) );
638 if ( need_size > (*total) )
640 int* buffer = ((int*)(hmm->obsProb)) - 3;
642 buffer = (int*)cvAlloc( need_size + 3 * sizeof(int));
643 buffer[0] = need_size;
644 buffer[1] = obs_info->obs_y;
645 buffer[2] = obs_info->obs_x;
647 hmm->obsProb = (float**)(buffer + 3);
654 int* obsx = ((int*)(hmm->obsProb)) - 1;
655 int* obsy = ((int*)(hmm->obsProb)) - 2;
657 assert( (*obsx > 0) && (*obsy > 0) );
659 /* is good distribution? */
660 if ( (obs_info->obs_x > (*obsx) ) || (obs_info->obs_y > (*obsy) ) )
664 /* if bad status - do reallocation actions */
665 assert( (status == 0) || (status == 3) );
669 float** tmp = hmm->obsProb;
672 /* distribute pointers of ehmm->obsProb */
673 for( i = 0; i < hmm->num_states; i++ )
675 hmm->u.ehmm[i].obsProb = tmp;
676 tmp += obs_info->obs_y;
681 /* distribute pointers of ehmm->obsProb[j] */
682 for( i = 0; i < hmm->num_states; i++ )
684 CvEHMM* ehmm = &( hmm->u.ehmm[i] );
686 for( j = 0; j < obs_info->obs_y; j++ )
688 ehmm->obsProb[j] = tmpf;
689 tmpf += ehmm->num_states * obs_info->obs_x;
692 }/* end of pointer distribution */
696 #define MAX_BUF_SIZE 1200
697 float local_log_mix_prob[MAX_BUF_SIZE];
698 double local_mix_prob[MAX_BUF_SIZE];
699 int vect_size = obs_info->obs_size;
700 CvStatus res = CV_NO_ERR;
702 float* log_mix_prob = local_log_mix_prob;
703 double* mix_prob = local_mix_prob;
706 int obs_x = obs_info->obs_x;
708 /* calculate temporary buffer size */
709 for( i = 0; i < hmm->num_states; i++ )
711 CvEHMM* ehmm = &(hmm->u.ehmm[i]);
712 CvEHMMState* state = ehmm->u.state;
715 for( j = 0; j < ehmm->num_states; j++ )
717 int t = state[j].num_mix;
718 if( max_mix < t ) max_mix = t;
720 max_mix *= ehmm->num_states;
721 if( max_size < max_mix ) max_size = max_mix;
724 max_size *= obs_x * vect_size;
726 /* allocate buffer */
727 if( max_size > MAX_BUF_SIZE )
729 log_mix_prob = (float*)cvAlloc( max_size*(sizeof(float) + sizeof(double)));
730 if( !log_mix_prob ) return CV_OUTOFMEM_ERR;
731 mix_prob = (double*)(log_mix_prob + max_size);
734 memset( log_mix_prob, 0, max_size*sizeof(float));
736 /*****************computing probabilities***********************/
738 /* loop through external states */
739 for( i = 0; i < hmm->num_states; i++ )
741 CvEHMM* ehmm = &(hmm->u.ehmm[i]);
742 CvEHMMState* state = ehmm->u.state;
745 int n_states = ehmm->num_states;
747 /* determine maximal number of mixtures (again) */
748 for( j = 0; j < ehmm->num_states; j++ )
750 int t = state[j].num_mix;
751 if( max_mix < t ) max_mix = t;
754 /* loop through rows of the observation matrix */
755 for( j = 0; j < obs_info->obs_y; j++ )
759 float* obs = obs_info->obs + j * obs_x * vect_size;
760 float* log_mp = max_mix > 1 ? log_mix_prob : ehmm->obsProb[j];
761 double* mp = mix_prob;
763 /* several passes are done below */
765 /* 1. calculate logarithms of probabilities for each mixture */
767 /* loop through mixtures */
768 for( m = 0; m < max_mix; m++ )
770 /* set pointer to first observation in the line */
773 /* cycles through obseravtions in the line */
774 for( n = 0; n < obs_x; n++, vect += vect_size, log_mp += n_states )
777 for( l = 0; l < n_states; l++ )
779 if( state[l].num_mix > m )
781 float* mu = state[l].mu + m*vect_size;
782 float* inv_var = state[l].inv_var + m*vect_size;
783 double prob = -state[l].log_var_val[m];
784 for( k = 0; k < vect_size; k++ )
786 double t = (vect[k] - mu[k])*inv_var[k];
789 log_mp[l] = MAX( (float)prob, -500 );
795 /* skip the rest if there is a single mixture */
796 if( max_mix == 1 ) continue;
798 /* 2. calculate exponent of log_mix_prob
799 (i.e. probability for each mixture) */
800 cvbFastExp( log_mix_prob, mix_prob, max_mix * obs_x * n_states );
802 /* 3. sum all mixtures with weights */
803 /* 3a. first mixture - simply scale by weight */
804 for( n = 0; n < obs_x; n++, mp += n_states )
807 for( l = 0; l < n_states; l++ )
809 mp[l] *= state[l].weight[0];
813 /* 3b. add other mixtures */
814 for( m = 1; m < max_mix; m++ )
816 int ofs = -m*obs_x*n_states;
817 for( n = 0; n < obs_x; n++, mp += n_states )
820 for( l = 0; l < n_states; l++ )
822 if( m < state[l].num_mix )
824 mp[l + ofs] += mp[l] * state[l].weight[m];
830 /* 4. Put logarithms of summary probabilities to the destination matrix */
831 cvbFastLog( mix_prob, ehmm->obsProb[j], obs_x * n_states );
835 if( log_mix_prob != local_log_mix_prob ) cvFree( &log_mix_prob );
840 for( i = 0; i < hmm->num_states; i++ )
842 CvEHMM* ehmm = &(hmm->u.ehmm[i]);
843 CvEHMMState* state = ehmm->u.state;
845 for( j = 0; j < obs_info->obs_y; j++ )
849 int obs_index = j * obs_info->obs_x;
851 float* B = ehmm->obsProb[j];
853 /* cycles through obs and states */
854 for( k = 0; k < obs_info->obs_x; k++ )
856 CvVect32f vect = (obs_info->obs) + (obs_index + k) * vect_size;
858 float* matr_line = B + k * ehmm->num_states;
860 for( m = 0; m < ehmm->num_states; m++ )
862 matr_line[m] = icvComputeGaussMixture( vect, state[m].mu, state[m].inv_var,
863 state[m].log_var_val, vect_size, state[m].weight,
873 /*F///////////////////////////////////////////////////////////////////////////////////////
874 // Name: EstimateTransProb
875 // Purpose: The function calculates the state and super state transition probabilities
876 // of the model given the images,
877 // the state segmentation and the input parameters
879 // Parameters: obs_info_array - array of pointers to image observations
880 // num_img - length of above array
881 // hmm - pointer to HMM structure
886 static CvStatus CV_STDCALL
887 icvEstimateTransProb( CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm )
891 CvEHMMState* first_state = hmm->u.ehmm->u.state;
892 /* as a counter we will use transP matrix */
897 icvSetZero_32f( hmm->transP, hmm->num_states, hmm->num_states );
898 for (i = 0; i < hmm->num_states; i++ )
900 icvSetZero_32f( hmm->u.ehmm[i].transP , hmm->u.ehmm[i].num_states, hmm->u.ehmm[i].num_states );
903 /* compute the counters */
904 for (i = 0; i < num_img; i++)
907 CvImgObsInfo* info = obs_info_array[i];
909 for (j = 0; j < info->obs_y; j++)
911 for (k = 0; k < info->obs_x; k++, counter++)
913 /* compute how many transitions from state to state
914 occured both in horizontal and vertical direction */
915 int superstate, state;
916 int nextsuperstate, nextstate;
919 superstate = info->state[2 * counter];
920 begin_ind = (int)(hmm->u.ehmm[superstate].u.state - first_state);
921 state = info->state[ 2 * counter + 1] - begin_ind;
923 if (j < info->obs_y - 1)
925 int transP_size = hmm->num_states;
927 nextsuperstate = info->state[ 2*(counter + info->obs_x) ];
929 hmm->transP[superstate * transP_size + nextsuperstate] += 1;
932 if (k < info->obs_x - 1)
934 int transP_size = hmm->u.ehmm[superstate].num_states;
936 nextstate = info->state[2*(counter+1) + 1] - begin_ind;
937 hmm->u.ehmm[superstate].transP[ state * transP_size + nextstate] += 1;
942 /* estimate superstate matrix */
943 for( i = 0; i < hmm->num_states; i++)
947 for( j = 0; j < hmm->num_states; j++)
949 total += hmm->transP[i * hmm->num_states + j];
953 inv_total = total ? 1.f/total : 0;
955 for( j = 0; j < hmm->num_states; j++)
957 hmm->transP[i * hmm->num_states + j] =
958 hmm->transP[i * hmm->num_states + j] ?
959 (float)log( hmm->transP[i * hmm->num_states + j] * inv_total ) : -BIG_FLT;
963 /* estimate other matrices */
964 for( k = 0; k < hmm->num_states; k++ )
966 CvEHMM* ehmm = &(hmm->u.ehmm[k]);
968 for( i = 0; i < ehmm->num_states; i++)
972 for( j = 0; j < ehmm->num_states; j++)
974 total += ehmm->transP[i*ehmm->num_states + j];
977 inv_total = total ? 1.f/total : 0;
979 for( j = 0; j < ehmm->num_states; j++)
981 ehmm->transP[i * ehmm->num_states + j] =
982 (ehmm->transP[i * ehmm->num_states + j]) ?
983 (float)log( ehmm->transP[i * ehmm->num_states + j] * inv_total) : -BIG_FLT ;
991 /*F///////////////////////////////////////////////////////////////////////////////////////
993 // Purpose: The function implements the mixture segmentation of the states of the
995 // Context: used with the Viterbi training of the embedded HMM
1005 static CvStatus CV_STDCALL
1006 icvMixSegmL2( CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm )
1010 CvEHMMState* state = hmm->u.ehmm[0].u.state;
1013 for (k = 0; k < num_img; k++)
1016 CvImgObsInfo* info = obs_info_array[k];
1018 for (i = 0; i < info->obs_y; i++)
1020 for (j = 0; j < info->obs_x; j++, counter++)
1022 int e_state = info->state[2 * counter + 1];
1025 min_dist = icvSquareDistance((info->obs) + (counter * info->obs_size),
1026 state[e_state].mu, info->obs_size);
1027 info->mix[counter] = 0;
1029 for (m = 1; m < state[e_state].num_mix; m++)
1031 float dist=icvSquareDistance( (info->obs) + (counter * info->obs_size),
1032 state[e_state].mu + m * info->obs_size,
1034 if (dist < min_dist)
1037 /* assign mixture with smallest distance */
1038 info->mix[counter] = m;
1048 CvStatus icvMixSegmProb(CvImgObsInfo* obs_info, int num_img, CvEHMM* hmm )
1052 CvEHMMState* state = hmm->ehmm[0].state_info;
1055 for (k = 0; k < num_img; k++)
1058 CvImgObsInfo* info = obs_info + k;
1060 for (i = 0; i < info->obs_y; i++)
1062 for (j = 0; j < info->obs_x; j++, counter++)
1064 int e_state = info->in_state[counter];
1067 max_prob = icvComputeUniModeGauss( info->obs[counter], state[e_state].mu[0],
1068 state[e_state].inv_var[0],
1069 state[e_state].log_var[0],
1071 info->mix[counter] = 0;
1073 for (m = 1; m < state[e_state].num_mix; m++)
1075 float prob=icvComputeUniModeGauss(info->obs[counter], state[e_state].mu[m],
1076 state[e_state].inv_var[m],
1077 state[e_state].log_var[m],
1079 if (prob > max_prob)
1082 // assign mixture with greatest probability.
1083 info->mix[counter] = m;
1093 static CvStatus CV_STDCALL
1094 icvViterbiSegmentation( int num_states, int /*num_obs*/, CvMatr32f transP,
1095 CvMatr32f B, int start_obs, int prob_type,
1096 int** q, int min_num_obs, int max_num_obs,
1099 // memory allocation
1101 int m_HMMType = _CV_ERGODIC; /* _CV_CAUSAL or _CV_ERGODIC */
1103 int m_ProbType = prob_type; /* _CV_LAST_STATE or _CV_BEST_STATE */
1105 int m_minNumObs = min_num_obs; /*??*/
1106 int m_maxNumObs = max_num_obs; /*??*/
1108 int m_numStates = num_states;
1110 float* m_pi = (float*)cvAlloc( num_states* sizeof(float) );
1111 CvMatr32f m_a = transP;
1113 // offset brobability matrix to starting observation
1114 CvMatr32f m_b = B + start_obs * num_states;
1115 //so m_xl will not be used more
1119 /* if (muDur != NULL){
1120 m_d = new int[m_numStates];
1121 m_l = new double[m_numStates];
1122 for (i = 0; i < m_numStates; i++){
1132 CvMatr32f m_Gamma = icvCreateMatrix_32f( num_states, m_maxNumObs );
1133 int* m_csi = (int*)cvAlloc( num_states * m_maxNumObs * sizeof(int) );
1135 //stores maximal result for every ending observation */
1136 CvVect32f m_MaxGamma = prob;
1139 // assert( m_xl + max_num_obs <= num_obs );
1141 /*??m_q = new int*[m_maxNumObs - m_minNumObs];
1142 ??for (i = 0; i < m_maxNumObs - m_minNumObs; i++)
1143 ?? m_q[i] = new int[m_minNumObs + i + 1];
1146 /******************************************************************/
1147 /* Viterbi initialization */
1148 /* set initial state probabilities, in logarithmic scale */
1149 for (i = 0; i < m_numStates; i++)
1155 for (i = 0; i < num_states; i++)
1157 m_Gamma[0 * num_states + i] = m_pi[i] + m_b[0 * num_states + i];
1158 m_csi[0 * num_states + i] = 0;
1161 /******************************************************************/
1162 /* Viterbi recursion */
1164 if ( m_HMMType == _CV_CAUSAL ) //causal model
1168 for (t = 1 ; t < m_maxNumObs; t++)
1170 // evaluate self-to-self transition for state 0
1171 m_Gamma[t * num_states + 0] = m_Gamma[(t-1) * num_states + 0] + m_a[0];
1172 m_csi[t * num_states + 0] = 0;
1174 for (j = 1; j < num_states; j++)
1176 float self = m_Gamma[ (t-1) * num_states + j] + m_a[ j * num_states + j];
1177 float prev = m_Gamma[ (t-1) * num_states +(j-1)] + m_a[ (j-1) * num_states + j];
1181 m_csi[t * num_states + j] = j-1;
1182 m_Gamma[t * num_states + j] = prev;
1186 m_csi[t * num_states + j] = j;
1187 m_Gamma[t * num_states + j] = self;
1190 m_Gamma[t * num_states + j] = m_Gamma[t * num_states + j] + m_b[t * num_states + j];
1194 else if ( m_HMMType == _CV_ERGODIC ) //ergodic model
1197 for (t = 1 ; t < m_maxNumObs; t++)
1199 for (j = 0; j < num_states; j++)
1202 m_Gamma[ t*num_states + j] = m_Gamma[(t-1) * num_states + 0] + m_a[0*num_states+j];
1203 m_csi[t *num_states + j] = 0;
1205 for (i = 1; i < num_states; i++)
1207 float currGamma = m_Gamma[(t-1) *num_states + i] + m_a[i *num_states + j];
1208 if (currGamma > m_Gamma[t *num_states + j])
1210 m_Gamma[t * num_states + j] = currGamma;
1211 m_csi[t * num_states + j] = i;
1214 m_Gamma[t *num_states + j] = m_Gamma[t *num_states + j] + m_b[t * num_states + j];
1219 for( last_obs = m_minNumObs-1, i = 0; last_obs < m_maxNumObs; last_obs++, i++ )
1223 /******************************************************************/
1224 /* Viterbi termination */
1226 if ( m_ProbType == _CV_LAST_STATE )
1228 m_MaxGamma[i] = m_Gamma[last_obs * num_states + num_states - 1];
1229 q[i][last_obs] = num_states - 1;
1231 else if( m_ProbType == _CV_BEST_STATE )
1235 m_MaxGamma[i] = m_Gamma[last_obs * num_states + 0];
1237 for(k = 1; k < num_states; k++)
1239 if ( m_Gamma[last_obs * num_states + k] > m_MaxGamma[i] )
1241 m_MaxGamma[i] = m_Gamma[last_obs * num_states + k];
1247 /******************************************************************/
1248 /* Viterbi backtracking */
1249 for (t = last_obs-1; t >= 0; t--)
1251 q[i][t] = m_csi[(t+1) * num_states + q[i][t+1] ];
1258 icvDeleteMatrix( m_Gamma );
1263 /*F///////////////////////////////////////////////////////////////////////////////////////
1264 // Name: icvEViterbi
1265 // Purpose: The function calculates the embedded Viterbi algorithm
1269 // obs_info - observations
1272 // Returns: the Embedded Viterbi probability (float)
1273 // and do state segmentation of observations
1277 static float CV_STDCALL icvEViterbi( CvImgObsInfo* obs_info, CvEHMM* hmm )
1280 float log_likelihood;
1282 float inv_obs_x = 1.f / obs_info->obs_x;
1284 CvEHMMState* first_state = hmm->u.ehmm->u.state;
1286 /* memory allocation for superB */
1287 CvMatr32f superB = icvCreateMatrix_32f(hmm->num_states, obs_info->obs_y );
1289 /* memory allocation for q */
1290 int*** q = (int***)cvAlloc( hmm->num_states * sizeof(int**) );
1291 int* super_q = (int*)cvAlloc( obs_info->obs_y * sizeof(int) );
1293 for (i = 0; i < hmm->num_states; i++)
1295 q[i] = (int**)cvAlloc( obs_info->obs_y * sizeof(int*) );
1297 for (j = 0; j < obs_info->obs_y ; j++)
1299 q[i][j] = (int*)cvAlloc( obs_info->obs_x * sizeof(int) );
1303 /* start Viterbi segmentation */
1304 for (i = 0; i < hmm->num_states; i++)
1306 CvEHMM* ehmm = &(hmm->u.ehmm[i]);
1308 for (j = 0; j < obs_info->obs_y; j++)
1312 /* 1D HMM Viterbi segmentation */
1313 icvViterbiSegmentation( ehmm->num_states, obs_info->obs_x,
1314 ehmm->transP, ehmm->obsProb[j], 0,
1315 _CV_LAST_STATE, &q[i][j], obs_info->obs_x,
1316 obs_info->obs_x, &max_gamma);
1318 superB[j * hmm->num_states + i] = max_gamma * inv_obs_x;
1322 /* perform global Viterbi segmentation (i.e. process higher-level HMM) */
1324 icvViterbiSegmentation( hmm->num_states, obs_info->obs_y,
1325 hmm->transP, superB, 0,
1326 _CV_LAST_STATE, &super_q, obs_info->obs_y,
1327 obs_info->obs_y, &log_likelihood );
1329 log_likelihood /= obs_info->obs_y ;
1333 /* assign new state to observation vectors */
1334 for (i = 0; i < obs_info->obs_y; i++)
1336 for (j = 0; j < obs_info->obs_x; j++, counter++)
1338 int superstate = super_q[i];
1339 int state = (int)(hmm->u.ehmm[superstate].u.state - first_state);
1341 obs_info->state[2 * counter] = superstate;
1342 obs_info->state[2 * counter + 1] = state + q[superstate][i][j];
1346 /* memory deallocation for superB */
1347 icvDeleteMatrix( superB );
1349 /*memory deallocation for q */
1350 for (i = 0; i < hmm->num_states; i++)
1352 for (j = 0; j < obs_info->obs_y ; j++)
1362 return log_likelihood;
1365 static CvStatus CV_STDCALL
1366 icvEstimateHMMStateParams( CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm )
1368 /* compute gamma, weights, means, vars */
1371 int vect_len = obs_info_array[0]->obs_size;
1373 float start_log_var_val = LN2PI * vect_len;
1375 CvVect32f tmp_vect = icvCreateVector_32f( vect_len );
1377 CvEHMMState* first_state = hmm->u.ehmm[0].u.state;
1379 assert( sizeof(float) == sizeof(int) );
1381 for(i = 0; i < hmm->num_states; i++ )
1383 total+= hmm->u.ehmm[i].num_states;
1386 /***************Gamma***********************/
1387 /* initialize gamma */
1388 for( i = 0; i < total; i++ )
1390 for (m = 0; m < first_state[i].num_mix; m++)
1392 ((int*)(first_state[i].weight))[m] = 0;
1396 /* maybe gamma must be computed in mixsegm process ?? */
1399 for (k = 0; k < num_img; k++)
1401 CvImgObsInfo* info = obs_info_array[k];
1402 int num_obs = info->obs_y * info->obs_x;
1404 for (i = 0; i < num_obs; i++)
1407 state = info->state[2*i + 1];
1408 mixture = info->mix[i];
1409 /* computes gamma - number of observations corresponding
1410 to every mixture of every state */
1411 ((int*)(first_state[state].weight))[mixture] += 1;
1414 /***************Mean and Var***********************/
1415 /* compute means and variances of every item */
1416 /* initially variance placed to inv_var */
1417 /* zero mean and variance */
1418 for (i = 0; i < total; i++)
1420 memset( (void*)first_state[i].mu, 0, first_state[i].num_mix * vect_len *
1422 memset( (void*)first_state[i].inv_var, 0, first_state[i].num_mix * vect_len *
1427 for (i = 0; i < num_img; i++)
1429 CvImgObsInfo* info = obs_info_array[i];
1430 int total_obs = info->obs_x * info->obs_y;
1432 float* vector = info->obs;
1434 for (j = 0; j < total_obs; j++, vector+=vect_len )
1436 int state = info->state[2 * j + 1];
1437 int mixture = info->mix[j];
1439 CvVect32f mean = first_state[state].mu + mixture * vect_len;
1440 CvVect32f mean2 = first_state[state].inv_var + mixture * vect_len;
1442 icvAddVector_32f( mean, vector, mean, vect_len );
1443 for( k = 0; k < vect_len; k++ )
1444 mean2[k] += vector[k]*vector[k];
1448 /*compute the means and variances */
1449 /* assume gamma already computed */
1450 for (i = 0; i < total; i++)
1452 CvEHMMState* state = &(first_state[i]);
1454 for (m = 0; m < state->num_mix; m++)
1457 CvVect32f mu = state->mu + m * vect_len;
1458 CvVect32f invar = state->inv_var + m * vect_len;
1460 if ( ((int*)state->weight)[m] > 1)
1462 float inv_gamma = 1.f/((int*)(state->weight))[m];
1464 icvScaleVector_32f( mu, mu, vect_len, inv_gamma);
1465 icvScaleVector_32f( invar, invar, vect_len, inv_gamma);
1468 icvMulVectors_32f(mu, mu, tmp_vect, vect_len);
1469 icvSubVector_32f( invar, tmp_vect, invar, vect_len);
1471 /* low bound of variance - 100 (Ara's experimental result) */
1472 for( k = 0; k < vect_len; k++ )
1474 invar[k] = (invar[k] > 100.f) ? invar[k] : 100.f;
1477 /* compute log_var */
1478 state->log_var_val[m] = start_log_var_val;
1479 for( k = 0; k < vect_len; k++ )
1481 state->log_var_val[m] += (float)log( invar[k] );
1484 /* SMOLI 27.10.2000 */
1485 state->log_var_val[m] *= 0.5;
1488 /* compute inv_var = 1/sqrt(2*variance) */
1489 icvScaleVector_32f(invar, invar, vect_len, 2.f );
1490 cvbInvSqrt( invar, invar, vect_len );
1494 /***************Weights***********************/
1495 /* normilize gammas - i.e. compute mixture weights */
1498 for (i = 0; i < total; i++)
1500 int gamma_total = 0;
1503 for (m = 0; m < first_state[i].num_mix; m++)
1505 gamma_total += ((int*)(first_state[i].weight))[m];
1508 norm = gamma_total ? (1.f/(float)gamma_total) : 0.f;
1510 for (m = 0; m < first_state[i].num_mix; m++)
1512 first_state[i].weight[m] = ((int*)(first_state[i].weight))[m] * norm;
1516 icvDeleteVector( tmp_vect);
1521 CvStatus icvLightingCorrection8uC1R( uchar* img, CvSize roi, int src_step )
1524 int width = roi.width;
1525 int height = roi.height;
1527 float x1, x2, y1, y2;
1528 int f[3] = {0, 0, 0};
1529 float a[3] = {0, 0, 0};
1536 float min = FLT_MAX;
1537 float max = -FLT_MAX;
1540 float* float_img = icvAlloc( width * height * sizeof(float) );
1542 x1 = width * (width + 1) / 2.0f; // Sum (1, ... , width)
1543 x2 = width * (width + 1 ) * (2 * width + 1) / 6.0f; // Sum (1^2, ... , width^2)
1544 y1 = height * (height + 1)/2.0f; // Sum (1, ... , width)
1545 y2 = height * (height + 1 ) * (2 * height + 1) / 6.0f; // Sum (1^2, ... , width^2)
1548 // extract grayvalues
1549 for (i = 0; i < height; i++)
1551 for (j = 0; j < width; j++)
1553 f[2] = f[2] + j * img[i*src_step + j];
1554 f[1] = f[1] + i * img[i*src_step + j];
1555 f[0] = f[0] + img[i*src_step + j];
1559 h1 = (float)f[0] * (float)x1 / (float)width;
1560 h2 = (float)f[0] * (float)y1 / (float)height;
1562 a[2] = ((float)f[2] - h1) / (float)(x2*height - x1*x1*height/(float)width);
1563 a[1] = ((float)f[1] - h2) / (float)(y2*width - y1*y1*width/(float)height);
1564 a[0] = (float)f[0]/(float)(width*height) - (float)y1*a[1]/(float)height -
1565 (float)x1*a[2]/(float)width;
1567 for (i = 0; i < height; i++)
1569 for (j = 0; j < width; j++)
1572 correction = a[0] + a[1]*(float)i + a[2]*(float)j;
1574 float_img[i*width + j] = img[i*src_step + j] - correction;
1576 if (float_img[i*width + j] < min) min = float_img[i*width+j];
1577 if (float_img[i*width + j] > max) max = float_img[i*width+j];
1581 //rescaling to the range 0:255
1586 c2 = 255.0f/(float)(max - min);
1588 c1 = (-(float)min)*c2;
1590 for (i = 0; i < height; i++)
1592 for (j = 0; j < width; j++)
1594 int value = (int)floor(c2*float_img[i*width + j] + c1);
1595 if (value < 0) value = 0;
1596 if (value > 255) value = 255;
1597 img[i*src_step + j] = (uchar)value;
1601 cvFree( &float_img );
1606 CvStatus icvLightingCorrection( icvImage* img )
1609 if ( img->type != IPL_DEPTH_8U || img->channels != 1 )
1610 return CV_BADFACTOR_ERR;
1612 roi = _cvSize( img->roi.width, img->roi.height );
1614 return _cvLightingCorrection8uC1R( img->data + img->roi.y * img->step + img->roi.x,
1622 cvCreate2DHMM( int *state_number, int *num_mix, int obs_size )
1626 CV_FUNCNAME( "cvCreate2DHMM" );
1630 IPPI_CALL( icvCreate2DHMM( &hmm, state_number, num_mix, obs_size ));
1638 cvRelease2DHMM( CvEHMM ** hmm )
1640 CV_FUNCNAME( "cvRelease2DHMM" );
1644 IPPI_CALL( icvRelease2DHMM( hmm ));
1648 CV_IMPL CvImgObsInfo*
1649 cvCreateObsInfo( CvSize num_obs, int obs_size )
1651 CvImgObsInfo *obs_info = 0;
1653 CV_FUNCNAME( "cvCreateObsInfo" );
1657 IPPI_CALL( icvCreateObsInfo( &obs_info, num_obs, obs_size ));
1665 cvReleaseObsInfo( CvImgObsInfo ** obs_info )
1667 CV_FUNCNAME( "cvReleaseObsInfo" );
1671 IPPI_CALL( icvReleaseObsInfo( obs_info ));
1678 cvUniformImgSegm( CvImgObsInfo * obs_info, CvEHMM * hmm )
1680 CV_FUNCNAME( "cvUniformImgSegm" );
1684 IPPI_CALL( icvUniformImgSegm( obs_info, hmm ));
1690 cvInitMixSegm( CvImgObsInfo ** obs_info_array, int num_img, CvEHMM * hmm )
1692 CV_FUNCNAME( "cvInitMixSegm" );
1696 IPPI_CALL( icvInitMixSegm( obs_info_array, num_img, hmm ));
1702 cvEstimateHMMStateParams( CvImgObsInfo ** obs_info_array, int num_img, CvEHMM * hmm )
1704 CV_FUNCNAME( "cvEstimateHMMStateParams" );
1708 IPPI_CALL( icvEstimateHMMStateParams( obs_info_array, num_img, hmm ));
1714 cvEstimateTransProb( CvImgObsInfo ** obs_info_array, int num_img, CvEHMM * hmm )
1716 CV_FUNCNAME( "cvEstimateTransProb" );
1720 IPPI_CALL( icvEstimateTransProb( obs_info_array, num_img, hmm ));
1727 cvEstimateObsProb( CvImgObsInfo * obs_info, CvEHMM * hmm )
1729 CV_FUNCNAME( "cvEstimateObsProb" );
1733 IPPI_CALL( icvEstimateObsProb( obs_info, hmm ));
1739 cvEViterbi( CvImgObsInfo * obs_info, CvEHMM * hmm )
1741 float result = FLT_MAX;
1743 CV_FUNCNAME( "cvEViterbi" );
1747 if( (obs_info == NULL) || (hmm == NULL) )
1748 CV_ERROR( CV_BadDataPtr, "Null pointer." );
1750 result = icvEViterbi( obs_info, hmm );
1758 cvMixSegmL2( CvImgObsInfo ** obs_info_array, int num_img, CvEHMM * hmm )
1760 CV_FUNCNAME( "cvMixSegmL2" );
1764 IPPI_CALL( icvMixSegmL2( obs_info_array, num_img, hmm ));