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#include "_cvaux.h" 42 43//*F/////////////////////////////////////////////////////////////////////////////////////// 44// Name: icvImgToObs_DCT_8u32f_C1R 45// Purpose: The function takes as input an image and returns the sequnce of observations 46// to be used with an embedded HMM; Each observation is top-left block of DCT 47// coefficient matrix. 48// Context: 49// Parameters: img - pointer to the original image ROI 50// imgStep - full row width of the image in bytes 51// roi - width and height of ROI in pixels 52// obs - pointer to resultant observation vectors 53// dctSize - size of the block for which DCT is calculated 54// obsSize - size of top-left block of DCT coeffs matrix, which is treated 55// as observation. Each observation vector consists of 56// obsSize.width * obsSize.height floats. 57// The following conditions should be satisfied: 58// 0 < objSize.width <= dctSize.width, 59// 0 < objSize.height <= dctSize.height. 60// delta - dctBlocks are overlapped and this parameter specifies horizontal 61// and vertical shift. 62// Returns: 63// CV_NO_ERR or error code 64// Notes: 65// The algorithm is following: 66// 1. First, number of observation vectors per row and per column are calculated: 67// 68// Nx = floor((roi.width - dctSize.width + delta.width)/delta.width); 69// Ny = floor((roi.height - dctSize.height + delta.height)/delta.height); 70// 71// So, total number of observation vectors is Nx*Ny, and total size of 72// array obs must be >= Nx*Ny*obsSize.width*obsSize.height*sizeof(float). 73// 2. Observation vectors are calculated in the following loop 74// ( actual implementation may be different ), where 75// I[x1:x2,y1:y2] means block of pixels from source image with 76// x1 <= x < x2, y1 <= y < y2, 77// D[x1:x2,y1:y2] means sub matrix of DCT matrix D. 78// O[x,y] means observation vector that corresponds to position 79// (x*delta.width,y*delta.height) in the source image 80// ( all indices are counted from 0 ). 81// 82// for( y = 0; y < Ny; y++ ) 83// { 84// for( x = 0; x < Nx; x++ ) 85// { 86// D = DCT(I[x*delta.width : x*delta.width + dctSize.width, 87// y*delta.height : y*delta.height + dctSize.height]); 88// O[x,y] = D[0:obsSize.width, 0:obsSize.height]; 89// } 90// } 91//F*/ 92 93/*comment out the following line to make DCT be calculated in floating-point arithmetics*/ 94//#define _CV_INT_DCT 95 96/* for integer DCT only */ 97#define DCT_SCALE 15 98 99#ifdef _CV_INT_DCT 100typedef int work_t; 101 102#define DESCALE CV_DESCALE 103#define SCALE(x) CV_FLT_TO_FIX((x),DCT_SCALE) 104#else 105typedef float work_t; 106 107#define DESCALE(x,n) (float)(x) 108#define SCALE(x) (float)(x) 109#endif 110 111/* calculate dct transform matrix */ 112static void icvCalcDCTMatrix( work_t * cfs, int n ); 113 114#define MAX_DCT_SIZE 32 115 116static CvStatus CV_STDCALL 117icvImgToObs_DCT_8u32f_C1R( uchar * img, int imgStep, CvSize roi, 118 float *obs, CvSize dctSize, 119 CvSize obsSize, CvSize delta ) 120{ 121 /* dct transform matrices: horizontal and vertical */ 122 work_t tab_x[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2]; 123 work_t tab_y[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2]; 124 125 /* temporary buffers for dct */ 126 work_t temp0[MAX_DCT_SIZE * 4]; 127 work_t temp1[MAX_DCT_SIZE * 4]; 128 work_t *buffer = 0; 129 work_t *buf_limit; 130 131 double s; 132 133 int y; 134 int Nx, Ny; 135 136 int n1 = dctSize.height, m1 = n1 / 2; 137 int n2 = dctSize.width, m2 = n2 / 2; 138 139 if( !img || !obs ) 140 return CV_NULLPTR_ERR; 141 142 if( roi.width <= 0 || roi.height <= 0 ) 143 return CV_BADSIZE_ERR; 144 145 if( delta.width <= 0 || delta.height <= 0 ) 146 return CV_BADRANGE_ERR; 147 148 if( obsSize.width <= 0 || dctSize.width < obsSize.width || 149 obsSize.height <= 0 || dctSize.height < obsSize.height ) 150 return CV_BADRANGE_ERR; 151 152 if( dctSize.width > MAX_DCT_SIZE || dctSize.height > MAX_DCT_SIZE ) 153 return CV_BADRANGE_ERR; 154 155 Nx = (roi.width - dctSize.width + delta.width) / delta.width; 156 Ny = (roi.height - dctSize.height + delta.height) / delta.height; 157 158 if( Nx <= 0 || Ny <= 0 ) 159 return CV_BADRANGE_ERR; 160 161 buffer = (work_t *)cvAlloc( roi.width * obsSize.height * sizeof( buffer[0] )); 162 if( !buffer ) 163 return CV_OUTOFMEM_ERR; 164 165 icvCalcDCTMatrix( tab_x, dctSize.width ); 166 icvCalcDCTMatrix( tab_y, dctSize.height ); 167 168 buf_limit = buffer + obsSize.height * roi.width; 169 170 for( y = 0; y < Ny; y++, img += delta.height * imgStep ) 171 { 172 int x, i, j, k; 173 work_t k0 = 0; 174 175 /* do transfroms for each column. Calc only first obsSize.height DCT coefficients */ 176 for( x = 0; x < roi.width; x++ ) 177 { 178 float is = 0; 179 work_t *buf = buffer + x; 180 work_t *tab = tab_y + 2; 181 182 if( n1 & 1 ) 183 { 184 is = img[x + m1 * imgStep]; 185 k0 = ((work_t) is) * tab[-1]; 186 } 187 188 /* first coefficient */ 189 for( j = 0; j < m1; j++ ) 190 { 191 float t0 = img[x + j * imgStep]; 192 float t1 = img[x + (n1 - 1 - j) * imgStep]; 193 float t2 = t0 + t1; 194 195 t0 -= t1; 196 temp0[j] = (work_t) t2; 197 is += t2; 198 temp1[j] = (work_t) t0; 199 } 200 201 buf[0] = DESCALE( is * tab[-2], PASS1_SHIFT ); 202 if( (buf += roi.width) >= buf_limit ) 203 continue; 204 205 /* other coefficients */ 206 for( ;; ) 207 { 208 s = 0; 209 210 for( k = 0; k < m1; k++ ) 211 s += temp1[k] * tab[k]; 212 213 buf[0] = DESCALE( s, PASS1_SHIFT ); 214 if( (buf += roi.width) >= buf_limit ) 215 break; 216 217 tab += m1; 218 s = 0; 219 220 if( n1 & 1 ) 221 { 222 k0 = -k0; 223 s = k0; 224 } 225 for( k = 0; k < m1; k++ ) 226 s += temp0[k] * tab[k]; 227 228 buf[0] = DESCALE( s, PASS1_SHIFT ); 229 tab += m1; 230 231 if( (buf += roi.width) >= buf_limit ) 232 break; 233 } 234 } 235 236 k0 = 0; 237 238 /* do transforms for rows. */ 239 for( x = 0; x + dctSize.width <= roi.width; x += delta.width ) 240 { 241 for( i = 0; i < obsSize.height; i++ ) 242 { 243 work_t *buf = buffer + x + roi.width * i; 244 work_t *tab = tab_x + 2; 245 float *obs_limit = obs + obsSize.width; 246 247 s = 0; 248 249 if( n2 & 1 ) 250 { 251 s = buf[m2]; 252 k0 = (work_t) (s * tab[-1]); 253 } 254 255 /* first coefficient */ 256 for( j = 0; j < m2; j++ ) 257 { 258 work_t t0 = buf[j]; 259 work_t t1 = buf[n2 - 1 - j]; 260 work_t t2 = t0 + t1; 261 262 t0 -= t1; 263 temp0[j] = (work_t) t2; 264 s += t2; 265 temp1[j] = (work_t) t0; 266 } 267 268 *obs++ = (float) DESCALE( s * tab[-2], PASS2_SHIFT ); 269 270 if( obs == obs_limit ) 271 continue; 272 273 /* other coefficients */ 274 for( ;; ) 275 { 276 s = 0; 277 278 for( k = 0; k < m2; k++ ) 279 s += temp1[k] * tab[k]; 280 281 obs[0] = (float) DESCALE( s, PASS2_SHIFT ); 282 if( ++obs == obs_limit ) 283 break; 284 285 tab += m2; 286 287 s = 0; 288 289 if( n2 & 1 ) 290 { 291 k0 = -k0; 292 s = k0; 293 } 294 for( k = 0; k < m2; k++ ) 295 s += temp0[k] * tab[k]; 296 obs[0] = (float) DESCALE( s, PASS2_SHIFT ); 297 298 tab += m2; 299 if( ++obs == obs_limit ) 300 break; 301 } 302 } 303 } 304 } 305 306 cvFree( &buffer ); 307 return CV_NO_ERR; 308} 309 310 311static CvStatus CV_STDCALL 312icvImgToObs_DCT_32f_C1R( float * img, int imgStep, CvSize roi, 313 float *obs, CvSize dctSize, 314 CvSize obsSize, CvSize delta ) 315{ 316 /* dct transform matrices: horizontal and vertical */ 317 work_t tab_x[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2]; 318 work_t tab_y[MAX_DCT_SIZE * MAX_DCT_SIZE / 2 + 2]; 319 320 /* temporary buffers for dct */ 321 work_t temp0[MAX_DCT_SIZE * 4]; 322 work_t temp1[MAX_DCT_SIZE * 4]; 323 work_t *buffer = 0; 324 work_t *buf_limit; 325 326 double s; 327 328 int y; 329 int Nx, Ny; 330 331 int n1 = dctSize.height, m1 = n1 / 2; 332 int n2 = dctSize.width, m2 = n2 / 2; 333 334 if( !img || !obs ) 335 return CV_NULLPTR_ERR; 336 337 if( roi.width <= 0 || roi.height <= 0 ) 338 return CV_BADSIZE_ERR; 339 340 if( delta.width <= 0 || delta.height <= 0 ) 341 return CV_BADRANGE_ERR; 342 343 if( obsSize.width <= 0 || dctSize.width < obsSize.width || 344 obsSize.height <= 0 || dctSize.height < obsSize.height ) 345 return CV_BADRANGE_ERR; 346 347 if( dctSize.width > MAX_DCT_SIZE || dctSize.height > MAX_DCT_SIZE ) 348 return CV_BADRANGE_ERR; 349 350 Nx = (roi.width - dctSize.width + delta.width) / delta.width; 351 Ny = (roi.height - dctSize.height + delta.height) / delta.height; 352 353 if( Nx <= 0 || Ny <= 0 ) 354 return CV_BADRANGE_ERR; 355 356 buffer = (work_t *)cvAlloc( roi.width * obsSize.height * sizeof( buffer[0] )); 357 if( !buffer ) 358 return CV_OUTOFMEM_ERR; 359 360 icvCalcDCTMatrix( tab_x, dctSize.width ); 361 icvCalcDCTMatrix( tab_y, dctSize.height ); 362 363 buf_limit = buffer + obsSize.height * roi.width; 364 365 imgStep /= sizeof(img[0]); 366 367 for( y = 0; y < Ny; y++, img += delta.height * imgStep ) 368 { 369 int x, i, j, k; 370 work_t k0 = 0; 371 372 /* do transfroms for each column. Calc only first obsSize.height DCT coefficients */ 373 for( x = 0; x < roi.width; x++ ) 374 { 375 float is = 0; 376 work_t *buf = buffer + x; 377 work_t *tab = tab_y + 2; 378 379 if( n1 & 1 ) 380 { 381 is = img[x + m1 * imgStep]; 382 k0 = ((work_t) is) * tab[-1]; 383 } 384 385 /* first coefficient */ 386 for( j = 0; j < m1; j++ ) 387 { 388 float t0 = img[x + j * imgStep]; 389 float t1 = img[x + (n1 - 1 - j) * imgStep]; 390 float t2 = t0 + t1; 391 392 t0 -= t1; 393 temp0[j] = (work_t) t2; 394 is += t2; 395 temp1[j] = (work_t) t0; 396 } 397 398 buf[0] = DESCALE( is * tab[-2], PASS1_SHIFT ); 399 if( (buf += roi.width) >= buf_limit ) 400 continue; 401 402 /* other coefficients */ 403 for( ;; ) 404 { 405 s = 0; 406 407 for( k = 0; k < m1; k++ ) 408 s += temp1[k] * tab[k]; 409 410 buf[0] = DESCALE( s, PASS1_SHIFT ); 411 if( (buf += roi.width) >= buf_limit ) 412 break; 413 414 tab += m1; 415 s = 0; 416 417 if( n1 & 1 ) 418 { 419 k0 = -k0; 420 s = k0; 421 } 422 for( k = 0; k < m1; k++ ) 423 s += temp0[k] * tab[k]; 424 425 buf[0] = DESCALE( s, PASS1_SHIFT ); 426 tab += m1; 427 428 if( (buf += roi.width) >= buf_limit ) 429 break; 430 } 431 } 432 433 k0 = 0; 434 435 /* do transforms for rows. */ 436 for( x = 0; x + dctSize.width <= roi.width; x += delta.width ) 437 { 438 for( i = 0; i < obsSize.height; i++ ) 439 { 440 work_t *buf = buffer + x + roi.width * i; 441 work_t *tab = tab_x + 2; 442 float *obs_limit = obs + obsSize.width; 443 444 s = 0; 445 446 if( n2 & 1 ) 447 { 448 s = buf[m2]; 449 k0 = (work_t) (s * tab[-1]); 450 } 451 452 /* first coefficient */ 453 for( j = 0; j < m2; j++ ) 454 { 455 work_t t0 = buf[j]; 456 work_t t1 = buf[n2 - 1 - j]; 457 work_t t2 = t0 + t1; 458 459 t0 -= t1; 460 temp0[j] = (work_t) t2; 461 s += t2; 462 temp1[j] = (work_t) t0; 463 } 464 465 *obs++ = (float) DESCALE( s * tab[-2], PASS2_SHIFT ); 466 467 if( obs == obs_limit ) 468 continue; 469 470 /* other coefficients */ 471 for( ;; ) 472 { 473 s = 0; 474 475 for( k = 0; k < m2; k++ ) 476 s += temp1[k] * tab[k]; 477 478 obs[0] = (float) DESCALE( s, PASS2_SHIFT ); 479 if( ++obs == obs_limit ) 480 break; 481 482 tab += m2; 483 484 s = 0; 485 486 if( n2 & 1 ) 487 { 488 k0 = -k0; 489 s = k0; 490 } 491 for( k = 0; k < m2; k++ ) 492 s += temp0[k] * tab[k]; 493 obs[0] = (float) DESCALE( s, PASS2_SHIFT ); 494 495 tab += m2; 496 if( ++obs == obs_limit ) 497 break; 498 } 499 } 500 } 501 } 502 503 cvFree( &buffer ); 504 return CV_NO_ERR; 505} 506 507 508static void 509icvCalcDCTMatrix( work_t * cfs, int n ) 510{ 511 static const double sqrt2 = 1.4142135623730950488016887242097; 512 static const double pi = 3.1415926535897932384626433832795; 513 514 static const double sincos[16 * 2] = { 515 1.00000000000000000, 0.00000000000000006, 516 0.70710678118654746, 0.70710678118654757, 517 0.49999999999999994, 0.86602540378443871, 518 0.38268343236508978, 0.92387953251128674, 519 0.30901699437494740, 0.95105651629515353, 520 0.25881904510252074, 0.96592582628906831, 521 0.22252093395631439, 0.97492791218182362, 522 0.19509032201612825, 0.98078528040323043, 523 0.17364817766693033, 0.98480775301220802, 524 0.15643446504023087, 0.98768834059513777, 525 0.14231483827328514, 0.98982144188093268, 526 0.13052619222005157, 0.99144486137381038, 527 0.12053668025532305, 0.99270887409805397, 528 0.11196447610330786, 0.99371220989324260, 529 0.10452846326765346, 0.99452189536827329, 530 0.09801714032956060, 0.99518472667219693, 531 }; 532 533#define ROTATE( c, s, dc, ds ) \ 534 { \ 535 t = c*dc - s*ds; \ 536 s = c*ds + s*dc; \ 537 c = t; \ 538 } 539 540#define WRITE2( j, a, b ) \ 541 { \ 542 cfs[j] = SCALE(a); \ 543 cfs2[j] = SCALE(b); \ 544 } 545 546 double t, scale = 1. / sqrt( (double)n ); 547 int i, j, m = n / 2; 548 549 cfs[0] = SCALE( scale ); 550 scale *= sqrt2; 551 cfs[1] = SCALE( scale ); 552 cfs += 2 - m; 553 554 if( n > 1 ) 555 { 556 double a0, b0; 557 double da0, db0; 558 work_t *cfs2 = cfs + m * n; 559 560 if( n <= 16 ) 561 { 562 da0 = a0 = sincos[2 * n - 1]; 563 db0 = b0 = sincos[2 * n - 2]; 564 } 565 else 566 { 567 t = pi / (2 * n); 568 da0 = a0 = cos( t ); 569 db0 = b0 = sin( t ); 570 } 571 572 /* other rows */ 573 for( i = 1; i <= m; i++ ) 574 { 575 double a = a0 * scale; 576 double b = b0 * scale; 577 double da = a0 * a0 - b0 * b0; 578 double db = a0 * b0 + a0 * b0; 579 580 cfs += m; 581 cfs2 -= m; 582 583 for( j = 0; j < m; j += 2 ) 584 { 585 WRITE2( j, a, b ); 586 ROTATE( a, b, da, db ); 587 if( j + 1 < m ) 588 { 589 WRITE2( j + 1, a, -b ); 590 ROTATE( a, b, da, db ); 591 } 592 } 593 594 ROTATE( a0, b0, da0, db0 ); 595 } 596 } 597#undef ROTATE 598#undef WRITE2 599} 600 601 602CV_IMPL void 603cvImgToObs_DCT( const void* arr, float *obs, CvSize dctSize, 604 CvSize obsSize, CvSize delta ) 605{ 606 CV_FUNCNAME( "cvImgToObs_DCT" ); 607 608 __BEGIN__; 609 610 CvMat stub, *mat = (CvMat*)arr; 611 612 CV_CALL( mat = cvGetMat( arr, &stub )); 613 614 switch( CV_MAT_TYPE( mat->type )) 615 { 616 case CV_8UC1: 617 IPPI_CALL( icvImgToObs_DCT_8u32f_C1R( mat->data.ptr, mat->step, 618 cvGetMatSize(mat), obs, 619 dctSize, obsSize, delta )); 620 break; 621 case CV_32FC1: 622 IPPI_CALL( icvImgToObs_DCT_32f_C1R( mat->data.fl, mat->step, 623 cvGetMatSize(mat), obs, 624 dctSize, obsSize, delta )); 625 break; 626 default: 627 CV_ERROR( CV_StsUnsupportedFormat, "" ); 628 } 629 630 __END__; 631} 632 633 634/* End of file. */ 635