1 /*M/////////////////////////////////////////////////////////////////////////////////////// 2 // 3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. 4 // 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. 8 // 9 // 10 // Intel License Agreement 11 // For Open Source Computer Vision Library 12 // 13 // Copyright (C) 2000, Intel Corporation, all rights reserved. 14 // Third party copyrights are property of their respective owners. 15 // 16 // Redistribution and use in source and binary forms, with or without modification, 17 // are permitted provided that the following conditions are met: 18 // 19 // * Redistribution's of source code must retain the above copyright notice, 20 // this list of conditions and the following disclaimer. 21 // 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. 25 // 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. 28 // 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. 39 // 40 //M*/ 41 42 #include "_cvaux.h" 43 44 #define LN2PI 1.837877f 45 #define BIG_FLT 1.e+10f 46 47 48 #define _CV_ERGODIC 1 49 #define _CV_CAUSAL 2 50 51 #define _CV_LAST_STATE 1 52 #define _CV_BEST_STATE 2 53 54 55 //*F/////////////////////////////////////////////////////////////////////////////////////// 56 // Name: _cvCreateObsInfo 57 // Purpose: The function allocates memory for CvImgObsInfo structure 58 // and its inner stuff 59 // Context: 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 64 // 65 // Returns: error status 66 // 67 // Notes: 68 //F*/ 69 static CvStatus CV_STDCALL icvCreateObsInfo( CvImgObsInfo** obs_info, 70 CvSize num_obs, int obs_size ) 71 { 72 int total = num_obs.height * num_obs.width; 73 74 CvImgObsInfo* obs = (CvImgObsInfo*)cvAlloc( sizeof( CvImgObsInfo) ); 75 76 obs->obs_x = num_obs.width; 77 obs->obs_y = num_obs.height; 78 79 obs->obs = (float*)cvAlloc( total * obs_size * sizeof(float) ); 80 81 obs->state = (int*)cvAlloc( 2 * total * sizeof(int) ); 82 obs->mix = (int*)cvAlloc( total * sizeof(int) ); 83 84 obs->obs_size = obs_size; 85 86 obs_info[0] = obs; 87 88 return CV_NO_ERR; 89 } 90 91 static CvStatus CV_STDCALL icvReleaseObsInfo( CvImgObsInfo** p_obs_info ) 92 { 93 CvImgObsInfo* obs_info = p_obs_info[0]; 94 95 cvFree( &(obs_info->obs) ); 96 cvFree( &(obs_info->mix) ); 97 cvFree( &(obs_info->state) ); 98 cvFree( &(obs_info) ); 99 100 p_obs_info[0] = NULL; 101 102 return CV_NO_ERR; 103 } 104 105 106 //*F/////////////////////////////////////////////////////////////////////////////////////// 107 // Name: icvCreate2DHMM 108 // Purpose: The function allocates memory for 2-dimensional embedded HMM model 109 // and its inner stuff 110 // Context: 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 116 // 117 // Returns: error status 118 // 119 // Notes: state_number[0] - number of states in external HMM. 120 // state_number[i] - number of states in embedded HMM 121 // 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// 124 // 125 //F*/ 126 static CvStatus CV_STDCALL icvCreate2DHMM( CvEHMM** this_hmm, 127 int* state_number, int* num_mix, int obs_size ) 128 { 129 int i; 130 int real_states = 0; 131 132 CvEHMMState* all_states; 133 CvEHMM* hmm; 134 int total_mix = 0; 135 float* pointers; 136 137 //compute total number of states of all level in 2d EHMM 138 for( i = 1; i <= state_number[0]; i++ ) 139 { 140 real_states += state_number[i]; 141 } 142 143 /* allocate memory for all hmms (from all levels) */ 144 hmm = (CvEHMM*)cvAlloc( (state_number[0] + 1) * sizeof(CvEHMM) ); 145 146 /* set number of superstates */ 147 hmm[0].num_states = state_number[0]; 148 hmm[0].level = 1; 149 150 /* allocate memory for all states */ 151 all_states = (CvEHMMState *)cvAlloc( real_states * sizeof( CvEHMMState ) ); 152 153 /* assign number of mixtures */ 154 for( i = 0; i < real_states; i++ ) 155 { 156 all_states[i].num_mix = num_mix[i]; 157 } 158 159 /* compute size of inner of all real states */ 160 for( i = 0; i < real_states; i++ ) 161 { 162 total_mix += num_mix[i]; 163 } 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) ); 167 168 /* organize memory */ 169 for( i = 0; i < real_states; i++ ) 170 { 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; 173 174 all_states[i].log_var_val = pointers; pointers += num_mix[i]; 175 all_states[i].weight = pointers; pointers += num_mix[i]; 176 } 177 178 /* set pointer to embedded hmm array */ 179 hmm->u.ehmm = hmm + 1; 180 181 for( i = 0; i < hmm[0].num_states; i++ ) 182 { 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]; 186 } 187 188 for( i = 0; i <= state_number[0]; i++ ) 189 { 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; 193 } 194 195 /* if all ok - return pointer */ 196 *this_hmm = hmm; 197 return CV_NO_ERR; 198 } 199 200 static CvStatus CV_STDCALL icvRelease2DHMM( CvEHMM** phmm ) 201 { 202 CvEHMM* hmm = phmm[0]; 203 int i; 204 for( i = 0; i < hmm[0].num_states + 1; i++ ) 205 { 206 icvDeleteMatrix( hmm[i].transP ); 207 } 208 209 if (hmm->obsProb != NULL) 210 { 211 int* tmp = ((int*)(hmm->obsProb)) - 3; 212 cvFree( &(tmp) ); 213 } 214 215 cvFree( &(hmm->u.ehmm->u.state->mu) ); 216 cvFree( &(hmm->u.ehmm->u.state) ); 217 218 219 /* free hmm structures */ 220 cvFree( phmm ); 221 222 phmm[0] = NULL; 223 224 return CV_NO_ERR; 225 } 226 227 /* distance between 2 vectors */ 228 static float icvSquareDistance( CvVect32f v1, CvVect32f v2, int len ) 229 { 230 int i; 231 double dist0 = 0; 232 double dist1 = 0; 233 234 for( i = 0; i <= len - 4; i += 4 ) 235 { 236 double t0 = v1[i] - v2[i]; 237 double t1 = v1[i+1] - v2[i+1]; 238 dist0 += t0*t0; 239 dist1 += t1*t1; 240 241 t0 = v1[i+2] - v2[i+2]; 242 t1 = v1[i+3] - v2[i+3]; 243 dist0 += t0*t0; 244 dist1 += t1*t1; 245 } 246 247 for( ; i < len; i++ ) 248 { 249 double t0 = v1[i] - v2[i]; 250 dist0 += t0*t0; 251 } 252 253 return (float)(dist0 + dist1); 254 } 255 256 /*can be used in CHMM & DHMM */ 257 static CvStatus CV_STDCALL 258 icvUniformImgSegm( CvImgObsInfo* obs_info, CvEHMM* hmm ) 259 { 260 #if 1 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; 266 267 /* check arguments */ 268 if ( !obs_info || !hmm ) return CV_NULLPTR_ERR; 269 270 first_state = hmm->u.ehmm->u.state; 271 272 for (i = 0; i < obs_info->obs_y; i++) 273 { 274 //bad line (division ) 275 int superstate = (int)((i * hmm->num_states)*inv_y);/* /obs_info->obs_y; */ 276 277 int index = (int)(hmm->u.ehmm[superstate].u.state - first_state); 278 279 for (j = 0; j < obs_info->obs_x; j++, counter++) 280 { 281 int state = (int)((j * hmm->u.ehmm[superstate].num_states)* inv_x); /* / obs_info->obs_x; */ 282 283 obs_info->state[2 * counter] = superstate; 284 obs_info->state[2 * counter + 1] = state + index; 285 } 286 } 287 #else 288 //this is not ready yet 289 290 int i,j,k,m; 291 CvEHMMState* first_state = hmm->u.ehmm->u.state; 292 293 /* check bad arguments */ 294 if ( hmm->num_states > obs_info->obs_y ) return CV_BADSIZE_ERR; 295 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 */ 299 300 //for every horizontal band compute subdivision 301 for( i = 0; i < hmm->num_states; i++ ) 302 { 303 CvEHMM* ehmm = &(hmm->u.ehmm[i]); 304 col_per_state[i] = (float)obs_info->obs_x / ehmm->num_states; 305 } 306 307 //compute state bounds 308 int ss_bound[1024]; 309 for( i = 0; i < hmm->num_states - 1; i++ ) 310 { 311 ss_bound[i] = floor( row_per_state * ( i+1 ) ); 312 } 313 ss_bound[hmm->num_states - 1] = obs_info->obs_y; 314 315 //work inside every superstate 316 317 int row = 0; 318 319 for( i = 0; i < hmm->num_states; i++ ) 320 { 321 CvEHMM* ehmm = &(hmm->u.ehmm[i]); 322 int index = ehmm->u.state - first_state; 323 324 //calc distribution in superstate 325 int es_bound[1024]; 326 for( j = 0; j < ehmm->num_states - 1; j++ ) 327 { 328 es_bound[j] = floor( col_per_state[i] * ( j+1 ) ); 329 } 330 es_bound[ehmm->num_states - 1] = obs_info->obs_x; 331 332 //assign states to first row of superstate 333 int col = 0; 334 for( j = 0; j < ehmm->num_states; j++ ) 335 { 336 for( k = col; k < es_bound[j]; k++, col++ ) 337 { 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; 340 } 341 col = es_bound[j]; 342 } 343 344 //copy the same to other rows of superstate 345 for( m = row; m < ss_bound[i]; m++ ) 346 { 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) ); 349 } 350 351 row = ss_bound[i]; 352 } 353 354 #endif 355 356 return CV_NO_ERR; 357 } 358 359 360 /*F/////////////////////////////////////////////////////////////////////////////////////// 361 // Name: InitMixSegm 362 // Purpose: The function implements the mixture segmentation of the states of the 363 // embedded HMM 364 // Context: used with the Viterbi training of the embedded HMM 365 // Function uses K-Means algorithm for clustering 366 // 367 // Parameters: obs_info_array - array of pointers to image observations 368 // num_img - length of above array 369 // hmm - pointer to HMM structure 370 // 371 // Returns: error status 372 // 373 // Notes: 374 //F*/ 375 static CvStatus CV_STDCALL 376 icvInitMixSegm( CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm ) 377 { 378 int k, i, j; 379 int* num_samples; /* number of observations in every state */ 380 int* counter; /* array of counters for every state */ 381 382 int** a_class; /* for every state - characteristic array */ 383 384 CvVect32f** samples; /* for every state - pointer to observation vectors */ 385 int*** samples_mix; /* for every state - array of pointers to vectors mixtures */ 386 387 CvTermCriteria criteria = cvTermCriteria( CV_TERMCRIT_EPS|CV_TERMCRIT_ITER, 388 1000, /* iter */ 389 0.01f ); /* eps */ 390 391 int total = 0; 392 393 CvEHMMState* first_state = hmm->u.ehmm->u.state; 394 395 for( i = 0 ; i < hmm->num_states; i++ ) 396 { 397 total += hmm->u.ehmm[i].num_states; 398 } 399 400 /* for every state integer is allocated - number of vectors in state */ 401 num_samples = (int*)cvAlloc( total * sizeof(int) ); 402 403 /* integer counter is allocated for every state */ 404 counter = (int*)cvAlloc( total * sizeof(int) ); 405 406 samples = (CvVect32f**)cvAlloc( total * sizeof(CvVect32f*) ); 407 samples_mix = (int***)cvAlloc( total * sizeof(int**) ); 408 409 /* clear */ 410 memset( num_samples, 0 , total*sizeof(int) ); 411 memset( counter, 0 , total*sizeof(int) ); 412 413 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++) 416 { 417 CvImgObsInfo* obs = obs_info_array[k]; 418 int count = 0; 419 420 for (i = 0; i < obs->obs_y; i++) 421 { 422 for (j = 0; j < obs->obs_x; j++, count++) 423 { 424 int state = obs->state[ 2 * count + 1]; 425 num_samples[state] += 1; 426 } 427 } 428 } 429 430 /* for every state int* is allocated */ 431 a_class = (int**)cvAlloc( total*sizeof(int*) ); 432 433 for (i = 0; i < total; i++) 434 { 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*) ); 438 } 439 440 /* for every state vectors which belong to state are gathered */ 441 for (k = 0; k < num_img; k++) 442 { 443 CvImgObsInfo* obs = obs_info_array[k]; 444 int num_obs = ( obs->obs_x ) * ( obs->obs_y ); 445 float* vector = obs->obs; 446 447 for (i = 0; i < num_obs; i++, vector+=obs->obs_size ) 448 { 449 int state = obs->state[2*i+1]; 450 451 samples[state][counter[state]] = vector; 452 samples_mix[state][counter[state]] = &(obs->mix[i]); 453 counter[state]++; 454 } 455 } 456 457 /* clear counters */ 458 memset( counter, 0, total*sizeof(int) ); 459 460 /* do the actual clustering using the K Means algorithm */ 461 for (i = 0; i < total; i++) 462 { 463 if ( first_state[i].num_mix == 1) 464 { 465 for (k = 0; k < num_samples[i]; k++) 466 { 467 /* all vectors belong to one mixture */ 468 a_class[i][k] = 0; 469 } 470 } 471 else if( num_samples[i] ) 472 { 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] ); 476 } 477 } 478 479 /* for every vector number of mixture is assigned */ 480 for( i = 0; i < total; i++ ) 481 { 482 for (j = 0; j < num_samples[i]; j++) 483 { 484 samples_mix[i][j][0] = a_class[i][j]; 485 } 486 } 487 488 for (i = 0; i < total; i++) 489 { 490 cvFree( &(a_class[i]) ); 491 cvFree( &(samples[i]) ); 492 cvFree( &(samples_mix[i]) ); 493 } 494 495 cvFree( &a_class ); 496 cvFree( &samples ); 497 cvFree( &samples_mix ); 498 cvFree( &counter ); 499 cvFree( &num_samples ); 500 501 return CV_NO_ERR; 502 } 503 504 /*F/////////////////////////////////////////////////////////////////////////////////////// 505 // Name: ComputeUniModeGauss 506 // Purpose: The function computes the Gaussian pdf for a sample vector 507 // Context: 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 512 // 513 // Returns: the pdf of the sample vector given the specified Gaussian 514 // 515 // Notes: 516 //F*/ 517 /*static float icvComputeUniModeGauss(CvVect32f vect, CvVect32f mu, 518 CvVect32f inv_var, float log_var_val, int vect_size) 519 { 520 int n; 521 double tmp; 522 double prob; 523 524 prob = -log_var_val; 525 526 for (n = 0; n < vect_size; n++) 527 { 528 tmp = (vect[n] - mu[n]) * inv_var[n]; 529 prob = prob - tmp * tmp; 530 } 531 //prob *= 0.5f; 532 533 return (float)prob; 534 }*/ 535 536 /*F/////////////////////////////////////////////////////////////////////////////////////// 537 // Name: ComputeGaussMixture 538 // Purpose: The function computes the mixture Gaussian pdf of a sample vector. 539 // Context: 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 550 // 551 // Returns: the pdf of the sample vector given the specified Gaussian mixture. 552 // 553 // Notes: 554 //F*/ 555 /* Calculate probability of observation at state in logarithmic scale*/ 556 /*static float 557 icvComputeGaussMixture( CvVect32f vect, float* mu, 558 float* inv_var, float* log_var_val, 559 int vect_size, float* weight, int num_mix ) 560 { 561 double prob, l_prob; 562 563 prob = 0.0f; 564 565 if (num_mix == 1) 566 { 567 return icvComputeUniModeGauss( vect, mu, inv_var, log_var_val[0], vect_size); 568 } 569 else 570 { 571 int m; 572 for (m = 0; m < num_mix; m++) 573 { 574 if ( weight[m] > 0.0) 575 { 576 l_prob = icvComputeUniModeGauss(vect, mu + m*vect_size, 577 inv_var + m * vect_size, 578 log_var_val[m], 579 vect_size); 580 581 prob = prob + weight[m]*exp((double)l_prob); 582 } 583 } 584 prob = log(prob); 585 } 586 return (float)prob; 587 }*/ 588 589 590 /*F/////////////////////////////////////////////////////////////////////////////////////// 591 // Name: EstimateObsProb 592 // Purpose: The function computes the probability of every observation in every state 593 // Context: 594 // Parameters: obs_info - observations 595 // hmm - hmm 596 // Returns: error status 597 // 598 // Notes: 599 //F*/ 600 static CvStatus CV_STDCALL icvEstimateObsProb( CvImgObsInfo* obs_info, CvEHMM* hmm ) 601 { 602 int i, j; 603 int total_states = 0; 604 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 */ 609 610 for( j = 0; j < hmm->num_states; j++ ) 611 { 612 total_states += hmm->u.ehmm[j].num_states; 613 } 614 615 if ( hmm->obsProb == NULL ) 616 { 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) ); 620 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); 626 status = 3; 627 628 } 629 else 630 { 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*)*/ ) ); 635 636 assert( sizeof(float*) == sizeof(int) ); 637 638 if ( need_size > (*total) ) 639 { 640 int* buffer = ((int*)(hmm->obsProb)) - 3; 641 cvFree( &buffer); 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; 646 647 hmm->obsProb = (float**)(buffer + 3); 648 649 status = 3; 650 } 651 } 652 if (!status) 653 { 654 int* obsx = ((int*)(hmm->obsProb)) - 1; 655 int* obsy = ((int*)(hmm->obsProb)) - 2; 656 657 assert( (*obsx > 0) && (*obsy > 0) ); 658 659 /* is good distribution? */ 660 if ( (obs_info->obs_x > (*obsx) ) || (obs_info->obs_y > (*obsy) ) ) 661 status = 3; 662 } 663 664 /* if bad status - do reallocation actions */ 665 assert( (status == 0) || (status == 3) ); 666 667 if ( status ) 668 { 669 float** tmp = hmm->obsProb; 670 float* tmpf; 671 672 /* distribute pointers of ehmm->obsProb */ 673 for( i = 0; i < hmm->num_states; i++ ) 674 { 675 hmm->u.ehmm[i].obsProb = tmp; 676 tmp += obs_info->obs_y; 677 } 678 679 tmpf = (float*)tmp; 680 681 /* distribute pointers of ehmm->obsProb[j] */ 682 for( i = 0; i < hmm->num_states; i++ ) 683 { 684 CvEHMM* ehmm = &( hmm->u.ehmm[i] ); 685 686 for( j = 0; j < obs_info->obs_y; j++ ) 687 { 688 ehmm->obsProb[j] = tmpf; 689 tmpf += ehmm->num_states * obs_info->obs_x; 690 } 691 } 692 }/* end of pointer distribution */ 693 694 #if 1 695 { 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; 701 702 float* log_mix_prob = local_log_mix_prob; 703 double* mix_prob = local_mix_prob; 704 705 int max_size = 0; 706 int obs_x = obs_info->obs_x; 707 708 /* calculate temporary buffer size */ 709 for( i = 0; i < hmm->num_states; i++ ) 710 { 711 CvEHMM* ehmm = &(hmm->u.ehmm[i]); 712 CvEHMMState* state = ehmm->u.state; 713 714 int max_mix = 0; 715 for( j = 0; j < ehmm->num_states; j++ ) 716 { 717 int t = state[j].num_mix; 718 if( max_mix < t ) max_mix = t; 719 } 720 max_mix *= ehmm->num_states; 721 if( max_size < max_mix ) max_size = max_mix; 722 } 723 724 max_size *= obs_x * vect_size; 725 726 /* allocate buffer */ 727 if( max_size > MAX_BUF_SIZE ) 728 { 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); 732 } 733 734 memset( log_mix_prob, 0, max_size*sizeof(float)); 735 736 /*****************computing probabilities***********************/ 737 738 /* loop through external states */ 739 for( i = 0; i < hmm->num_states; i++ ) 740 { 741 CvEHMM* ehmm = &(hmm->u.ehmm[i]); 742 CvEHMMState* state = ehmm->u.state; 743 744 int max_mix = 0; 745 int n_states = ehmm->num_states; 746 747 /* determine maximal number of mixtures (again) */ 748 for( j = 0; j < ehmm->num_states; j++ ) 749 { 750 int t = state[j].num_mix; 751 if( max_mix < t ) max_mix = t; 752 } 753 754 /* loop through rows of the observation matrix */ 755 for( j = 0; j < obs_info->obs_y; j++ ) 756 { 757 int m, n; 758 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; 762 763 /* several passes are done below */ 764 765 /* 1. calculate logarithms of probabilities for each mixture */ 766 767 /* loop through mixtures */ 768 for( m = 0; m < max_mix; m++ ) 769 { 770 /* set pointer to first observation in the line */ 771 float* vect = obs; 772 773 /* cycles through obseravtions in the line */ 774 for( n = 0; n < obs_x; n++, vect += vect_size, log_mp += n_states ) 775 { 776 int k, l; 777 for( l = 0; l < n_states; l++ ) 778 { 779 if( state[l].num_mix > m ) 780 { 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++ ) 785 { 786 double t = (vect[k] - mu[k])*inv_var[k]; 787 prob -= t*t; 788 } 789 log_mp[l] = MAX( (float)prob, -500 ); 790 } 791 } 792 } 793 } 794 795 /* skip the rest if there is a single mixture */ 796 if( max_mix == 1 ) continue; 797 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 ); 801 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 ) 805 { 806 int l; 807 for( l = 0; l < n_states; l++ ) 808 { 809 mp[l] *= state[l].weight[0]; 810 } 811 } 812 813 /* 3b. add other mixtures */ 814 for( m = 1; m < max_mix; m++ ) 815 { 816 int ofs = -m*obs_x*n_states; 817 for( n = 0; n < obs_x; n++, mp += n_states ) 818 { 819 int l; 820 for( l = 0; l < n_states; l++ ) 821 { 822 if( m < state[l].num_mix ) 823 { 824 mp[l + ofs] += mp[l] * state[l].weight[m]; 825 } 826 } 827 } 828 } 829 830 /* 4. Put logarithms of summary probabilities to the destination matrix */ 831 cvbFastLog( mix_prob, ehmm->obsProb[j], obs_x * n_states ); 832 } 833 } 834 835 if( log_mix_prob != local_log_mix_prob ) cvFree( &log_mix_prob ); 836 return res; 837 #undef MAX_BUF_SIZE 838 } 839 #else 840 for( i = 0; i < hmm->num_states; i++ ) 841 { 842 CvEHMM* ehmm = &(hmm->u.ehmm[i]); 843 CvEHMMState* state = ehmm->u.state; 844 845 for( j = 0; j < obs_info->obs_y; j++ ) 846 { 847 int k,m; 848 849 int obs_index = j * obs_info->obs_x; 850 851 float* B = ehmm->obsProb[j]; 852 853 /* cycles through obs and states */ 854 for( k = 0; k < obs_info->obs_x; k++ ) 855 { 856 CvVect32f vect = (obs_info->obs) + (obs_index + k) * vect_size; 857 858 float* matr_line = B + k * ehmm->num_states; 859 860 for( m = 0; m < ehmm->num_states; m++ ) 861 { 862 matr_line[m] = icvComputeGaussMixture( vect, state[m].mu, state[m].inv_var, 863 state[m].log_var_val, vect_size, state[m].weight, 864 state[m].num_mix ); 865 } 866 } 867 } 868 } 869 #endif 870 } 871 872 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 878 // Context: 879 // Parameters: obs_info_array - array of pointers to image observations 880 // num_img - length of above array 881 // hmm - pointer to HMM structure 882 // Returns: void 883 // 884 // Notes: 885 //F*/ 886 static CvStatus CV_STDCALL 887 icvEstimateTransProb( CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm ) 888 { 889 int i, j, k; 890 891 CvEHMMState* first_state = hmm->u.ehmm->u.state; 892 /* as a counter we will use transP matrix */ 893 894 /* initialization */ 895 896 /* clear transP */ 897 icvSetZero_32f( hmm->transP, hmm->num_states, hmm->num_states ); 898 for (i = 0; i < hmm->num_states; i++ ) 899 { 900 icvSetZero_32f( hmm->u.ehmm[i].transP , hmm->u.ehmm[i].num_states, hmm->u.ehmm[i].num_states ); 901 } 902 903 /* compute the counters */ 904 for (i = 0; i < num_img; i++) 905 { 906 int counter = 0; 907 CvImgObsInfo* info = obs_info_array[i]; 908 909 for (j = 0; j < info->obs_y; j++) 910 { 911 for (k = 0; k < info->obs_x; k++, counter++) 912 { 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; 917 int begin_ind; 918 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; 922 923 if (j < info->obs_y - 1) 924 { 925 int transP_size = hmm->num_states; 926 927 nextsuperstate = info->state[ 2*(counter + info->obs_x) ]; 928 929 hmm->transP[superstate * transP_size + nextsuperstate] += 1; 930 } 931 932 if (k < info->obs_x - 1) 933 { 934 int transP_size = hmm->u.ehmm[superstate].num_states; 935 936 nextstate = info->state[2*(counter+1) + 1] - begin_ind; 937 hmm->u.ehmm[superstate].transP[ state * transP_size + nextstate] += 1; 938 } 939 } 940 } 941 } 942 /* estimate superstate matrix */ 943 for( i = 0; i < hmm->num_states; i++) 944 { 945 float total = 0; 946 float inv_total; 947 for( j = 0; j < hmm->num_states; j++) 948 { 949 total += hmm->transP[i * hmm->num_states + j]; 950 } 951 //assert( total ); 952 953 inv_total = total ? 1.f/total : 0; 954 955 for( j = 0; j < hmm->num_states; j++) 956 { 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; 960 } 961 } 962 963 /* estimate other matrices */ 964 for( k = 0; k < hmm->num_states; k++ ) 965 { 966 CvEHMM* ehmm = &(hmm->u.ehmm[k]); 967 968 for( i = 0; i < ehmm->num_states; i++) 969 { 970 float total = 0; 971 float inv_total; 972 for( j = 0; j < ehmm->num_states; j++) 973 { 974 total += ehmm->transP[i*ehmm->num_states + j]; 975 } 976 //assert( total ); 977 inv_total = total ? 1.f/total : 0; 978 979 for( j = 0; j < ehmm->num_states; j++) 980 { 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 ; 984 } 985 } 986 } 987 return CV_NO_ERR; 988 } 989 990 991 /*F/////////////////////////////////////////////////////////////////////////////////////// 992 // Name: MixSegmL2 993 // Purpose: The function implements the mixture segmentation of the states of the 994 // embedded HMM 995 // Context: used with the Viterbi training of the embedded HMM 996 // 997 // Parameters: 998 // obs_info_array 999 // num_img 1000 // hmm 1001 // Returns: void 1002 // 1003 // Notes: 1004 //F*/ 1005 static CvStatus CV_STDCALL 1006 icvMixSegmL2( CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm ) 1007 { 1008 int k, i, j, m; 1009 1010 CvEHMMState* state = hmm->u.ehmm[0].u.state; 1011 1012 1013 for (k = 0; k < num_img; k++) 1014 { 1015 int counter = 0; 1016 CvImgObsInfo* info = obs_info_array[k]; 1017 1018 for (i = 0; i < info->obs_y; i++) 1019 { 1020 for (j = 0; j < info->obs_x; j++, counter++) 1021 { 1022 int e_state = info->state[2 * counter + 1]; 1023 float min_dist; 1024 1025 min_dist = icvSquareDistance((info->obs) + (counter * info->obs_size), 1026 state[e_state].mu, info->obs_size); 1027 info->mix[counter] = 0; 1028 1029 for (m = 1; m < state[e_state].num_mix; m++) 1030 { 1031 float dist=icvSquareDistance( (info->obs) + (counter * info->obs_size), 1032 state[e_state].mu + m * info->obs_size, 1033 info->obs_size); 1034 if (dist < min_dist) 1035 { 1036 min_dist = dist; 1037 /* assign mixture with smallest distance */ 1038 info->mix[counter] = m; 1039 } 1040 } 1041 } 1042 } 1043 } 1044 return CV_NO_ERR; 1045 } 1046 1047 /* 1048 CvStatus icvMixSegmProb(CvImgObsInfo* obs_info, int num_img, CvEHMM* hmm ) 1049 { 1050 int k, i, j, m; 1051 1052 CvEHMMState* state = hmm->ehmm[0].state_info; 1053 1054 1055 for (k = 0; k < num_img; k++) 1056 { 1057 int counter = 0; 1058 CvImgObsInfo* info = obs_info + k; 1059 1060 for (i = 0; i < info->obs_y; i++) 1061 { 1062 for (j = 0; j < info->obs_x; j++, counter++) 1063 { 1064 int e_state = info->in_state[counter]; 1065 float max_prob; 1066 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], 1070 info->obs_size ); 1071 info->mix[counter] = 0; 1072 1073 for (m = 1; m < state[e_state].num_mix; m++) 1074 { 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], 1078 info->obs_size); 1079 if (prob > max_prob) 1080 { 1081 max_prob = prob; 1082 // assign mixture with greatest probability. 1083 info->mix[counter] = m; 1084 } 1085 } 1086 } 1087 } 1088 } 1089 1090 return CV_NO_ERR; 1091 } 1092 */ 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, 1097 float* prob ) 1098 { 1099 // memory allocation 1100 int i, j, last_obs; 1101 int m_HMMType = _CV_ERGODIC; /* _CV_CAUSAL or _CV_ERGODIC */ 1102 1103 int m_ProbType = prob_type; /* _CV_LAST_STATE or _CV_BEST_STATE */ 1104 1105 int m_minNumObs = min_num_obs; /*??*/ 1106 int m_maxNumObs = max_num_obs; /*??*/ 1107 1108 int m_numStates = num_states; 1109 1110 float* m_pi = (float*)cvAlloc( num_states* sizeof(float) ); 1111 CvMatr32f m_a = transP; 1112 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 1116 1117 //m_xl = start_obs; 1118 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++){ 1123 m_l[i] = muDur[i]; 1124 } 1125 } 1126 else{ 1127 m_d = NULL; 1128 m_l = NULL; 1129 } 1130 */ 1131 1132 CvMatr32f m_Gamma = icvCreateMatrix_32f( num_states, m_maxNumObs ); 1133 int* m_csi = (int*)cvAlloc( num_states * m_maxNumObs * sizeof(int) ); 1134 1135 //stores maximal result for every ending observation */ 1136 CvVect32f m_MaxGamma = prob; 1137 1138 1139 // assert( m_xl + max_num_obs <= num_obs ); 1140 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]; 1144 */ 1145 1146 /******************************************************************/ 1147 /* Viterbi initialization */ 1148 /* set initial state probabilities, in logarithmic scale */ 1149 for (i = 0; i < m_numStates; i++) 1150 { 1151 m_pi[i] = -BIG_FLT; 1152 } 1153 m_pi[0] = 0.0f; 1154 1155 for (i = 0; i < num_states; i++) 1156 { 1157 m_Gamma[0 * num_states + i] = m_pi[i] + m_b[0 * num_states + i]; 1158 m_csi[0 * num_states + i] = 0; 1159 } 1160 1161 /******************************************************************/ 1162 /* Viterbi recursion */ 1163 1164 if ( m_HMMType == _CV_CAUSAL ) //causal model 1165 { 1166 int t,j; 1167 1168 for (t = 1 ; t < m_maxNumObs; t++) 1169 { 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; 1173 1174 for (j = 1; j < num_states; j++) 1175 { 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]; 1178 1179 if ( prev > self ) 1180 { 1181 m_csi[t * num_states + j] = j-1; 1182 m_Gamma[t * num_states + j] = prev; 1183 } 1184 else 1185 { 1186 m_csi[t * num_states + j] = j; 1187 m_Gamma[t * num_states + j] = self; 1188 } 1189 1190 m_Gamma[t * num_states + j] = m_Gamma[t * num_states + j] + m_b[t * num_states + j]; 1191 } 1192 } 1193 } 1194 else if ( m_HMMType == _CV_ERGODIC ) //ergodic model 1195 { 1196 int t; 1197 for (t = 1 ; t < m_maxNumObs; t++) 1198 { 1199 for (j = 0; j < num_states; j++) 1200 { 1201 int i; 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; 1204 1205 for (i = 1; i < num_states; i++) 1206 { 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]) 1209 { 1210 m_Gamma[t * num_states + j] = currGamma; 1211 m_csi[t * num_states + j] = i; 1212 } 1213 } 1214 m_Gamma[t *num_states + j] = m_Gamma[t *num_states + j] + m_b[t * num_states + j]; 1215 } 1216 } 1217 } 1218 1219 for( last_obs = m_minNumObs-1, i = 0; last_obs < m_maxNumObs; last_obs++, i++ ) 1220 { 1221 int t; 1222 1223 /******************************************************************/ 1224 /* Viterbi termination */ 1225 1226 if ( m_ProbType == _CV_LAST_STATE ) 1227 { 1228 m_MaxGamma[i] = m_Gamma[last_obs * num_states + num_states - 1]; 1229 q[i][last_obs] = num_states - 1; 1230 } 1231 else if( m_ProbType == _CV_BEST_STATE ) 1232 { 1233 int k; 1234 q[i][last_obs] = 0; 1235 m_MaxGamma[i] = m_Gamma[last_obs * num_states + 0]; 1236 1237 for(k = 1; k < num_states; k++) 1238 { 1239 if ( m_Gamma[last_obs * num_states + k] > m_MaxGamma[i] ) 1240 { 1241 m_MaxGamma[i] = m_Gamma[last_obs * num_states + k]; 1242 q[i][last_obs] = k; 1243 } 1244 } 1245 } 1246 1247 /******************************************************************/ 1248 /* Viterbi backtracking */ 1249 for (t = last_obs-1; t >= 0; t--) 1250 { 1251 q[i][t] = m_csi[(t+1) * num_states + q[i][t+1] ]; 1252 } 1253 } 1254 1255 /* memory free */ 1256 cvFree( &m_pi ); 1257 cvFree( &m_csi ); 1258 icvDeleteMatrix( m_Gamma ); 1259 1260 return CV_NO_ERR; 1261 } 1262 1263 /*F/////////////////////////////////////////////////////////////////////////////////////// 1264 // Name: icvEViterbi 1265 // Purpose: The function calculates the embedded Viterbi algorithm 1266 // for 1 image 1267 // Context: 1268 // Parameters: 1269 // obs_info - observations 1270 // hmm - HMM 1271 // 1272 // Returns: the Embedded Viterbi probability (float) 1273 // and do state segmentation of observations 1274 // 1275 // Notes: 1276 //F*/ 1277 static float CV_STDCALL icvEViterbi( CvImgObsInfo* obs_info, CvEHMM* hmm ) 1278 { 1279 int i, j, counter; 1280 float log_likelihood; 1281 1282 float inv_obs_x = 1.f / obs_info->obs_x; 1283 1284 CvEHMMState* first_state = hmm->u.ehmm->u.state; 1285 1286 /* memory allocation for superB */ 1287 CvMatr32f superB = icvCreateMatrix_32f(hmm->num_states, obs_info->obs_y ); 1288 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) ); 1292 1293 for (i = 0; i < hmm->num_states; i++) 1294 { 1295 q[i] = (int**)cvAlloc( obs_info->obs_y * sizeof(int*) ); 1296 1297 for (j = 0; j < obs_info->obs_y ; j++) 1298 { 1299 q[i][j] = (int*)cvAlloc( obs_info->obs_x * sizeof(int) ); 1300 } 1301 } 1302 1303 /* start Viterbi segmentation */ 1304 for (i = 0; i < hmm->num_states; i++) 1305 { 1306 CvEHMM* ehmm = &(hmm->u.ehmm[i]); 1307 1308 for (j = 0; j < obs_info->obs_y; j++) 1309 { 1310 float max_gamma; 1311 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); 1317 1318 superB[j * hmm->num_states + i] = max_gamma * inv_obs_x; 1319 } 1320 } 1321 1322 /* perform global Viterbi segmentation (i.e. process higher-level HMM) */ 1323 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 ); 1328 1329 log_likelihood /= obs_info->obs_y ; 1330 1331 1332 counter = 0; 1333 /* assign new state to observation vectors */ 1334 for (i = 0; i < obs_info->obs_y; i++) 1335 { 1336 for (j = 0; j < obs_info->obs_x; j++, counter++) 1337 { 1338 int superstate = super_q[i]; 1339 int state = (int)(hmm->u.ehmm[superstate].u.state - first_state); 1340 1341 obs_info->state[2 * counter] = superstate; 1342 obs_info->state[2 * counter + 1] = state + q[superstate][i][j]; 1343 } 1344 } 1345 1346 /* memory deallocation for superB */ 1347 icvDeleteMatrix( superB ); 1348 1349 /*memory deallocation for q */ 1350 for (i = 0; i < hmm->num_states; i++) 1351 { 1352 for (j = 0; j < obs_info->obs_y ; j++) 1353 { 1354 cvFree( &q[i][j] ); 1355 } 1356 cvFree( &q[i] ); 1357 } 1358 1359 cvFree( &q ); 1360 cvFree( &super_q ); 1361 1362 return log_likelihood; 1363 } 1364 1365 static CvStatus CV_STDCALL 1366 icvEstimateHMMStateParams( CvImgObsInfo** obs_info_array, int num_img, CvEHMM* hmm ) 1367 { 1368 /* compute gamma, weights, means, vars */ 1369 int k, i, j, m; 1370 int total = 0; 1371 int vect_len = obs_info_array[0]->obs_size; 1372 1373 float start_log_var_val = LN2PI * vect_len; 1374 1375 CvVect32f tmp_vect = icvCreateVector_32f( vect_len ); 1376 1377 CvEHMMState* first_state = hmm->u.ehmm[0].u.state; 1378 1379 assert( sizeof(float) == sizeof(int) ); 1380 1381 for(i = 0; i < hmm->num_states; i++ ) 1382 { 1383 total+= hmm->u.ehmm[i].num_states; 1384 } 1385 1386 /***************Gamma***********************/ 1387 /* initialize gamma */ 1388 for( i = 0; i < total; i++ ) 1389 { 1390 for (m = 0; m < first_state[i].num_mix; m++) 1391 { 1392 ((int*)(first_state[i].weight))[m] = 0; 1393 } 1394 } 1395 1396 /* maybe gamma must be computed in mixsegm process ?? */ 1397 1398 /* compute gamma */ 1399 for (k = 0; k < num_img; k++) 1400 { 1401 CvImgObsInfo* info = obs_info_array[k]; 1402 int num_obs = info->obs_y * info->obs_x; 1403 1404 for (i = 0; i < num_obs; i++) 1405 { 1406 int state, mixture; 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; 1412 } 1413 } 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++) 1419 { 1420 memset( (void*)first_state[i].mu, 0, first_state[i].num_mix * vect_len * 1421 sizeof(float) ); 1422 memset( (void*)first_state[i].inv_var, 0, first_state[i].num_mix * vect_len * 1423 sizeof(float) ); 1424 } 1425 1426 /* compute sums */ 1427 for (i = 0; i < num_img; i++) 1428 { 1429 CvImgObsInfo* info = obs_info_array[i]; 1430 int total_obs = info->obs_x * info->obs_y; 1431 1432 float* vector = info->obs; 1433 1434 for (j = 0; j < total_obs; j++, vector+=vect_len ) 1435 { 1436 int state = info->state[2 * j + 1]; 1437 int mixture = info->mix[j]; 1438 1439 CvVect32f mean = first_state[state].mu + mixture * vect_len; 1440 CvVect32f mean2 = first_state[state].inv_var + mixture * vect_len; 1441 1442 icvAddVector_32f( mean, vector, mean, vect_len ); 1443 for( k = 0; k < vect_len; k++ ) 1444 mean2[k] += vector[k]*vector[k]; 1445 } 1446 } 1447 1448 /*compute the means and variances */ 1449 /* assume gamma already computed */ 1450 for (i = 0; i < total; i++) 1451 { 1452 CvEHMMState* state = &(first_state[i]); 1453 1454 for (m = 0; m < state->num_mix; m++) 1455 { 1456 int k; 1457 CvVect32f mu = state->mu + m * vect_len; 1458 CvVect32f invar = state->inv_var + m * vect_len; 1459 1460 if ( ((int*)state->weight)[m] > 1) 1461 { 1462 float inv_gamma = 1.f/((int*)(state->weight))[m]; 1463 1464 icvScaleVector_32f( mu, mu, vect_len, inv_gamma); 1465 icvScaleVector_32f( invar, invar, vect_len, inv_gamma); 1466 } 1467 1468 icvMulVectors_32f(mu, mu, tmp_vect, vect_len); 1469 icvSubVector_32f( invar, tmp_vect, invar, vect_len); 1470 1471 /* low bound of variance - 100 (Ara's experimental result) */ 1472 for( k = 0; k < vect_len; k++ ) 1473 { 1474 invar[k] = (invar[k] > 100.f) ? invar[k] : 100.f; 1475 } 1476 1477 /* compute log_var */ 1478 state->log_var_val[m] = start_log_var_val; 1479 for( k = 0; k < vect_len; k++ ) 1480 { 1481 state->log_var_val[m] += (float)log( invar[k] ); 1482 } 1483 1484 /* SMOLI 27.10.2000 */ 1485 state->log_var_val[m] *= 0.5; 1486 1487 1488 /* compute inv_var = 1/sqrt(2*variance) */ 1489 icvScaleVector_32f(invar, invar, vect_len, 2.f ); 1490 cvbInvSqrt( invar, invar, vect_len ); 1491 } 1492 } 1493 1494 /***************Weights***********************/ 1495 /* normilize gammas - i.e. compute mixture weights */ 1496 1497 //compute weights 1498 for (i = 0; i < total; i++) 1499 { 1500 int gamma_total = 0; 1501 float norm; 1502 1503 for (m = 0; m < first_state[i].num_mix; m++) 1504 { 1505 gamma_total += ((int*)(first_state[i].weight))[m]; 1506 } 1507 1508 norm = gamma_total ? (1.f/(float)gamma_total) : 0.f; 1509 1510 for (m = 0; m < first_state[i].num_mix; m++) 1511 { 1512 first_state[i].weight[m] = ((int*)(first_state[i].weight))[m] * norm; 1513 } 1514 } 1515 1516 icvDeleteVector( tmp_vect); 1517 return CV_NO_ERR; 1518 } 1519 1520 /* 1521 CvStatus icvLightingCorrection8uC1R( uchar* img, CvSize roi, int src_step ) 1522 { 1523 int i, j; 1524 int width = roi.width; 1525 int height = roi.height; 1526 1527 float x1, x2, y1, y2; 1528 int f[3] = {0, 0, 0}; 1529 float a[3] = {0, 0, 0}; 1530 1531 float h1; 1532 float h2; 1533 1534 float c1,c2; 1535 1536 float min = FLT_MAX; 1537 float max = -FLT_MAX; 1538 float correction; 1539 1540 float* float_img = icvAlloc( width * height * sizeof(float) ); 1541 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) 1546 1547 1548 // extract grayvalues 1549 for (i = 0; i < height; i++) 1550 { 1551 for (j = 0; j < width; j++) 1552 { 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]; 1556 } 1557 } 1558 1559 h1 = (float)f[0] * (float)x1 / (float)width; 1560 h2 = (float)f[0] * (float)y1 / (float)height; 1561 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; 1566 1567 for (i = 0; i < height; i++) 1568 { 1569 for (j = 0; j < width; j++) 1570 { 1571 1572 correction = a[0] + a[1]*(float)i + a[2]*(float)j; 1573 1574 float_img[i*width + j] = img[i*src_step + j] - correction; 1575 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]; 1578 } 1579 } 1580 1581 //rescaling to the range 0:255 1582 c2 = 0; 1583 if (max == min) 1584 c2 = 255.0f; 1585 else 1586 c2 = 255.0f/(float)(max - min); 1587 1588 c1 = (-(float)min)*c2; 1589 1590 for (i = 0; i < height; i++) 1591 { 1592 for (j = 0; j < width; j++) 1593 { 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; 1598 } 1599 } 1600 1601 cvFree( &float_img ); 1602 return CV_NO_ERR; 1603 } 1604 1605 1606 CvStatus icvLightingCorrection( icvImage* img ) 1607 { 1608 CvSize roi; 1609 if ( img->type != IPL_DEPTH_8U || img->channels != 1 ) 1610 return CV_BADFACTOR_ERR; 1611 1612 roi = _cvSize( img->roi.width, img->roi.height ); 1613 1614 return _cvLightingCorrection8uC1R( img->data + img->roi.y * img->step + img->roi.x, 1615 roi, img->step ); 1616 1617 } 1618 1619 */ 1620 1621 CV_IMPL CvEHMM* 1622 cvCreate2DHMM( int *state_number, int *num_mix, int obs_size ) 1623 { 1624 CvEHMM* hmm = 0; 1625 1626 CV_FUNCNAME( "cvCreate2DHMM" ); 1627 1628 __BEGIN__; 1629 1630 IPPI_CALL( icvCreate2DHMM( &hmm, state_number, num_mix, obs_size )); 1631 1632 __END__; 1633 1634 return hmm; 1635 } 1636 1637 CV_IMPL void 1638 cvRelease2DHMM( CvEHMM ** hmm ) 1639 { 1640 CV_FUNCNAME( "cvRelease2DHMM" ); 1641 1642 __BEGIN__; 1643 1644 IPPI_CALL( icvRelease2DHMM( hmm )); 1645 __END__; 1646 } 1647 1648 CV_IMPL CvImgObsInfo* 1649 cvCreateObsInfo( CvSize num_obs, int obs_size ) 1650 { 1651 CvImgObsInfo *obs_info = 0; 1652 1653 CV_FUNCNAME( "cvCreateObsInfo" ); 1654 1655 __BEGIN__; 1656 1657 IPPI_CALL( icvCreateObsInfo( &obs_info, num_obs, obs_size )); 1658 1659 __END__; 1660 1661 return obs_info; 1662 } 1663 1664 CV_IMPL void 1665 cvReleaseObsInfo( CvImgObsInfo ** obs_info ) 1666 { 1667 CV_FUNCNAME( "cvReleaseObsInfo" ); 1668 1669 __BEGIN__; 1670 1671 IPPI_CALL( icvReleaseObsInfo( obs_info )); 1672 1673 __END__; 1674 } 1675 1676 1677 CV_IMPL void 1678 cvUniformImgSegm( CvImgObsInfo * obs_info, CvEHMM * hmm ) 1679 { 1680 CV_FUNCNAME( "cvUniformImgSegm" ); 1681 1682 __BEGIN__; 1683 1684 IPPI_CALL( icvUniformImgSegm( obs_info, hmm )); 1685 __CLEANUP__; 1686 __END__; 1687 } 1688 1689 CV_IMPL void 1690 cvInitMixSegm( CvImgObsInfo ** obs_info_array, int num_img, CvEHMM * hmm ) 1691 { 1692 CV_FUNCNAME( "cvInitMixSegm" ); 1693 1694 __BEGIN__; 1695 1696 IPPI_CALL( icvInitMixSegm( obs_info_array, num_img, hmm )); 1697 1698 __END__; 1699 } 1700 1701 CV_IMPL void 1702 cvEstimateHMMStateParams( CvImgObsInfo ** obs_info_array, int num_img, CvEHMM * hmm ) 1703 { 1704 CV_FUNCNAME( "cvEstimateHMMStateParams" ); 1705 1706 __BEGIN__; 1707 1708 IPPI_CALL( icvEstimateHMMStateParams( obs_info_array, num_img, hmm )); 1709 1710 __END__; 1711 } 1712 1713 CV_IMPL void 1714 cvEstimateTransProb( CvImgObsInfo ** obs_info_array, int num_img, CvEHMM * hmm ) 1715 { 1716 CV_FUNCNAME( "cvEstimateTransProb" ); 1717 1718 __BEGIN__; 1719 1720 IPPI_CALL( icvEstimateTransProb( obs_info_array, num_img, hmm )); 1721 1722 __END__; 1723 1724 } 1725 1726 CV_IMPL void 1727 cvEstimateObsProb( CvImgObsInfo * obs_info, CvEHMM * hmm ) 1728 { 1729 CV_FUNCNAME( "cvEstimateObsProb" ); 1730 1731 __BEGIN__; 1732 1733 IPPI_CALL( icvEstimateObsProb( obs_info, hmm )); 1734 1735 __END__; 1736 } 1737 1738 CV_IMPL float 1739 cvEViterbi( CvImgObsInfo * obs_info, CvEHMM * hmm ) 1740 { 1741 float result = FLT_MAX; 1742 1743 CV_FUNCNAME( "cvEViterbi" ); 1744 1745 __BEGIN__; 1746 1747 if( (obs_info == NULL) || (hmm == NULL) ) 1748 CV_ERROR( CV_BadDataPtr, "Null pointer." ); 1749 1750 result = icvEViterbi( obs_info, hmm ); 1751 1752 __END__; 1753 1754 return result; 1755 } 1756 1757 CV_IMPL void 1758 cvMixSegmL2( CvImgObsInfo ** obs_info_array, int num_img, CvEHMM * hmm ) 1759 { 1760 CV_FUNCNAME( "cvMixSegmL2" ); 1761 1762 __BEGIN__; 1763 1764 IPPI_CALL( icvMixSegmL2( obs_info_array, num_img, hmm )); 1765 1766 __END__; 1767 } 1768 1769 /* End of file */ 1770 1771