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 "_cxcore.h" 43 44// On Win64 (IA64) optimized versions of DFT and DCT fail the tests 45#if defined WIN64 && !defined EM64T 46#pragma optimize("", off) 47#endif 48 49icvDFTInitAlloc_C_32fc_t icvDFTInitAlloc_C_32fc_p = 0; 50icvDFTFree_C_32fc_t icvDFTFree_C_32fc_p = 0; 51icvDFTGetBufSize_C_32fc_t icvDFTGetBufSize_C_32fc_p = 0; 52icvDFTFwd_CToC_32fc_t icvDFTFwd_CToC_32fc_p = 0; 53icvDFTInv_CToC_32fc_t icvDFTInv_CToC_32fc_p = 0; 54 55icvDFTInitAlloc_C_64fc_t icvDFTInitAlloc_C_64fc_p = 0; 56icvDFTFree_C_64fc_t icvDFTFree_C_64fc_p = 0; 57icvDFTGetBufSize_C_64fc_t icvDFTGetBufSize_C_64fc_p = 0; 58icvDFTFwd_CToC_64fc_t icvDFTFwd_CToC_64fc_p = 0; 59icvDFTInv_CToC_64fc_t icvDFTInv_CToC_64fc_p = 0; 60 61icvDFTInitAlloc_R_32f_t icvDFTInitAlloc_R_32f_p = 0; 62icvDFTFree_R_32f_t icvDFTFree_R_32f_p = 0; 63icvDFTGetBufSize_R_32f_t icvDFTGetBufSize_R_32f_p = 0; 64icvDFTFwd_RToPack_32f_t icvDFTFwd_RToPack_32f_p = 0; 65icvDFTInv_PackToR_32f_t icvDFTInv_PackToR_32f_p = 0; 66 67icvDFTInitAlloc_R_64f_t icvDFTInitAlloc_R_64f_p = 0; 68icvDFTFree_R_64f_t icvDFTFree_R_64f_p = 0; 69icvDFTGetBufSize_R_64f_t icvDFTGetBufSize_R_64f_p = 0; 70icvDFTFwd_RToPack_64f_t icvDFTFwd_RToPack_64f_p = 0; 71icvDFTInv_PackToR_64f_t icvDFTInv_PackToR_64f_p = 0; 72 73/*icvDCTFwdInitAlloc_32f_t icvDCTFwdInitAlloc_32f_p = 0; 74icvDCTFwdFree_32f_t icvDCTFwdFree_32f_p = 0; 75icvDCTFwdGetBufSize_32f_t icvDCTFwdGetBufSize_32f_p = 0; 76icvDCTFwd_32f_t icvDCTFwd_32f_p = 0; 77 78icvDCTInvInitAlloc_32f_t icvDCTInvInitAlloc_32f_p = 0; 79icvDCTInvFree_32f_t icvDCTInvFree_32f_p = 0; 80icvDCTInvGetBufSize_32f_t icvDCTInvGetBufSize_32f_p = 0; 81icvDCTInv_32f_t icvDCTInv_32f_p = 0; 82 83icvDCTFwdInitAlloc_64f_t icvDCTFwdInitAlloc_64f_p = 0; 84icvDCTFwdFree_64f_t icvDCTFwdFree_64f_p = 0; 85icvDCTFwdGetBufSize_64f_t icvDCTFwdGetBufSize_64f_p = 0; 86icvDCTFwd_64f_t icvDCTFwd_64f_p = 0; 87 88icvDCTInvInitAlloc_64f_t icvDCTInvInitAlloc_64f_p = 0; 89icvDCTInvFree_64f_t icvDCTInvFree_64f_p = 0; 90icvDCTInvGetBufSize_64f_t icvDCTInvGetBufSize_64f_p = 0; 91icvDCTInv_64f_t icvDCTInv_64f_p = 0;*/ 92 93/****************************************************************************************\ 94 Discrete Fourier Transform 95\****************************************************************************************/ 96 97#define CV_MAX_LOCAL_DFT_SIZE (1 << 15) 98 99static const uchar log2tab[] = { 0, 0, 1, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0 }; 100static int icvlog2( int n ) 101{ 102 int m = 0; 103 int f = (n >= (1 << 16))*16; 104 n >>= f; 105 m += f; 106 f = (n >= (1 << 8))*8; 107 n >>= f; 108 m += f; 109 f = (n >= (1 << 4))*4; 110 n >>= f; 111 return m + f + log2tab[n]; 112} 113 114static unsigned char icvRevTable[] = 115{ 116 0x00,0x80,0x40,0xc0,0x20,0xa0,0x60,0xe0,0x10,0x90,0x50,0xd0,0x30,0xb0,0x70,0xf0, 117 0x08,0x88,0x48,0xc8,0x28,0xa8,0x68,0xe8,0x18,0x98,0x58,0xd8,0x38,0xb8,0x78,0xf8, 118 0x04,0x84,0x44,0xc4,0x24,0xa4,0x64,0xe4,0x14,0x94,0x54,0xd4,0x34,0xb4,0x74,0xf4, 119 0x0c,0x8c,0x4c,0xcc,0x2c,0xac,0x6c,0xec,0x1c,0x9c,0x5c,0xdc,0x3c,0xbc,0x7c,0xfc, 120 0x02,0x82,0x42,0xc2,0x22,0xa2,0x62,0xe2,0x12,0x92,0x52,0xd2,0x32,0xb2,0x72,0xf2, 121 0x0a,0x8a,0x4a,0xca,0x2a,0xaa,0x6a,0xea,0x1a,0x9a,0x5a,0xda,0x3a,0xba,0x7a,0xfa, 122 0x06,0x86,0x46,0xc6,0x26,0xa6,0x66,0xe6,0x16,0x96,0x56,0xd6,0x36,0xb6,0x76,0xf6, 123 0x0e,0x8e,0x4e,0xce,0x2e,0xae,0x6e,0xee,0x1e,0x9e,0x5e,0xde,0x3e,0xbe,0x7e,0xfe, 124 0x01,0x81,0x41,0xc1,0x21,0xa1,0x61,0xe1,0x11,0x91,0x51,0xd1,0x31,0xb1,0x71,0xf1, 125 0x09,0x89,0x49,0xc9,0x29,0xa9,0x69,0xe9,0x19,0x99,0x59,0xd9,0x39,0xb9,0x79,0xf9, 126 0x05,0x85,0x45,0xc5,0x25,0xa5,0x65,0xe5,0x15,0x95,0x55,0xd5,0x35,0xb5,0x75,0xf5, 127 0x0d,0x8d,0x4d,0xcd,0x2d,0xad,0x6d,0xed,0x1d,0x9d,0x5d,0xdd,0x3d,0xbd,0x7d,0xfd, 128 0x03,0x83,0x43,0xc3,0x23,0xa3,0x63,0xe3,0x13,0x93,0x53,0xd3,0x33,0xb3,0x73,0xf3, 129 0x0b,0x8b,0x4b,0xcb,0x2b,0xab,0x6b,0xeb,0x1b,0x9b,0x5b,0xdb,0x3b,0xbb,0x7b,0xfb, 130 0x07,0x87,0x47,0xc7,0x27,0xa7,0x67,0xe7,0x17,0x97,0x57,0xd7,0x37,0xb7,0x77,0xf7, 131 0x0f,0x8f,0x4f,0xcf,0x2f,0xaf,0x6f,0xef,0x1f,0x9f,0x5f,0xdf,0x3f,0xbf,0x7f,0xff 132}; 133 134static const double icvDxtTab[][2] = 135{ 136{ 1.00000000000000000, 0.00000000000000000 }, 137{-1.00000000000000000, 0.00000000000000000 }, 138{ 0.00000000000000000, 1.00000000000000000 }, 139{ 0.70710678118654757, 0.70710678118654746 }, 140{ 0.92387953251128674, 0.38268343236508978 }, 141{ 0.98078528040323043, 0.19509032201612825 }, 142{ 0.99518472667219693, 0.09801714032956060 }, 143{ 0.99879545620517241, 0.04906767432741802 }, 144{ 0.99969881869620425, 0.02454122852291229 }, 145{ 0.99992470183914450, 0.01227153828571993 }, 146{ 0.99998117528260111, 0.00613588464915448 }, 147{ 0.99999529380957619, 0.00306795676296598 }, 148{ 0.99999882345170188, 0.00153398018628477 }, 149{ 0.99999970586288223, 0.00076699031874270 }, 150{ 0.99999992646571789, 0.00038349518757140 }, 151{ 0.99999998161642933, 0.00019174759731070 }, 152{ 0.99999999540410733, 0.00009587379909598 }, 153{ 0.99999999885102686, 0.00004793689960307 }, 154{ 0.99999999971275666, 0.00002396844980842 }, 155{ 0.99999999992818922, 0.00001198422490507 }, 156{ 0.99999999998204725, 0.00000599211245264 }, 157{ 0.99999999999551181, 0.00000299605622633 }, 158{ 0.99999999999887801, 0.00000149802811317 }, 159{ 0.99999999999971945, 0.00000074901405658 }, 160{ 0.99999999999992983, 0.00000037450702829 }, 161{ 0.99999999999998246, 0.00000018725351415 }, 162{ 0.99999999999999567, 0.00000009362675707 }, 163{ 0.99999999999999889, 0.00000004681337854 }, 164{ 0.99999999999999978, 0.00000002340668927 }, 165{ 0.99999999999999989, 0.00000001170334463 }, 166{ 1.00000000000000000, 0.00000000585167232 }, 167{ 1.00000000000000000, 0.00000000292583616 } 168}; 169 170#define icvBitRev(i,shift) \ 171 ((int)((((unsigned)icvRevTable[(i)&255] << 24)+ \ 172 ((unsigned)icvRevTable[((i)>> 8)&255] << 16)+ \ 173 ((unsigned)icvRevTable[((i)>>16)&255] << 8)+ \ 174 ((unsigned)icvRevTable[((i)>>24)])) >> (shift))) 175 176static int 177icvDFTFactorize( int n, int* factors ) 178{ 179 int nf = 0, f, i, j; 180 181 if( n <= 5 ) 182 { 183 factors[0] = n; 184 return 1; 185 } 186 187 f = (((n - 1)^n)+1) >> 1; 188 if( f > 1 ) 189 { 190 factors[nf++] = f; 191 n = f == n ? 1 : n/f; 192 } 193 194 for( f = 3; n > 1; ) 195 { 196 int d = n/f; 197 if( d*f == n ) 198 { 199 factors[nf++] = f; 200 n = d; 201 } 202 else 203 { 204 f += 2; 205 if( f*f > n ) 206 break; 207 } 208 } 209 210 if( n > 1 ) 211 factors[nf++] = n; 212 213 f = (factors[0] & 1) == 0; 214 for( i = f; i < (nf+f)/2; i++ ) 215 CV_SWAP( factors[i], factors[nf-i-1+f], j ); 216 217 return nf; 218} 219 220static void 221icvDFTInit( int n0, int nf, int* factors, int* itab, int elem_size, void* _wave, int inv_itab ) 222{ 223 int digits[34], radix[34]; 224 int n = factors[0], m = 0; 225 int* itab0 = itab; 226 int i, j, k; 227 CvComplex64f w, w1; 228 double t; 229 230 if( n0 <= 5 ) 231 { 232 itab[0] = 0; 233 itab[n0-1] = n0-1; 234 235 if( n0 != 4 ) 236 { 237 for( i = 1; i < n0-1; i++ ) 238 itab[i] = i; 239 } 240 else 241 { 242 itab[1] = 2; 243 itab[2] = 1; 244 } 245 if( n0 == 5 ) 246 { 247 if( elem_size == sizeof(CvComplex64f) ) 248 ((CvComplex64f*)_wave)[0] = CvComplex64f(1.,0.); 249 else 250 ((CvComplex32f*)_wave)[0] = CvComplex32f(1.f,0.f); 251 } 252 if( n0 != 4 ) 253 return; 254 m = 2; 255 } 256 else 257 { 258 // radix[] is initialized from index 'nf' down to zero 259 assert (nf < 34); 260 radix[nf] = 1; 261 digits[nf] = 0; 262 for( i = 0; i < nf; i++ ) 263 { 264 digits[i] = 0; 265 radix[nf-i-1] = radix[nf-i]*factors[nf-i-1]; 266 } 267 268 if( inv_itab && factors[0] != factors[nf-1] ) 269 itab = (int*)_wave; 270 271 if( (n & 1) == 0 ) 272 { 273 int a = radix[1], na2 = n*a>>1, na4 = na2 >> 1; 274 m = icvlog2(n); 275 276 if( n <= 2 ) 277 { 278 itab[0] = 0; 279 itab[1] = na2; 280 } 281 else if( n <= 256 ) 282 { 283 int shift = 10 - m; 284 for( i = 0; i <= n - 4; i += 4 ) 285 { 286 j = (icvRevTable[i>>2]>>shift)*a; 287 itab[i] = j; 288 itab[i+1] = j + na2; 289 itab[i+2] = j + na4; 290 itab[i+3] = j + na2 + na4; 291 } 292 } 293 else 294 { 295 int shift = 34 - m; 296 for( i = 0; i < n; i += 4 ) 297 { 298 int i4 = i >> 2; 299 j = icvBitRev(i4,shift)*a; 300 itab[i] = j; 301 itab[i+1] = j + na2; 302 itab[i+2] = j + na4; 303 itab[i+3] = j + na2 + na4; 304 } 305 } 306 307 digits[1]++; 308 309 if( nf >= 2 ) 310 { 311 for( i = n, j = radix[2]; i < n0; ) 312 { 313 for( k = 0; k < n; k++ ) 314 itab[i+k] = itab[k] + j; 315 if( (i += n) >= n0 ) 316 break; 317 j += radix[2]; 318 for( k = 1; ++digits[k] >= factors[k]; k++ ) 319 { 320 digits[k] = 0; 321 j += radix[k+2] - radix[k]; 322 } 323 } 324 } 325 } 326 else 327 { 328 for( i = 0, j = 0;; ) 329 { 330 itab[i] = j; 331 if( ++i >= n0 ) 332 break; 333 j += radix[1]; 334 for( k = 0; ++digits[k] >= factors[k]; k++ ) 335 { 336 digits[k] = 0; 337 j += radix[k+2] - radix[k]; 338 } 339 } 340 } 341 342 if( itab != itab0 ) 343 { 344 itab0[0] = 0; 345 for( i = n0 & 1; i < n0; i += 2 ) 346 { 347 int k0 = itab[i]; 348 int k1 = itab[i+1]; 349 itab0[k0] = i; 350 itab0[k1] = i+1; 351 } 352 } 353 } 354 355 if( (n0 & (n0-1)) == 0 ) 356 { 357 w.re = w1.re = icvDxtTab[m][0]; 358 w.im = w1.im = -icvDxtTab[m][1]; 359 } 360 else 361 { 362 t = -CV_PI*2/n0; 363 w.im = w1.im = sin(t); 364 w.re = w1.re = sqrt(1. - w1.im*w1.im); 365 } 366 n = (n0+1)/2; 367 368 if( elem_size == sizeof(CvComplex64f) ) 369 { 370 CvComplex64f* wave = (CvComplex64f*)_wave; 371 372 wave[0].re = 1.; 373 wave[0].im = 0.; 374 375 if( (n0 & 1) == 0 ) 376 { 377 wave[n].re = -1.; 378 wave[n].im = 0; 379 } 380 381 for( i = 1; i < n; i++ ) 382 { 383 wave[i] = w; 384 wave[n0-i].re = w.re; 385 wave[n0-i].im = -w.im; 386 387 t = w.re*w1.re - w.im*w1.im; 388 w.im = w.re*w1.im + w.im*w1.re; 389 w.re = t; 390 } 391 } 392 else 393 { 394 CvComplex32f* wave = (CvComplex32f*)_wave; 395 assert( elem_size == sizeof(CvComplex32f) ); 396 397 wave[0].re = 1.f; 398 wave[0].im = 0.f; 399 400 if( (n0 & 1) == 0 ) 401 { 402 wave[n].re = -1.f; 403 wave[n].im = 0.f; 404 } 405 406 for( i = 1; i < n; i++ ) 407 { 408 wave[i].re = (float)w.re; 409 wave[i].im = (float)w.im; 410 wave[n0-i].re = (float)w.re; 411 wave[n0-i].im = (float)-w.im; 412 413 t = w.re*w1.re - w.im*w1.im; 414 w.im = w.re*w1.im + w.im*w1.re; 415 w.re = t; 416 } 417 } 418} 419 420 421static const double icv_sin_120 = 0.86602540378443864676372317075294; 422static const double icv_sin_45 = 0.70710678118654752440084436210485; 423static const double icv_fft5_2 = 0.559016994374947424102293417182819; 424static const double icv_fft5_3 = -0.951056516295153572116439333379382; 425static const double icv_fft5_4 = -1.538841768587626701285145288018455; 426static const double icv_fft5_5 = 0.363271264002680442947733378740309; 427 428#define ICV_DFT_NO_PERMUTE 2 429#define ICV_DFT_COMPLEX_INPUT_OR_OUTPUT 4 430 431// mixed-radix complex discrete Fourier transform: double-precision version 432static CvStatus CV_STDCALL 433icvDFT_64fc( const CvComplex64f* src, CvComplex64f* dst, int n, 434 int nf, int* factors, const int* itab, 435 const CvComplex64f* wave, int tab_size, 436 const void* spec, CvComplex64f* buf, 437 int flags, double scale ) 438{ 439 int n0 = n, f_idx, nx; 440 int inv = flags & CV_DXT_INVERSE; 441 int dw0 = tab_size, dw; 442 int i, j, k; 443 CvComplex64f t; 444 int tab_step; 445 446 if( spec ) 447 { 448 assert( icvDFTFwd_CToC_64fc_p != 0 && icvDFTInv_CToC_64fc_p != 0 ); 449 return !inv ? 450 icvDFTFwd_CToC_64fc_p( src, dst, spec, buf ): 451 icvDFTInv_CToC_64fc_p( src, dst, spec, buf ); 452 } 453 454 tab_step = tab_size == n ? 1 : tab_size == n*2 ? 2 : tab_size/n; 455 456 // 0. shuffle data 457 if( dst != src ) 458 { 459 assert( (flags & ICV_DFT_NO_PERMUTE) == 0 ); 460 if( !inv ) 461 { 462 for( i = 0; i <= n - 2; i += 2, itab += 2*tab_step ) 463 { 464 int k0 = itab[0], k1 = itab[tab_step]; 465 assert( (unsigned)k0 < (unsigned)n && (unsigned)k1 < (unsigned)n ); 466 dst[i] = src[k0]; dst[i+1] = src[k1]; 467 } 468 469 if( i < n ) 470 dst[n-1] = src[n-1]; 471 } 472 else 473 { 474 for( i = 0; i <= n - 2; i += 2, itab += 2*tab_step ) 475 { 476 int k0 = itab[0], k1 = itab[tab_step]; 477 assert( (unsigned)k0 < (unsigned)n && (unsigned)k1 < (unsigned)n ); 478 t.re = src[k0].re; t.im = -src[k0].im; 479 dst[i] = t; 480 t.re = src[k1].re; t.im = -src[k1].im; 481 dst[i+1] = t; 482 } 483 484 if( i < n ) 485 { 486 t.re = src[n-1].re; t.im = -src[n-1].im; 487 dst[i] = t; 488 } 489 } 490 } 491 else 492 { 493 if( (flags & ICV_DFT_NO_PERMUTE) == 0 ) 494 { 495 if( factors[0] != factors[nf-1] ) 496 return CV_INPLACE_NOT_SUPPORTED_ERR; 497 if( nf == 1 ) 498 { 499 if( (n & 3) == 0 ) 500 { 501 int n2 = n/2; 502 CvComplex64f* dsth = dst + n2; 503 504 for( i = 0; i < n2; i += 2, itab += tab_step*2 ) 505 { 506 j = itab[0]; 507 assert( (unsigned)j < (unsigned)n2 ); 508 509 CV_SWAP(dst[i+1], dsth[j], t); 510 if( j > i ) 511 { 512 CV_SWAP(dst[i], dst[j], t); 513 CV_SWAP(dsth[i+1], dsth[j+1], t); 514 } 515 } 516 } 517 // else do nothing 518 } 519 else 520 { 521 for( i = 0; i < n; i++, itab += tab_step ) 522 { 523 j = itab[0]; 524 assert( (unsigned)j < (unsigned)n ); 525 if( j > i ) 526 CV_SWAP(dst[i], dst[j], t); 527 } 528 } 529 } 530 531 if( inv ) 532 { 533 for( i = 0; i <= n - 2; i += 2 ) 534 { 535 double t0 = -dst[i].im; 536 double t1 = -dst[i+1].im; 537 dst[i].im = t0; dst[i+1].im = t1; 538 } 539 540 if( i < n ) 541 dst[n-1].im = -dst[n-1].im; 542 } 543 } 544 545 n = 1; 546 // 1. power-2 transforms 547 if( (factors[0] & 1) == 0 ) 548 { 549 // radix-4 transform 550 for( ; n*4 <= factors[0]; ) 551 { 552 nx = n; 553 n *= 4; 554 dw0 /= 4; 555 556 for( i = 0; i < n0; i += n ) 557 { 558 CvComplex64f* v0; 559 CvComplex64f* v1; 560 double r0, i0, r1, i1, r2, i2, r3, i3, r4, i4; 561 562 v0 = dst + i; 563 v1 = v0 + nx*2; 564 565 r2 = v0[0].re; i2 = v0[0].im; 566 r1 = v0[nx].re; i1 = v0[nx].im; 567 568 r0 = r1 + r2; i0 = i1 + i2; 569 r2 -= r1; i2 -= i1; 570 571 i3 = v1[nx].re; r3 = v1[nx].im; 572 i4 = v1[0].re; r4 = v1[0].im; 573 574 r1 = i4 + i3; i1 = r4 + r3; 575 r3 = r4 - r3; i3 = i3 - i4; 576 577 v0[0].re = r0 + r1; v0[0].im = i0 + i1; 578 v1[0].re = r0 - r1; v1[0].im = i0 - i1; 579 v0[nx].re = r2 + r3; v0[nx].im = i2 + i3; 580 v1[nx].re = r2 - r3; v1[nx].im = i2 - i3; 581 582 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 ) 583 { 584 v0 = dst + i + j; 585 v1 = v0 + nx*2; 586 587 r2 = v0[nx].re*wave[dw*2].re - v0[nx].im*wave[dw*2].im; 588 i2 = v0[nx].re*wave[dw*2].im + v0[nx].im*wave[dw*2].re; 589 r0 = v1[0].re*wave[dw].im + v1[0].im*wave[dw].re; 590 i0 = v1[0].re*wave[dw].re - v1[0].im*wave[dw].im; 591 r3 = v1[nx].re*wave[dw*3].im + v1[nx].im*wave[dw*3].re; 592 i3 = v1[nx].re*wave[dw*3].re - v1[nx].im*wave[dw*3].im; 593 594 r1 = i0 + i3; i1 = r0 + r3; 595 r3 = r0 - r3; i3 = i3 - i0; 596 r4 = v0[0].re; i4 = v0[0].im; 597 598 r0 = r4 + r2; i0 = i4 + i2; 599 r2 = r4 - r2; i2 = i4 - i2; 600 601 v0[0].re = r0 + r1; v0[0].im = i0 + i1; 602 v1[0].re = r0 - r1; v1[0].im = i0 - i1; 603 v0[nx].re = r2 + r3; v0[nx].im = i2 + i3; 604 v1[nx].re = r2 - r3; v1[nx].im = i2 - i3; 605 } 606 } 607 } 608 609 for( ; n < factors[0]; ) 610 { 611 // do the remaining radix-2 transform 612 nx = n; 613 n *= 2; 614 dw0 /= 2; 615 616 for( i = 0; i < n0; i += n ) 617 { 618 CvComplex64f* v = dst + i; 619 double r0 = v[0].re + v[nx].re; 620 double i0 = v[0].im + v[nx].im; 621 double r1 = v[0].re - v[nx].re; 622 double i1 = v[0].im - v[nx].im; 623 v[0].re = r0; v[0].im = i0; 624 v[nx].re = r1; v[nx].im = i1; 625 626 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 ) 627 { 628 v = dst + i + j; 629 r1 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im; 630 i1 = v[nx].im*wave[dw].re + v[nx].re*wave[dw].im; 631 r0 = v[0].re; i0 = v[0].im; 632 633 v[0].re = r0 + r1; v[0].im = i0 + i1; 634 v[nx].re = r0 - r1; v[nx].im = i0 - i1; 635 } 636 } 637 } 638 } 639 640 // 2. all the other transforms 641 for( f_idx = (factors[0]&1) ? 0 : 1; f_idx < nf; f_idx++ ) 642 { 643 int factor = factors[f_idx]; 644 nx = n; 645 n *= factor; 646 dw0 /= factor; 647 648 if( factor == 3 ) 649 { 650 // radix-3 651 for( i = 0; i < n0; i += n ) 652 { 653 CvComplex64f* v = dst + i; 654 655 double r1 = v[nx].re + v[nx*2].re; 656 double i1 = v[nx].im + v[nx*2].im; 657 double r0 = v[0].re; 658 double i0 = v[0].im; 659 double r2 = icv_sin_120*(v[nx].im - v[nx*2].im); 660 double i2 = icv_sin_120*(v[nx*2].re - v[nx].re); 661 v[0].re = r0 + r1; v[0].im = i0 + i1; 662 r0 -= 0.5*r1; i0 -= 0.5*i1; 663 v[nx].re = r0 + r2; v[nx].im = i0 + i2; 664 v[nx*2].re = r0 - r2; v[nx*2].im = i0 - i2; 665 666 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 ) 667 { 668 v = dst + i + j; 669 r0 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im; 670 i0 = v[nx].re*wave[dw].im + v[nx].im*wave[dw].re; 671 i2 = v[nx*2].re*wave[dw*2].re - v[nx*2].im*wave[dw*2].im; 672 r2 = v[nx*2].re*wave[dw*2].im + v[nx*2].im*wave[dw*2].re; 673 r1 = r0 + i2; i1 = i0 + r2; 674 675 r2 = icv_sin_120*(i0 - r2); i2 = icv_sin_120*(i2 - r0); 676 r0 = v[0].re; i0 = v[0].im; 677 v[0].re = r0 + r1; v[0].im = i0 + i1; 678 r0 -= 0.5*r1; i0 -= 0.5*i1; 679 v[nx].re = r0 + r2; v[nx].im = i0 + i2; 680 v[nx*2].re = r0 - r2; v[nx*2].im = i0 - i2; 681 } 682 } 683 } 684 else if( factor == 5 ) 685 { 686 // radix-5 687 for( i = 0; i < n0; i += n ) 688 { 689 for( j = 0, dw = 0; j < nx; j++, dw += dw0 ) 690 { 691 CvComplex64f* v0 = dst + i + j; 692 CvComplex64f* v1 = v0 + nx*2; 693 CvComplex64f* v2 = v1 + nx*2; 694 695 double r0, i0, r1, i1, r2, i2, r3, i3, r4, i4, r5, i5; 696 697 r3 = v0[nx].re*wave[dw].re - v0[nx].im*wave[dw].im; 698 i3 = v0[nx].re*wave[dw].im + v0[nx].im*wave[dw].re; 699 r2 = v2[0].re*wave[dw*4].re - v2[0].im*wave[dw*4].im; 700 i2 = v2[0].re*wave[dw*4].im + v2[0].im*wave[dw*4].re; 701 702 r1 = r3 + r2; i1 = i3 + i2; 703 r3 -= r2; i3 -= i2; 704 705 r4 = v1[nx].re*wave[dw*3].re - v1[nx].im*wave[dw*3].im; 706 i4 = v1[nx].re*wave[dw*3].im + v1[nx].im*wave[dw*3].re; 707 r0 = v1[0].re*wave[dw*2].re - v1[0].im*wave[dw*2].im; 708 i0 = v1[0].re*wave[dw*2].im + v1[0].im*wave[dw*2].re; 709 710 r2 = r4 + r0; i2 = i4 + i0; 711 r4 -= r0; i4 -= i0; 712 713 r0 = v0[0].re; i0 = v0[0].im; 714 r5 = r1 + r2; i5 = i1 + i2; 715 716 v0[0].re = r0 + r5; v0[0].im = i0 + i5; 717 718 r0 -= 0.25*r5; i0 -= 0.25*i5; 719 r1 = icv_fft5_2*(r1 - r2); i1 = icv_fft5_2*(i1 - i2); 720 r2 = -icv_fft5_3*(i3 + i4); i2 = icv_fft5_3*(r3 + r4); 721 722 i3 *= -icv_fft5_5; r3 *= icv_fft5_5; 723 i4 *= -icv_fft5_4; r4 *= icv_fft5_4; 724 725 r5 = r2 + i3; i5 = i2 + r3; 726 r2 -= i4; i2 -= r4; 727 728 r3 = r0 + r1; i3 = i0 + i1; 729 r0 -= r1; i0 -= i1; 730 731 v0[nx].re = r3 + r2; v0[nx].im = i3 + i2; 732 v2[0].re = r3 - r2; v2[0].im = i3 - i2; 733 734 v1[0].re = r0 + r5; v1[0].im = i0 + i5; 735 v1[nx].re = r0 - r5; v1[nx].im = i0 - i5; 736 } 737 } 738 } 739 else 740 { 741 // radix-"factor" - an odd number 742 int p, q, factor2 = (factor - 1)/2; 743 int d, dd, dw_f = tab_size/factor; 744 CvComplex64f* a = buf; 745 CvComplex64f* b = buf + factor2; 746 747 for( i = 0; i < n0; i += n ) 748 { 749 for( j = 0, dw = 0; j < nx; j++, dw += dw0 ) 750 { 751 CvComplex64f* v = dst + i + j; 752 CvComplex64f v_0 = v[0]; 753 CvComplex64f vn_0 = v_0; 754 755 if( j == 0 ) 756 { 757 for( p = 1, k = nx; p <= factor2; p++, k += nx ) 758 { 759 double r0 = v[k].re + v[n-k].re; 760 double i0 = v[k].im - v[n-k].im; 761 double r1 = v[k].re - v[n-k].re; 762 double i1 = v[k].im + v[n-k].im; 763 764 vn_0.re += r0; vn_0.im += i1; 765 a[p-1].re = r0; a[p-1].im = i0; 766 b[p-1].re = r1; b[p-1].im = i1; 767 } 768 } 769 else 770 { 771 const CvComplex64f* wave_ = wave + dw*factor; 772 d = dw; 773 774 for( p = 1, k = nx; p <= factor2; p++, k += nx, d += dw ) 775 { 776 double r2 = v[k].re*wave[d].re - v[k].im*wave[d].im; 777 double i2 = v[k].re*wave[d].im + v[k].im*wave[d].re; 778 779 double r1 = v[n-k].re*wave_[-d].re - v[n-k].im*wave_[-d].im; 780 double i1 = v[n-k].re*wave_[-d].im + v[n-k].im*wave_[-d].re; 781 782 double r0 = r2 + r1; 783 double i0 = i2 - i1; 784 r1 = r2 - r1; 785 i1 = i2 + i1; 786 787 vn_0.re += r0; vn_0.im += i1; 788 a[p-1].re = r0; a[p-1].im = i0; 789 b[p-1].re = r1; b[p-1].im = i1; 790 } 791 } 792 793 v[0] = vn_0; 794 795 for( p = 1, k = nx; p <= factor2; p++, k += nx ) 796 { 797 CvComplex64f s0 = v_0, s1 = v_0; 798 d = dd = dw_f*p; 799 800 for( q = 0; q < factor2; q++ ) 801 { 802 double r0 = wave[d].re * a[q].re; 803 double i0 = wave[d].im * a[q].im; 804 double r1 = wave[d].re * b[q].im; 805 double i1 = wave[d].im * b[q].re; 806 807 s1.re += r0 + i0; s0.re += r0 - i0; 808 s1.im += r1 - i1; s0.im += r1 + i1; 809 810 d += dd; 811 d -= -(d >= tab_size) & tab_size; 812 } 813 814 v[k] = s0; 815 v[n-k] = s1; 816 } 817 } 818 } 819 } 820 } 821 822 if( fabs(scale - 1.) > DBL_EPSILON ) 823 { 824 double re_scale = scale, im_scale = scale; 825 if( inv ) 826 im_scale = -im_scale; 827 828 for( i = 0; i < n0; i++ ) 829 { 830 double t0 = dst[i].re*re_scale; 831 double t1 = dst[i].im*im_scale; 832 dst[i].re = t0; 833 dst[i].im = t1; 834 } 835 } 836 else if( inv ) 837 { 838 for( i = 0; i <= n0 - 2; i += 2 ) 839 { 840 double t0 = -dst[i].im; 841 double t1 = -dst[i+1].im; 842 dst[i].im = t0; 843 dst[i+1].im = t1; 844 } 845 846 if( i < n0 ) 847 dst[n0-1].im = -dst[n0-1].im; 848 } 849 850 return CV_OK; 851} 852 853 854// mixed-radix complex discrete Fourier transform: single-precision version 855static CvStatus CV_STDCALL 856icvDFT_32fc( const CvComplex32f* src, CvComplex32f* dst, int n, 857 int nf, int* factors, const int* itab, 858 const CvComplex32f* wave, int tab_size, 859 const void* spec, CvComplex32f* buf, 860 int flags, double scale ) 861{ 862 int n0 = n, f_idx, nx; 863 int inv = flags & CV_DXT_INVERSE; 864 int dw0 = tab_size, dw; 865 int i, j, k; 866 CvComplex32f t; 867 int tab_step = tab_size == n ? 1 : tab_size == n*2 ? 2 : tab_size/n; 868 869 if( spec ) 870 { 871 assert( icvDFTFwd_CToC_32fc_p != 0 && icvDFTInv_CToC_32fc_p != 0 ); 872 return !inv ? 873 icvDFTFwd_CToC_32fc_p( src, dst, spec, buf ): 874 icvDFTInv_CToC_32fc_p( src, dst, spec, buf ); 875 } 876 877 // 0. shuffle data 878 if( dst != src ) 879 { 880 assert( (flags & ICV_DFT_NO_PERMUTE) == 0 ); 881 if( !inv ) 882 { 883 for( i = 0; i <= n - 2; i += 2, itab += 2*tab_step ) 884 { 885 int k0 = itab[0], k1 = itab[tab_step]; 886 assert( (unsigned)k0 < (unsigned)n && (unsigned)k1 < (unsigned)n ); 887 dst[i] = src[k0]; dst[i+1] = src[k1]; 888 } 889 890 if( i < n ) 891 dst[n-1] = src[n-1]; 892 } 893 else 894 { 895 for( i = 0; i <= n - 2; i += 2, itab += 2*tab_step ) 896 { 897 int k0 = itab[0], k1 = itab[tab_step]; 898 assert( (unsigned)k0 < (unsigned)n && (unsigned)k1 < (unsigned)n ); 899 t.re = src[k0].re; t.im = -src[k0].im; 900 dst[i] = t; 901 t.re = src[k1].re; t.im = -src[k1].im; 902 dst[i+1] = t; 903 } 904 905 if( i < n ) 906 { 907 t.re = src[n-1].re; t.im = -src[n-1].im; 908 dst[i] = t; 909 } 910 } 911 } 912 else 913 { 914 if( (flags & ICV_DFT_NO_PERMUTE) == 0 ) 915 { 916 if( factors[0] != factors[nf-1] ) 917 return CV_INPLACE_NOT_SUPPORTED_ERR; 918 if( nf == 1 ) 919 { 920 if( (n & 3) == 0 ) 921 { 922 int n2 = n/2; 923 CvComplex32f* dsth = dst + n2; 924 925 for( i = 0; i < n2; i += 2, itab += tab_step*2 ) 926 { 927 j = itab[0]; 928 assert( (unsigned)j < (unsigned)n2 ); 929 930 CV_SWAP(dst[i+1], dsth[j], t); 931 if( j > i ) 932 { 933 CV_SWAP(dst[i], dst[j], t); 934 CV_SWAP(dsth[i+1], dsth[j+1], t); 935 } 936 } 937 } 938 // else do nothing 939 } 940 else 941 { 942 for( i = 0; 943 i < n; 944 i++) 945 { 946 j = itab[0]; 947 assert( (unsigned)j < (unsigned)n ); 948 if( j > i ) 949 CV_SWAP(dst[i], dst[j], t); 950 itab += tab_step; 951 } 952 } 953 } 954 955 if( inv ) 956 { 957 // conjugate the vector - i.e. invert sign of the imaginary part 958 int* idst = (int*)dst; 959 for( i = 0; i <= n - 2; i += 2 ) 960 { 961 int t0 = idst[i*2+1] ^ 0x80000000; 962 int t1 = idst[i*2+3] ^ 0x80000000; 963 idst[i*2+1] = t0; idst[i*2+3] = t1; 964 } 965 966 if( i < n ) 967 idst[2*i+1] ^= 0x80000000; 968 } 969 } 970 971 n = 1; 972 // 1. power-2 transforms 973 if( (factors[0] & 1) == 0 ) 974 { 975 // radix-4 transform 976 for( ; n*4 <= factors[0]; ) 977 { 978 nx = n; 979 n *= 4; 980 dw0 /= 4; 981 982 for( i = 0; i < n0; i += n ) 983 { 984 CvComplex32f* v0; 985 CvComplex32f* v1; 986 double r0, i0, r1, i1, r2, i2, r3, i3, r4, i4; 987 988 v0 = dst + i; 989 v1 = v0 + nx*2; 990 991 r2 = v0[0].re; i2 = v0[0].im; 992 r1 = v0[nx].re; i1 = v0[nx].im; 993 994 r0 = r1 + r2; i0 = i1 + i2; 995 r2 -= r1; i2 -= i1; 996 997 i3 = v1[nx].re; r3 = v1[nx].im; 998 i4 = v1[0].re; r4 = v1[0].im; 999 1000 r1 = i4 + i3; i1 = r4 + r3; 1001 r3 = r4 - r3; i3 = i3 - i4; 1002 1003 v0[0].re = (float)(r0 + r1); v0[0].im = (float)(i0 + i1); 1004 v1[0].re = (float)(r0 - r1); v1[0].im = (float)(i0 - i1); 1005 v0[nx].re = (float)(r2 + r3); v0[nx].im = (float)(i2 + i3); 1006 v1[nx].re = (float)(r2 - r3); v1[nx].im = (float)(i2 - i3); 1007 1008 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 ) 1009 { 1010 v0 = dst + i + j; 1011 v1 = v0 + nx*2; 1012 1013 r2 = v0[nx].re*wave[dw*2].re - v0[nx].im*wave[dw*2].im; 1014 i2 = v0[nx].re*wave[dw*2].im + v0[nx].im*wave[dw*2].re; 1015 r0 = v1[0].re*wave[dw].im + v1[0].im*wave[dw].re; 1016 i0 = v1[0].re*wave[dw].re - v1[0].im*wave[dw].im; 1017 r3 = v1[nx].re*wave[dw*3].im + v1[nx].im*wave[dw*3].re; 1018 i3 = v1[nx].re*wave[dw*3].re - v1[nx].im*wave[dw*3].im; 1019 1020 r1 = i0 + i3; i1 = r0 + r3; 1021 r3 = r0 - r3; i3 = i3 - i0; 1022 r4 = v0[0].re; i4 = v0[0].im; 1023 1024 r0 = r4 + r2; i0 = i4 + i2; 1025 r2 = r4 - r2; i2 = i4 - i2; 1026 1027 v0[0].re = (float)(r0 + r1); v0[0].im = (float)(i0 + i1); 1028 v1[0].re = (float)(r0 - r1); v1[0].im = (float)(i0 - i1); 1029 v0[nx].re = (float)(r2 + r3); v0[nx].im = (float)(i2 + i3); 1030 v1[nx].re = (float)(r2 - r3); v1[nx].im = (float)(i2 - i3); 1031 } 1032 } 1033 } 1034 1035 for( ; n < factors[0]; ) 1036 { 1037 // do the remaining radix-2 transform 1038 nx = n; 1039 n *= 2; 1040 dw0 /= 2; 1041 1042 for( i = 0; i < n0; i += n ) 1043 { 1044 CvComplex32f* v = dst + i; 1045 double r0 = v[0].re + v[nx].re; 1046 double i0 = v[0].im + v[nx].im; 1047 double r1 = v[0].re - v[nx].re; 1048 double i1 = v[0].im - v[nx].im; 1049 v[0].re = (float)r0; v[0].im = (float)i0; 1050 v[nx].re = (float)r1; v[nx].im = (float)i1; 1051 1052 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 ) 1053 { 1054 v = dst + i + j; 1055 r1 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im; 1056 i1 = v[nx].im*wave[dw].re + v[nx].re*wave[dw].im; 1057 r0 = v[0].re; i0 = v[0].im; 1058 1059 v[0].re = (float)(r0 + r1); v[0].im = (float)(i0 + i1); 1060 v[nx].re = (float)(r0 - r1); v[nx].im = (float)(i0 - i1); 1061 } 1062 } 1063 } 1064 } 1065 1066 // 2. all the other transforms 1067 for( f_idx = (factors[0]&1) ? 0 : 1; f_idx < nf; f_idx++ ) 1068 { 1069 int factor = factors[f_idx]; 1070 nx = n; 1071 n *= factor; 1072 dw0 /= factor; 1073 1074 if( factor == 3 ) 1075 { 1076 // radix-3 1077 for( i = 0; i < n0; i += n ) 1078 { 1079 CvComplex32f* v = dst + i; 1080 1081 double r1 = v[nx].re + v[nx*2].re; 1082 double i1 = v[nx].im + v[nx*2].im; 1083 double r0 = v[0].re; 1084 double i0 = v[0].im; 1085 double r2 = icv_sin_120*(v[nx].im - v[nx*2].im); 1086 double i2 = icv_sin_120*(v[nx*2].re - v[nx].re); 1087 v[0].re = (float)(r0 + r1); v[0].im = (float)(i0 + i1); 1088 r0 -= 0.5*r1; i0 -= 0.5*i1; 1089 v[nx].re = (float)(r0 + r2); v[nx].im = (float)(i0 + i2); 1090 v[nx*2].re = (float)(r0 - r2); v[nx*2].im = (float)(i0 - i2); 1091 1092 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 ) 1093 { 1094 v = dst + i + j; 1095 r0 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im; 1096 i0 = v[nx].re*wave[dw].im + v[nx].im*wave[dw].re; 1097 i2 = v[nx*2].re*wave[dw*2].re - v[nx*2].im*wave[dw*2].im; 1098 r2 = v[nx*2].re*wave[dw*2].im + v[nx*2].im*wave[dw*2].re; 1099 r1 = r0 + i2; i1 = i0 + r2; 1100 1101 r2 = icv_sin_120*(i0 - r2); i2 = icv_sin_120*(i2 - r0); 1102 r0 = v[0].re; i0 = v[0].im; 1103 v[0].re = (float)(r0 + r1); v[0].im = (float)(i0 + i1); 1104 r0 -= 0.5*r1; i0 -= 0.5*i1; 1105 v[nx].re = (float)(r0 + r2); v[nx].im = (float)(i0 + i2); 1106 v[nx*2].re = (float)(r0 - r2); v[nx*2].im = (float)(i0 - i2); 1107 } 1108 } 1109 } 1110 else if( factor == 5 ) 1111 { 1112 // radix-5 1113 for( i = 0; i < n0; i += n ) 1114 { 1115 for( j = 0, dw = 0; j < nx; j++, dw += dw0 ) 1116 { 1117 CvComplex32f* v0 = dst + i + j; 1118 CvComplex32f* v1 = v0 + nx*2; 1119 CvComplex32f* v2 = v1 + nx*2; 1120 1121 double r0, i0, r1, i1, r2, i2, r3, i3, r4, i4, r5, i5; 1122 1123 r3 = v0[nx].re*wave[dw].re - v0[nx].im*wave[dw].im; 1124 i3 = v0[nx].re*wave[dw].im + v0[nx].im*wave[dw].re; 1125 r2 = v2[0].re*wave[dw*4].re - v2[0].im*wave[dw*4].im; 1126 i2 = v2[0].re*wave[dw*4].im + v2[0].im*wave[dw*4].re; 1127 1128 r1 = r3 + r2; i1 = i3 + i2; 1129 r3 -= r2; i3 -= i2; 1130 1131 r4 = v1[nx].re*wave[dw*3].re - v1[nx].im*wave[dw*3].im; 1132 i4 = v1[nx].re*wave[dw*3].im + v1[nx].im*wave[dw*3].re; 1133 r0 = v1[0].re*wave[dw*2].re - v1[0].im*wave[dw*2].im; 1134 i0 = v1[0].re*wave[dw*2].im + v1[0].im*wave[dw*2].re; 1135 1136 r2 = r4 + r0; i2 = i4 + i0; 1137 r4 -= r0; i4 -= i0; 1138 1139 r0 = v0[0].re; i0 = v0[0].im; 1140 r5 = r1 + r2; i5 = i1 + i2; 1141 1142 v0[0].re = (float)(r0 + r5); v0[0].im = (float)(i0 + i5); 1143 1144 r0 -= 0.25*r5; i0 -= 0.25*i5; 1145 r1 = icv_fft5_2*(r1 - r2); i1 = icv_fft5_2*(i1 - i2); 1146 r2 = -icv_fft5_3*(i3 + i4); i2 = icv_fft5_3*(r3 + r4); 1147 1148 i3 *= -icv_fft5_5; r3 *= icv_fft5_5; 1149 i4 *= -icv_fft5_4; r4 *= icv_fft5_4; 1150 1151 r5 = r2 + i3; i5 = i2 + r3; 1152 r2 -= i4; i2 -= r4; 1153 1154 r3 = r0 + r1; i3 = i0 + i1; 1155 r0 -= r1; i0 -= i1; 1156 1157 v0[nx].re = (float)(r3 + r2); v0[nx].im = (float)(i3 + i2); 1158 v2[0].re = (float)(r3 - r2); v2[0].im = (float)(i3 - i2); 1159 1160 v1[0].re = (float)(r0 + r5); v1[0].im = (float)(i0 + i5); 1161 v1[nx].re = (float)(r0 - r5); v1[nx].im = (float)(i0 - i5); 1162 } 1163 } 1164 } 1165 else 1166 { 1167 // radix-"factor" - an odd number 1168 int p, q, factor2 = (factor - 1)/2; 1169 int d, dd, dw_f = tab_size/factor; 1170 CvComplex32f* a = buf; 1171 CvComplex32f* b = buf + factor2; 1172 1173 for( i = 0; i < n0; i += n ) 1174 { 1175 for( j = 0, dw = 0; j < nx; j++, dw += dw0 ) 1176 { 1177 CvComplex32f* v = dst + i + j; 1178 CvComplex32f v_0 = v[0]; 1179 CvComplex64f vn_0(v_0); 1180 1181 if( j == 0 ) 1182 { 1183 for( p = 1, k = nx; p <= factor2; p++, k += nx ) 1184 { 1185 double r0 = v[k].re + v[n-k].re; 1186 double i0 = v[k].im - v[n-k].im; 1187 double r1 = v[k].re - v[n-k].re; 1188 double i1 = v[k].im + v[n-k].im; 1189 1190 vn_0.re += r0; vn_0.im += i1; 1191 a[p-1].re = (float)r0; a[p-1].im = (float)i0; 1192 b[p-1].re = (float)r1; b[p-1].im = (float)i1; 1193 } 1194 } 1195 else 1196 { 1197 const CvComplex32f* wave_ = wave + dw*factor; 1198 d = dw; 1199 1200 for( p = 1, k = nx; p <= factor2; p++, k += nx, d += dw ) 1201 { 1202 double r2 = v[k].re*wave[d].re - v[k].im*wave[d].im; 1203 double i2 = v[k].re*wave[d].im + v[k].im*wave[d].re; 1204 1205 double r1 = v[n-k].re*wave_[-d].re - v[n-k].im*wave_[-d].im; 1206 double i1 = v[n-k].re*wave_[-d].im + v[n-k].im*wave_[-d].re; 1207 1208 double r0 = r2 + r1; 1209 double i0 = i2 - i1; 1210 r1 = r2 - r1; 1211 i1 = i2 + i1; 1212 1213 vn_0.re += r0; vn_0.im += i1; 1214 a[p-1].re = (float)r0; a[p-1].im = (float)i0; 1215 b[p-1].re = (float)r1; b[p-1].im = (float)i1; 1216 } 1217 } 1218 1219 v[0].re = (float)vn_0.re; 1220 v[0].im = (float)vn_0.im; 1221 1222 for( p = 1, k = nx; p <= factor2; p++, k += nx ) 1223 { 1224 CvComplex64f s0(v_0), s1 = s0; 1225 d = dd = dw_f*p; 1226 1227 for( q = 0; q < factor2; q++ ) 1228 { 1229 double r0 = wave[d].re * a[q].re; 1230 double i0 = wave[d].im * a[q].im; 1231 double r1 = wave[d].re * b[q].im; 1232 double i1 = wave[d].im * b[q].re; 1233 1234 s1.re += r0 + i0; s0.re += r0 - i0; 1235 s1.im += r1 - i1; s0.im += r1 + i1; 1236 1237 d += dd; 1238 d -= -(d >= tab_size) & tab_size; 1239 } 1240 1241 v[k].re = (float)s0.re; 1242 v[k].im = (float)s0.im; 1243 v[n-k].re = (float)s1.re; 1244 v[n-k].im = (float)s1.im; 1245 } 1246 } 1247 } 1248 } 1249 } 1250 1251 if( fabs(scale - 1.) > DBL_EPSILON ) 1252 { 1253 double re_scale = scale, im_scale = scale; 1254 if( inv ) 1255 im_scale = -im_scale; 1256 1257 for( i = 0; i < n0; i++ ) 1258 { 1259 double t0 = dst[i].re*re_scale; 1260 double t1 = dst[i].im*im_scale; 1261 dst[i].re = (float)t0; 1262 dst[i].im = (float)t1; 1263 } 1264 } 1265 else if( inv ) 1266 { 1267 for( i = 0; i <= n0 - 2; i += 2 ) 1268 { 1269 double t0 = -dst[i].im; 1270 double t1 = -dst[i+1].im; 1271 dst[i].im = (float)t0; 1272 dst[i+1].im = (float)t1; 1273 } 1274 1275 if( i < n0 ) 1276 dst[n0-1].im = -dst[n0-1].im; 1277 } 1278 1279 return CV_OK; 1280} 1281 1282 1283/* FFT of real vector 1284 output vector format: 1285 re(0), re(1), im(1), ... , re(n/2-1), im((n+1)/2-1) [, re((n+1)/2)] OR ... 1286 re(0), 0, re(1), im(1), ..., re(n/2-1), im((n+1)/2-1) [, re((n+1)/2), 0] */ 1287#define ICV_REAL_DFT( flavor, datatype ) \ 1288static CvStatus CV_STDCALL \ 1289icvRealDFT_##flavor( const datatype* src, datatype* dst, \ 1290 int n, int nf, int* factors, const int* itab, \ 1291 const CvComplex##flavor* wave, int tab_size, \ 1292 const void* spec, CvComplex##flavor* buf, \ 1293 int flags, double scale ) \ 1294{ \ 1295 int complex_output = (flags & ICV_DFT_COMPLEX_INPUT_OR_OUTPUT) != 0;\ 1296 int j, n2 = n >> 1; \ 1297 dst += complex_output; \ 1298 \ 1299 if( spec ) \ 1300 { \ 1301 icvDFTFwd_RToPack_##flavor##_p( src, dst, spec, buf ); \ 1302 goto finalize; \ 1303 } \ 1304 \ 1305 assert( tab_size == n ); \ 1306 \ 1307 if( n == 1 ) \ 1308 { \ 1309 dst[0] = (datatype)(src[0]*scale); \ 1310 } \ 1311 else if( n == 2 ) \ 1312 { \ 1313 double t = (src[0] + src[1])*scale; \ 1314 dst[1] = (datatype)((src[0] - src[1])*scale); \ 1315 dst[0] = (datatype)t; \ 1316 } \ 1317 else if( n & 1 ) \ 1318 { \ 1319 dst -= complex_output; \ 1320 CvComplex##flavor* _dst = (CvComplex##flavor*)dst; \ 1321 _dst[0].re = (datatype)(src[0]*scale); \ 1322 _dst[0].im = 0; \ 1323 for( j = 1; j < n; j += 2 ) \ 1324 { \ 1325 double t0 = src[itab[j]]*scale; \ 1326 double t1 = src[itab[j+1]]*scale; \ 1327 _dst[j].re = (datatype)t0; \ 1328 _dst[j].im = 0; \ 1329 _dst[j+1].re = (datatype)t1; \ 1330 _dst[j+1].im = 0; \ 1331 } \ 1332 icvDFT_##flavor##c( _dst, _dst, n, nf, factors, itab, wave, \ 1333 tab_size, 0, buf, ICV_DFT_NO_PERMUTE, 1. ); \ 1334 if( !complex_output ) \ 1335 dst[1] = dst[0]; \ 1336 return CV_OK; \ 1337 } \ 1338 else \ 1339 { \ 1340 double t0, t; \ 1341 double h1_re, h1_im, h2_re, h2_im; \ 1342 double scale2 = scale*0.5; \ 1343 factors[0] >>= 1; \ 1344 \ 1345 icvDFT_##flavor##c( (CvComplex##flavor*)src, \ 1346 (CvComplex##flavor*)dst, n2, \ 1347 nf - (factors[0] == 1), \ 1348 factors + (factors[0] == 1), \ 1349 itab, wave, tab_size, 0, buf, 0, 1. ); \ 1350 factors[0] <<= 1; \ 1351 \ 1352 t = dst[0] - dst[1]; \ 1353 dst[0] = (datatype)((dst[0] + dst[1])*scale); \ 1354 dst[1] = (datatype)(t*scale); \ 1355 \ 1356 t0 = dst[n2]; \ 1357 t = dst[n-1]; \ 1358 dst[n-1] = dst[1]; \ 1359 \ 1360 for( j = 2, wave++; j < n2; j += 2, wave++ ) \ 1361 { \ 1362 /* calc odd */ \ 1363 h2_re = scale2*(dst[j+1] + t); \ 1364 h2_im = scale2*(dst[n-j] - dst[j]); \ 1365 \ 1366 /* calc even */ \ 1367 h1_re = scale2*(dst[j] + dst[n-j]); \ 1368 h1_im = scale2*(dst[j+1] - t); \ 1369 \ 1370 /* rotate */ \ 1371 t = h2_re*wave->re - h2_im*wave->im; \ 1372 h2_im = h2_re*wave->im + h2_im*wave->re; \ 1373 h2_re = t; \ 1374 t = dst[n-j-1]; \ 1375 \ 1376 dst[j-1] = (datatype)(h1_re + h2_re); \ 1377 dst[n-j-1] = (datatype)(h1_re - h2_re); \ 1378 dst[j] = (datatype)(h1_im + h2_im); \ 1379 dst[n-j] = (datatype)(h2_im - h1_im); \ 1380 } \ 1381 \ 1382 if( j <= n2 ) \ 1383 { \ 1384 dst[n2-1] = (datatype)(t0*scale); \ 1385 dst[n2] = (datatype)(-t*scale); \ 1386 } \ 1387 } \ 1388 \ 1389finalize: \ 1390 if( complex_output ) \ 1391 { \ 1392 dst[-1] = dst[0]; \ 1393 dst[0] = 0; \ 1394 if( (n & 1) == 0 ) \ 1395 dst[n] = 0; \ 1396 } \ 1397 \ 1398 return CV_OK; \ 1399} 1400 1401 1402/* Inverse FFT of complex conjugate-symmetric vector 1403 input vector format: 1404 re[0], re[1], im[1], ... , re[n/2-1], im[n/2-1], re[n/2] OR 1405 re(0), 0, re(1), im(1), ..., re(n/2-1), im((n+1)/2-1) [, re((n+1)/2), 0] */ 1406#define ICV_CCS_IDFT( flavor, datatype ) \ 1407static CvStatus CV_STDCALL \ 1408icvCCSIDFT_##flavor( const datatype* src, datatype* dst, \ 1409 int n, int nf, int* factors, const int* itab, \ 1410 const CvComplex##flavor* wave, int tab_size, \ 1411 const void* spec, CvComplex##flavor* buf, \ 1412 int flags, double scale ) \ 1413{ \ 1414 int complex_input = (flags & ICV_DFT_COMPLEX_INPUT_OR_OUTPUT) != 0; \ 1415 int j, k, n2 = (n+1) >> 1; \ 1416 double save_s1 = 0.; \ 1417 double t0, t1, t2, t3, t; \ 1418 \ 1419 assert( tab_size == n ); \ 1420 \ 1421 if( complex_input ) \ 1422 { \ 1423 assert( src != dst ); \ 1424 save_s1 = src[1]; \ 1425 ((datatype*)src)[1] = src[0]; \ 1426 src++; \ 1427 } \ 1428 \ 1429 if( spec ) \ 1430 { \ 1431 icvDFTInv_PackToR_##flavor##_p( src, dst, spec, buf ); \ 1432 goto finalize; \ 1433 } \ 1434 \ 1435 if( n == 1 ) \ 1436 { \ 1437 dst[0] = (datatype)(src[0]*scale); \ 1438 } \ 1439 else if( n == 2 ) \ 1440 { \ 1441 t = (src[0] + src[1])*scale; \ 1442 dst[1] = (datatype)((src[0] - src[1])*scale); \ 1443 dst[0] = (datatype)t; \ 1444 } \ 1445 else if( n & 1 ) \ 1446 { \ 1447 CvComplex##flavor* _src = (CvComplex##flavor*)(src-1); \ 1448 CvComplex##flavor* _dst = (CvComplex##flavor*)dst; \ 1449 \ 1450 _dst[0].re = src[0]; \ 1451 _dst[0].im = 0; \ 1452 for( j = 1; j < n2; j++ ) \ 1453 { \ 1454 int k0 = itab[j], k1 = itab[n-j]; \ 1455 t0 = _src[j].re; t1 = _src[j].im; \ 1456 _dst[k0].re = (datatype)t0; _dst[k0].im = (datatype)-t1; \ 1457 _dst[k1].re = (datatype)t0; _dst[k1].im = (datatype)t1; \ 1458 } \ 1459 \ 1460 icvDFT_##flavor##c( _dst, _dst, n, nf, factors, itab, wave, \ 1461 tab_size, 0, buf, ICV_DFT_NO_PERMUTE, 1. ); \ 1462 dst[0] = (datatype)(dst[0]*scale); \ 1463 for( j = 1; j < n; j += 2 ) \ 1464 { \ 1465 t0 = dst[j*2]*scale; \ 1466 t1 = dst[j*2+2]*scale; \ 1467 dst[j] = (datatype)t0; \ 1468 dst[j+1] = (datatype)t1; \ 1469 } \ 1470 } \ 1471 else \ 1472 { \ 1473 int inplace = src == dst; \ 1474 const CvComplex##flavor* w = wave; \ 1475 \ 1476 t = src[1]; \ 1477 t0 = (src[0] + src[n-1]); \ 1478 t1 = (src[n-1] - src[0]); \ 1479 dst[0] = (datatype)t0; \ 1480 dst[1] = (datatype)t1; \ 1481 \ 1482 for( j = 2, w++; j < n2; j += 2, w++ ) \ 1483 { \ 1484 double h1_re, h1_im, h2_re, h2_im; \ 1485 \ 1486 h1_re = (t + src[n-j-1]); \ 1487 h1_im = (src[j] - src[n-j]); \ 1488 \ 1489 h2_re = (t - src[n-j-1]); \ 1490 h2_im = (src[j] + src[n-j]); \ 1491 \ 1492 t = h2_re*w->re + h2_im*w->im; \ 1493 h2_im = h2_im*w->re - h2_re*w->im; \ 1494 h2_re = t; \ 1495 \ 1496 t = src[j+1]; \ 1497 t0 = h1_re - h2_im; \ 1498 t1 = -h1_im - h2_re; \ 1499 t2 = h1_re + h2_im; \ 1500 t3 = h1_im - h2_re; \ 1501 \ 1502 if( inplace ) \ 1503 { \ 1504 dst[j] = (datatype)t0; \ 1505 dst[j+1] = (datatype)t1; \ 1506 dst[n-j] = (datatype)t2; \ 1507 dst[n-j+1]= (datatype)t3; \ 1508 } \ 1509 else \ 1510 { \ 1511 int j2 = j >> 1; \ 1512 k = itab[j2]; \ 1513 dst[k] = (datatype)t0; \ 1514 dst[k+1] = (datatype)t1; \ 1515 k = itab[n2-j2]; \ 1516 dst[k] = (datatype)t2; \ 1517 dst[k+1]= (datatype)t3; \ 1518 } \ 1519 } \ 1520 \ 1521 if( j <= n2 ) \ 1522 { \ 1523 t0 = t*2; \ 1524 t1 = src[n2]*2; \ 1525 \ 1526 if( inplace ) \ 1527 { \ 1528 dst[n2] = (datatype)t0; \ 1529 dst[n2+1] = (datatype)t1; \ 1530 } \ 1531 else \ 1532 { \ 1533 k = itab[n2]; \ 1534 dst[k*2] = (datatype)t0; \ 1535 dst[k*2+1] = (datatype)t1; \ 1536 } \ 1537 } \ 1538 \ 1539 factors[0] >>= 1; \ 1540 icvDFT_##flavor##c( (CvComplex##flavor*)dst, \ 1541 (CvComplex##flavor*)dst, n2, \ 1542 nf - (factors[0] == 1), \ 1543 factors + (factors[0] == 1), itab, \ 1544 wave, tab_size, 0, buf, \ 1545 inplace ? 0 : ICV_DFT_NO_PERMUTE, 1. ); \ 1546 factors[0] <<= 1; \ 1547 \ 1548 for( j = 0; j < n; j += 2 ) \ 1549 { \ 1550 t0 = dst[j]*scale; \ 1551 t1 = dst[j+1]*(-scale); \ 1552 dst[j] = (datatype)t0; \ 1553 dst[j+1] = (datatype)t1; \ 1554 } \ 1555 } \ 1556 \ 1557finalize: \ 1558 if( complex_input ) \ 1559 ((datatype*)src)[0] = (datatype)save_s1; \ 1560 \ 1561 return CV_OK; \ 1562} 1563 1564 1565ICV_REAL_DFT( 64f, double ) 1566ICV_CCS_IDFT( 64f, double ) 1567ICV_REAL_DFT( 32f, float ) 1568ICV_CCS_IDFT( 32f, float ) 1569 1570 1571static void 1572icvCopyColumn( const uchar* _src, int src_step, 1573 uchar* _dst, int dst_step, 1574 int len, int elem_size ) 1575{ 1576 int i, t0, t1; 1577 const int* src = (const int*)_src; 1578 int* dst = (int*)_dst; 1579 src_step /= sizeof(src[0]); 1580 dst_step /= sizeof(dst[0]); 1581 1582 if( elem_size == sizeof(int) ) 1583 { 1584 for( i = 0; i < len; i++, src += src_step, dst += dst_step ) 1585 dst[0] = src[0]; 1586 } 1587 else if( elem_size == sizeof(int)*2 ) 1588 { 1589 for( i = 0; i < len; i++, src += src_step, dst += dst_step ) 1590 { 1591 t0 = src[0]; t1 = src[1]; 1592 dst[0] = t0; dst[1] = t1; 1593 } 1594 } 1595 else if( elem_size == sizeof(int)*4 ) 1596 { 1597 for( i = 0; i < len; i++, src += src_step, dst += dst_step ) 1598 { 1599 t0 = src[0]; t1 = src[1]; 1600 dst[0] = t0; dst[1] = t1; 1601 t0 = src[2]; t1 = src[3]; 1602 dst[2] = t0; dst[3] = t1; 1603 } 1604 } 1605} 1606 1607 1608static void 1609icvCopyFrom2Columns( const uchar* _src, int src_step, 1610 uchar* _dst0, uchar* _dst1, 1611 int len, int elem_size ) 1612{ 1613 int i, t0, t1; 1614 const int* src = (const int*)_src; 1615 int* dst0 = (int*)_dst0; 1616 int* dst1 = (int*)_dst1; 1617 src_step /= sizeof(src[0]); 1618 1619 if( elem_size == sizeof(int) ) 1620 { 1621 for( i = 0; i < len; i++, src += src_step ) 1622 { 1623 t0 = src[0]; t1 = src[1]; 1624 dst0[i] = t0; dst1[i] = t1; 1625 } 1626 } 1627 else if( elem_size == sizeof(int)*2 ) 1628 { 1629 for( i = 0; i < len*2; i += 2, src += src_step ) 1630 { 1631 t0 = src[0]; t1 = src[1]; 1632 dst0[i] = t0; dst0[i+1] = t1; 1633 t0 = src[2]; t1 = src[3]; 1634 dst1[i] = t0; dst1[i+1] = t1; 1635 } 1636 } 1637 else if( elem_size == sizeof(int)*4 ) 1638 { 1639 for( i = 0; i < len*4; i += 4, src += src_step ) 1640 { 1641 t0 = src[0]; t1 = src[1]; 1642 dst0[i] = t0; dst0[i+1] = t1; 1643 t0 = src[2]; t1 = src[3]; 1644 dst0[i+2] = t0; dst0[i+3] = t1; 1645 t0 = src[4]; t1 = src[5]; 1646 dst1[i] = t0; dst1[i+1] = t1; 1647 t0 = src[6]; t1 = src[7]; 1648 dst1[i+2] = t0; dst1[i+3] = t1; 1649 } 1650 } 1651} 1652 1653 1654static void 1655icvCopyTo2Columns( const uchar* _src0, const uchar* _src1, 1656 uchar* _dst, int dst_step, 1657 int len, int elem_size ) 1658{ 1659 int i, t0, t1; 1660 const int* src0 = (const int*)_src0; 1661 const int* src1 = (const int*)_src1; 1662 int* dst = (int*)_dst; 1663 dst_step /= sizeof(dst[0]); 1664 1665 if( elem_size == sizeof(int) ) 1666 { 1667 for( i = 0; i < len; i++, dst += dst_step ) 1668 { 1669 t0 = src0[i]; t1 = src1[i]; 1670 dst[0] = t0; dst[1] = t1; 1671 } 1672 } 1673 else if( elem_size == sizeof(int)*2 ) 1674 { 1675 for( i = 0; i < len*2; i += 2, dst += dst_step ) 1676 { 1677 t0 = src0[i]; t1 = src0[i+1]; 1678 dst[0] = t0; dst[1] = t1; 1679 t0 = src1[i]; t1 = src1[i+1]; 1680 dst[2] = t0; dst[3] = t1; 1681 } 1682 } 1683 else if( elem_size == sizeof(int)*4 ) 1684 { 1685 for( i = 0; i < len*4; i += 4, dst += dst_step ) 1686 { 1687 t0 = src0[i]; t1 = src0[i+1]; 1688 dst[0] = t0; dst[1] = t1; 1689 t0 = src0[i+2]; t1 = src0[i+3]; 1690 dst[2] = t0; dst[3] = t1; 1691 t0 = src1[i]; t1 = src1[i+1]; 1692 dst[4] = t0; dst[5] = t1; 1693 t0 = src1[i+2]; t1 = src1[i+3]; 1694 dst[6] = t0; dst[7] = t1; 1695 } 1696 } 1697} 1698 1699 1700static void 1701icvExpandCCS( uchar* _ptr, int len, int elem_size ) 1702{ 1703 int i; 1704 _ptr -= elem_size; 1705 memcpy( _ptr, _ptr + elem_size, elem_size ); 1706 memset( _ptr + elem_size, 0, elem_size ); 1707 if( (len & 1) == 0 ) 1708 memset( _ptr + (len+1)*elem_size, 0, elem_size ); 1709 1710 if( elem_size == sizeof(float) ) 1711 { 1712 CvComplex32f* ptr = (CvComplex32f*)_ptr; 1713 1714 for( i = 1; i < (len+1)/2; i++ ) 1715 { 1716 CvComplex32f t; 1717 t.re = ptr[i].re; 1718 t.im = -ptr[i].im; 1719 ptr[len-i] = t; 1720 } 1721 } 1722 else 1723 { 1724 CvComplex64f* ptr = (CvComplex64f*)_ptr; 1725 1726 for( i = 1; i < (len+1)/2; i++ ) 1727 { 1728 CvComplex64f t; 1729 t.re = ptr[i].re; 1730 t.im = -ptr[i].im; 1731 ptr[len-i] = t; 1732 } 1733 } 1734} 1735 1736 1737typedef CvStatus (CV_STDCALL *CvDFTFunc)( 1738 const void* src, void* dst, int n, int nf, int* factors, 1739 const int* itab, const void* wave, int tab_size, 1740 const void* spec, void* buf, int inv, double scale ); 1741 1742CV_IMPL void 1743cvDFT( const CvArr* srcarr, CvArr* dstarr, int flags, int nonzero_rows ) 1744{ 1745 static CvDFTFunc dft_tbl[6]; 1746 static int inittab = 0; 1747 1748 void* buffer = 0; 1749 int local_alloc = 1; 1750 int depth = -1; 1751 void *spec_c = 0, *spec_r = 0, *spec = 0; 1752 1753 CV_FUNCNAME( "cvDFT" ); 1754 1755 __BEGIN__; 1756 1757 int prev_len = 0, buf_size = 0, stage = 0; 1758 int nf = 0, inv = (flags & CV_DXT_INVERSE) != 0; 1759 int real_transform = 0; 1760 CvMat *src = (CvMat*)srcarr, *dst = (CvMat*)dstarr; 1761 CvMat srcstub, dststub, *src0; 1762 int complex_elem_size, elem_size; 1763 int factors[34], inplace_transform = 0; 1764 int ipp_norm_flag = 0; 1765 1766 if( !inittab ) 1767 { 1768 dft_tbl[0] = (CvDFTFunc)icvDFT_32fc; 1769 dft_tbl[1] = (CvDFTFunc)icvRealDFT_32f; 1770 dft_tbl[2] = (CvDFTFunc)icvCCSIDFT_32f; 1771 dft_tbl[3] = (CvDFTFunc)icvDFT_64fc; 1772 dft_tbl[4] = (CvDFTFunc)icvRealDFT_64f; 1773 dft_tbl[5] = (CvDFTFunc)icvCCSIDFT_64f; 1774 inittab = 1; 1775 } 1776 1777 if( !CV_IS_MAT( src )) 1778 { 1779 int coi = 0; 1780 CV_CALL( src = cvGetMat( src, &srcstub, &coi )); 1781 1782 if( coi != 0 ) 1783 CV_ERROR( CV_BadCOI, "" ); 1784 } 1785 1786 if( !CV_IS_MAT( dst )) 1787 { 1788 int coi = 0; 1789 CV_CALL( dst = cvGetMat( dst, &dststub, &coi )); 1790 1791 if( coi != 0 ) 1792 CV_ERROR( CV_BadCOI, "" ); 1793 } 1794 1795 src0 = src; 1796 elem_size = CV_ELEM_SIZE1(src->type); 1797 complex_elem_size = elem_size*2; 1798 1799 // check types and sizes 1800 if( !CV_ARE_DEPTHS_EQ(src, dst) ) 1801 CV_ERROR( CV_StsUnmatchedFormats, "" ); 1802 1803 depth = CV_MAT_DEPTH(src->type); 1804 if( depth < CV_32F ) 1805 CV_ERROR( CV_StsUnsupportedFormat, 1806 "Only 32fC1, 32fC2, 64fC1 and 64fC2 formats are supported" ); 1807 1808 if( CV_ARE_CNS_EQ(src, dst) ) 1809 { 1810 if( CV_MAT_CN(src->type) > 2 ) 1811 CV_ERROR( CV_StsUnsupportedFormat, 1812 "Only 32fC1, 32fC2, 64fC1 and 64fC2 formats are supported" ); 1813 1814 if( !CV_ARE_SIZES_EQ(src, dst) ) 1815 CV_ERROR( CV_StsUnmatchedSizes, "" ); 1816 real_transform = CV_MAT_CN(src->type) == 1; 1817 if( !real_transform ) 1818 elem_size = complex_elem_size; 1819 } 1820 else if( !inv && CV_MAT_CN(src->type) == 1 && CV_MAT_CN(dst->type) == 2 ) 1821 { 1822 if( (src->cols != 1 || dst->cols != 1 || 1823 (src->rows/2+1 != dst->rows && src->rows != dst->rows)) && 1824 (src->cols/2+1 != dst->cols || src->rows != dst->rows) ) 1825 CV_ERROR( CV_StsUnmatchedSizes, "" ); 1826 real_transform = 1; 1827 } 1828 else if( inv && CV_MAT_CN(src->type) == 2 && CV_MAT_CN(dst->type) == 1 ) 1829 { 1830 if( (src->cols != 1 || dst->cols != 1 || 1831 (dst->rows/2+1 != src->rows && src->rows != dst->rows)) && 1832 (dst->cols/2+1 != src->cols || src->rows != dst->rows) ) 1833 CV_ERROR( CV_StsUnmatchedSizes, "" ); 1834 real_transform = 1; 1835 } 1836 else 1837 CV_ERROR( CV_StsUnmatchedFormats, 1838 "Incorrect or unsupported combination of input & output formats" ); 1839 1840 if( src->cols == 1 && nonzero_rows > 0 ) 1841 CV_ERROR( CV_StsNotImplemented, 1842 "This mode (using nonzero_rows with a single-column matrix) breaks the function logic, so it is prohibited.\n" 1843 "For fast convolution/correlation use 2-column matrix or single-row matrix instead" ); 1844 1845 // determine, which transform to do first - row-wise 1846 // (stage 0) or column-wise (stage 1) transform 1847 if( !(flags & CV_DXT_ROWS) && src->rows > 1 && 1848 ((src->cols == 1 && !CV_IS_MAT_CONT(src->type & dst->type)) || 1849 (src->cols > 1 && inv && real_transform)) ) 1850 stage = 1; 1851 1852 ipp_norm_flag = !(flags & CV_DXT_SCALE) ? 8 : (flags & CV_DXT_INVERSE) ? 2 : 1; 1853 1854 for(;;) 1855 { 1856 double scale = 1; 1857 uchar* wave = 0; 1858 int* itab = 0; 1859 uchar* ptr; 1860 int i, len, count, sz = 0; 1861 int use_buf = 0, odd_real = 0; 1862 CvDFTFunc dft_func; 1863 1864 if( stage == 0 ) // row-wise transform 1865 { 1866 len = !inv ? src->cols : dst->cols; 1867 count = src->rows; 1868 if( len == 1 && !(flags & CV_DXT_ROWS) ) 1869 { 1870 len = !inv ? src->rows : dst->rows; 1871 count = 1; 1872 } 1873 odd_real = real_transform && (len & 1); 1874 } 1875 else 1876 { 1877 len = dst->rows; 1878 count = !inv ? src0->cols : dst->cols; 1879 sz = 2*len*complex_elem_size; 1880 } 1881 1882 spec = 0; 1883 if( len*count >= 64 && icvDFTInitAlloc_R_32f_p != 0 ) // use IPP DFT if available 1884 { 1885 int ipp_sz = 0; 1886 1887 if( real_transform && stage == 0 ) 1888 { 1889 if( depth == CV_32F ) 1890 { 1891 if( spec_r ) 1892 IPPI_CALL( icvDFTFree_R_32f_p( spec_r )); 1893 IPPI_CALL( icvDFTInitAlloc_R_32f_p( 1894 &spec_r, len, ipp_norm_flag, cvAlgHintNone )); 1895 IPPI_CALL( icvDFTGetBufSize_R_32f_p( spec_r, &ipp_sz )); 1896 } 1897 else 1898 { 1899 if( spec_r ) 1900 IPPI_CALL( icvDFTFree_R_64f_p( spec_r )); 1901 IPPI_CALL( icvDFTInitAlloc_R_64f_p( 1902 &spec_r, len, ipp_norm_flag, cvAlgHintNone )); 1903 IPPI_CALL( icvDFTGetBufSize_R_64f_p( spec_r, &ipp_sz )); 1904 } 1905 spec = spec_r; 1906 } 1907 else 1908 { 1909 if( depth == CV_32F ) 1910 { 1911 if( spec_c ) 1912 IPPI_CALL( icvDFTFree_C_32fc_p( spec_c )); 1913 IPPI_CALL( icvDFTInitAlloc_C_32fc_p( 1914 &spec_c, len, ipp_norm_flag, cvAlgHintNone )); 1915 IPPI_CALL( icvDFTGetBufSize_C_32fc_p( spec_c, &ipp_sz )); 1916 } 1917 else 1918 { 1919 if( spec_c ) 1920 IPPI_CALL( icvDFTFree_C_64fc_p( spec_c )); 1921 IPPI_CALL( icvDFTInitAlloc_C_64fc_p( 1922 &spec_c, len, ipp_norm_flag, cvAlgHintNone )); 1923 IPPI_CALL( icvDFTGetBufSize_C_64fc_p( spec_c, &ipp_sz )); 1924 } 1925 spec = spec_c; 1926 } 1927 1928 sz += ipp_sz; 1929 } 1930 else 1931 { 1932 if( len != prev_len ) 1933 nf = icvDFTFactorize( len, factors ); 1934 1935 inplace_transform = factors[0] == factors[nf-1]; 1936 sz += len*(complex_elem_size + sizeof(int)); 1937 i = nf > 1 && (factors[0] & 1) == 0; 1938 if( (factors[i] & 1) != 0 && factors[i] > 5 ) 1939 sz += (factors[i]+1)*complex_elem_size; 1940 1941 if( (stage == 0 && ((src->data.ptr == dst->data.ptr && !inplace_transform) || odd_real)) || 1942 (stage == 1 && !inplace_transform) ) 1943 { 1944 use_buf = 1; 1945 sz += len*complex_elem_size; 1946 } 1947 } 1948 1949 if( sz > buf_size ) 1950 { 1951 prev_len = 0; // because we release the buffer, 1952 // force recalculation of 1953 // twiddle factors and permutation table 1954 if( !local_alloc && buffer ) 1955 cvFree( &buffer ); 1956 if( sz <= CV_MAX_LOCAL_DFT_SIZE ) 1957 { 1958 buf_size = sz = CV_MAX_LOCAL_DFT_SIZE; 1959 buffer = cvStackAlloc(sz + 32); 1960 local_alloc = 1; 1961 } 1962 else 1963 { 1964 CV_CALL( buffer = cvAlloc(sz + 32) ); 1965 buf_size = sz; 1966 local_alloc = 0; 1967 } 1968 } 1969 1970 ptr = (uchar*)buffer; 1971 if( !spec ) 1972 { 1973 wave = ptr; 1974 ptr += len*complex_elem_size; 1975 itab = (int*)ptr; 1976 ptr = (uchar*)cvAlignPtr( ptr + len*sizeof(int), 16 ); 1977 1978 if( len != prev_len || (!inplace_transform && inv && real_transform)) 1979 icvDFTInit( len, nf, factors, itab, complex_elem_size, 1980 wave, stage == 0 && inv && real_transform ); 1981 // otherwise reuse the tables calculated on the previous stage 1982 } 1983 1984 if( stage == 0 ) 1985 { 1986 uchar* tmp_buf = 0; 1987 int dptr_offset = 0; 1988 int dst_full_len = len*elem_size; 1989 int _flags = inv + (CV_MAT_CN(src->type) != CV_MAT_CN(dst->type) ? 1990 ICV_DFT_COMPLEX_INPUT_OR_OUTPUT : 0); 1991 if( use_buf ) 1992 { 1993 tmp_buf = ptr; 1994 ptr += len*complex_elem_size; 1995 if( odd_real && !inv && len > 1 && 1996 !(_flags & ICV_DFT_COMPLEX_INPUT_OR_OUTPUT)) 1997 dptr_offset = elem_size; 1998 } 1999 2000 if( !inv && (_flags & ICV_DFT_COMPLEX_INPUT_OR_OUTPUT) ) 2001 dst_full_len += (len & 1) ? elem_size : complex_elem_size; 2002 2003 dft_func = dft_tbl[(!real_transform ? 0 : !inv ? 1 : 2) + (depth == CV_64F)*3]; 2004 2005 if( count > 1 && !(flags & CV_DXT_ROWS) && (!inv || !real_transform) ) 2006 stage = 1; 2007 else if( flags & CV_DXT_SCALE ) 2008 scale = 1./(len * (flags & CV_DXT_ROWS ? 1 : count)); 2009 2010 if( nonzero_rows <= 0 || nonzero_rows > count ) 2011 nonzero_rows = count; 2012 2013 for( i = 0; i < nonzero_rows; i++ ) 2014 { 2015 uchar* sptr = src->data.ptr + i*src->step; 2016 uchar* dptr0 = dst->data.ptr + i*dst->step; 2017 uchar* dptr = dptr0; 2018 2019 if( tmp_buf ) 2020 dptr = tmp_buf; 2021 2022 dft_func( sptr, dptr, len, nf, factors, itab, wave, len, spec, ptr, _flags, scale ); 2023 if( dptr != dptr0 ) 2024 memcpy( dptr0, dptr + dptr_offset, dst_full_len ); 2025 } 2026 2027 for( ; i < count; i++ ) 2028 { 2029 uchar* dptr0 = dst->data.ptr + i*dst->step; 2030 memset( dptr0, 0, dst_full_len ); 2031 } 2032 2033 if( stage != 1 ) 2034 break; 2035 src = dst; 2036 } 2037 else 2038 { 2039 int a = 0, b = count; 2040 uchar *buf0, *buf1, *dbuf0, *dbuf1; 2041 uchar* sptr0 = src->data.ptr; 2042 uchar* dptr0 = dst->data.ptr; 2043 buf0 = ptr; 2044 ptr += len*complex_elem_size; 2045 buf1 = ptr; 2046 ptr += len*complex_elem_size; 2047 dbuf0 = buf0, dbuf1 = buf1; 2048 2049 if( use_buf ) 2050 { 2051 dbuf1 = ptr; 2052 dbuf0 = buf1; 2053 ptr += len*complex_elem_size; 2054 } 2055 2056 dft_func = dft_tbl[(depth == CV_64F)*3]; 2057 2058 if( real_transform && inv && src->cols > 1 ) 2059 stage = 0; 2060 else if( flags & CV_DXT_SCALE ) 2061 scale = 1./(len * count); 2062 2063 if( real_transform ) 2064 { 2065 int even; 2066 a = 1; 2067 even = (count & 1) == 0; 2068 b = (count+1)/2; 2069 if( !inv ) 2070 { 2071 memset( buf0, 0, len*complex_elem_size ); 2072 icvCopyColumn( sptr0, src->step, buf0, complex_elem_size, len, elem_size ); 2073 sptr0 += CV_MAT_CN(dst->type)*elem_size; 2074 if( even ) 2075 { 2076 memset( buf1, 0, len*complex_elem_size ); 2077 icvCopyColumn( sptr0 + (count-2)*elem_size, src->step, 2078 buf1, complex_elem_size, len, elem_size ); 2079 } 2080 } 2081 else if( CV_MAT_CN(src->type) == 1 ) 2082 { 2083 icvCopyColumn( sptr0, src->step, buf0 + elem_size, elem_size, len, elem_size ); 2084 icvExpandCCS( buf0 + elem_size, len, elem_size ); 2085 if( even ) 2086 { 2087 icvCopyColumn( sptr0 + (count-1)*elem_size, src->step, 2088 buf1 + elem_size, elem_size, len, elem_size ); 2089 icvExpandCCS( buf1 + elem_size, len, elem_size ); 2090 } 2091 sptr0 += elem_size; 2092 } 2093 else 2094 { 2095 icvCopyColumn( sptr0, src->step, buf0, complex_elem_size, len, complex_elem_size ); 2096 //memcpy( buf0 + elem_size, buf0, elem_size ); 2097 //icvExpandCCS( buf0 + elem_size, len, elem_size ); 2098 if( even ) 2099 { 2100 icvCopyColumn( sptr0 + b*complex_elem_size, src->step, 2101 buf1, complex_elem_size, len, complex_elem_size ); 2102 //memcpy( buf0 + elem_size, buf0, elem_size ); 2103 //icvExpandCCS( buf0 + elem_size, len, elem_size ); 2104 } 2105 sptr0 += complex_elem_size; 2106 } 2107 2108 if( even ) 2109 IPPI_CALL( dft_func( buf1, dbuf1, len, nf, factors, itab, 2110 wave, len, spec, ptr, inv, scale )); 2111 IPPI_CALL( dft_func( buf0, dbuf0, len, nf, factors, itab, 2112 wave, len, spec, ptr, inv, scale )); 2113 2114 if( CV_MAT_CN(dst->type) == 1 ) 2115 { 2116 if( !inv ) 2117 { 2118 // copy the half of output vector to the first/last column. 2119 // before doing that, defgragment the vector 2120 memcpy( dbuf0 + elem_size, dbuf0, elem_size ); 2121 icvCopyColumn( dbuf0 + elem_size, elem_size, dptr0, 2122 dst->step, len, elem_size ); 2123 if( even ) 2124 { 2125 memcpy( dbuf1 + elem_size, dbuf1, elem_size ); 2126 icvCopyColumn( dbuf1 + elem_size, elem_size, 2127 dptr0 + (count-1)*elem_size, 2128 dst->step, len, elem_size ); 2129 } 2130 dptr0 += elem_size; 2131 } 2132 else 2133 { 2134 // copy the real part of the complex vector to the first/last column 2135 icvCopyColumn( dbuf0, complex_elem_size, dptr0, dst->step, len, elem_size ); 2136 if( even ) 2137 icvCopyColumn( dbuf1, complex_elem_size, dptr0 + (count-1)*elem_size, 2138 dst->step, len, elem_size ); 2139 dptr0 += elem_size; 2140 } 2141 } 2142 else 2143 { 2144 assert( !inv ); 2145 icvCopyColumn( dbuf0, complex_elem_size, dptr0, 2146 dst->step, len, complex_elem_size ); 2147 if( even ) 2148 icvCopyColumn( dbuf1, complex_elem_size, 2149 dptr0 + b*complex_elem_size, 2150 dst->step, len, complex_elem_size ); 2151 dptr0 += complex_elem_size; 2152 } 2153 } 2154 2155 for( i = a; i < b; i += 2 ) 2156 { 2157 if( i+1 < b ) 2158 { 2159 icvCopyFrom2Columns( sptr0, src->step, buf0, buf1, len, complex_elem_size ); 2160 IPPI_CALL( dft_func( buf1, dbuf1, len, nf, factors, itab, 2161 wave, len, spec, ptr, inv, scale )); 2162 } 2163 else 2164 icvCopyColumn( sptr0, src->step, buf0, complex_elem_size, len, complex_elem_size ); 2165 2166 IPPI_CALL( dft_func( buf0, dbuf0, len, nf, factors, itab, 2167 wave, len, spec, ptr, inv, scale )); 2168 2169 if( i+1 < b ) 2170 icvCopyTo2Columns( dbuf0, dbuf1, dptr0, dst->step, len, complex_elem_size ); 2171 else 2172 icvCopyColumn( dbuf0, complex_elem_size, dptr0, dst->step, len, complex_elem_size ); 2173 sptr0 += 2*complex_elem_size; 2174 dptr0 += 2*complex_elem_size; 2175 } 2176 2177 if( stage != 0 ) 2178 break; 2179 src = dst; 2180 } 2181 } 2182 2183 __END__; 2184 2185 if( buffer && !local_alloc ) 2186 cvFree( &buffer ); 2187 2188 if( spec_c ) 2189 { 2190 if( depth == CV_32F ) 2191 icvDFTFree_C_32fc_p( spec_c ); 2192 else 2193 icvDFTFree_C_64fc_p( spec_c ); 2194 } 2195 2196 if( spec_r ) 2197 { 2198 if( depth == CV_32F ) 2199 icvDFTFree_R_32f_p( spec_r ); 2200 else 2201 icvDFTFree_R_64f_p( spec_r ); 2202 } 2203} 2204 2205 2206CV_IMPL void 2207cvMulSpectrums( const CvArr* srcAarr, const CvArr* srcBarr, 2208 CvArr* dstarr, int flags ) 2209{ 2210 CV_FUNCNAME( "cvMulSpectrums" ); 2211 2212 __BEGIN__; 2213 2214 CvMat stubA, *srcA = (CvMat*)srcAarr; 2215 CvMat stubB, *srcB = (CvMat*)srcBarr; 2216 CvMat dststub, *dst = (CvMat*)dstarr; 2217 int stepA, stepB, stepC; 2218 int type, cn, is_1d; 2219 int j, j0, j1, k, rows, cols, ncols; 2220 2221 if( !CV_IS_MAT(srcA)) 2222 CV_CALL( srcA = cvGetMat( srcA, &stubA, 0 )); 2223 2224 if( !CV_IS_MAT(srcB)) 2225 CV_CALL( srcB = cvGetMat( srcB, &stubB, 0 )); 2226 2227 if( !CV_IS_MAT(dst)) 2228 CV_CALL( dst = cvGetMat( dst, &dststub, 0 )); 2229 2230 if( !CV_ARE_TYPES_EQ( srcA, srcB ) || !CV_ARE_TYPES_EQ( srcA, dst )) 2231 CV_ERROR( CV_StsUnmatchedFormats, "" ); 2232 2233 if( !CV_ARE_SIZES_EQ( srcA, srcB ) || !CV_ARE_SIZES_EQ( srcA, dst )) 2234 CV_ERROR( CV_StsUnmatchedSizes, "" ); 2235 2236 type = CV_MAT_TYPE( dst->type ); 2237 cn = CV_MAT_CN(type); 2238 rows = srcA->rows; 2239 cols = srcA->cols; 2240 is_1d = (flags & CV_DXT_ROWS) || 2241 (rows == 1 || (cols == 1 && 2242 CV_IS_MAT_CONT( srcA->type & srcB->type & dst->type ))); 2243 2244 if( is_1d && !(flags & CV_DXT_ROWS) ) 2245 cols = cols + rows - 1, rows = 1; 2246 ncols = cols*cn; 2247 j0 = cn == 1; 2248 j1 = ncols - (cols % 2 == 0 && cn == 1); 2249 2250 if( CV_MAT_DEPTH(type) == CV_32F ) 2251 { 2252 float* dataA = srcA->data.fl; 2253 float* dataB = srcB->data.fl; 2254 float* dataC = dst->data.fl; 2255 2256 stepA = srcA->step/sizeof(dataA[0]); 2257 stepB = srcB->step/sizeof(dataB[0]); 2258 stepC = dst->step/sizeof(dataC[0]); 2259 2260 if( !is_1d && cn == 1 ) 2261 { 2262 for( k = 0; k < (cols % 2 ? 1 : 2); k++ ) 2263 { 2264 if( k == 1 ) 2265 dataA += cols - 1, dataB += cols - 1, dataC += cols - 1; 2266 dataC[0] = dataA[0]*dataB[0]; 2267 if( rows % 2 == 0 ) 2268 dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA]*dataB[(rows-1)*stepB]; 2269 if( !(flags & CV_DXT_MUL_CONJ) ) 2270 for( j = 1; j <= rows - 2; j += 2 ) 2271 { 2272 double re = (double)dataA[j*stepA]*dataB[j*stepB] - 2273 (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; 2274 double im = (double)dataA[j*stepA]*dataB[(j+1)*stepB] + 2275 (double)dataA[(j+1)*stepA]*dataB[j*stepB]; 2276 dataC[j*stepC] = (float)re; dataC[(j+1)*stepC] = (float)im; 2277 } 2278 else 2279 for( j = 1; j <= rows - 2; j += 2 ) 2280 { 2281 double re = (double)dataA[j*stepA]*dataB[j*stepB] + 2282 (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; 2283 double im = (double)dataA[(j+1)*stepA]*dataB[j*stepB] - 2284 (double)dataA[j*stepA]*dataB[(j+1)*stepB]; 2285 dataC[j*stepC] = (float)re; dataC[(j+1)*stepC] = (float)im; 2286 } 2287 if( k == 1 ) 2288 dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1; 2289 } 2290 } 2291 2292 for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC ) 2293 { 2294 if( is_1d && cn == 1 ) 2295 { 2296 dataC[0] = dataA[0]*dataB[0]; 2297 if( cols % 2 == 0 ) 2298 dataC[j1] = dataA[j1]*dataB[j1]; 2299 } 2300 2301 if( !(flags & CV_DXT_MUL_CONJ) ) 2302 for( j = j0; j < j1; j += 2 ) 2303 { 2304 double re = (double)dataA[j]*dataB[j] - (double)dataA[j+1]*dataB[j+1]; 2305 double im = (double)dataA[j+1]*dataB[j] + (double)dataA[j]*dataB[j+1]; 2306 dataC[j] = (float)re; dataC[j+1] = (float)im; 2307 } 2308 else 2309 for( j = j0; j < j1; j += 2 ) 2310 { 2311 double re = (double)dataA[j]*dataB[j] + (double)dataA[j+1]*dataB[j+1]; 2312 double im = (double)dataA[j+1]*dataB[j] - (double)dataA[j]*dataB[j+1]; 2313 dataC[j] = (float)re; dataC[j+1] = (float)im; 2314 } 2315 } 2316 } 2317 else if( CV_MAT_DEPTH(type) == CV_64F ) 2318 { 2319 double* dataA = srcA->data.db; 2320 double* dataB = srcB->data.db; 2321 double* dataC = dst->data.db; 2322 2323 stepA = srcA->step/sizeof(dataA[0]); 2324 stepB = srcB->step/sizeof(dataB[0]); 2325 stepC = dst->step/sizeof(dataC[0]); 2326 2327 if( !is_1d && cn == 1 ) 2328 { 2329 for( k = 0; k < (cols % 2 ? 1 : 2); k++ ) 2330 { 2331 if( k == 1 ) 2332 dataA += cols - 1, dataB += cols - 1, dataC += cols - 1; 2333 dataC[0] = dataA[0]*dataB[0]; 2334 if( rows % 2 == 0 ) 2335 dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA]*dataB[(rows-1)*stepB]; 2336 if( !(flags & CV_DXT_MUL_CONJ) ) 2337 for( j = 1; j <= rows - 2; j += 2 ) 2338 { 2339 double re = dataA[j*stepA]*dataB[j*stepB] - 2340 dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; 2341 double im = dataA[j*stepA]*dataB[(j+1)*stepB] + 2342 dataA[(j+1)*stepA]*dataB[j*stepB]; 2343 dataC[j*stepC] = re; dataC[(j+1)*stepC] = im; 2344 } 2345 else 2346 for( j = 1; j <= rows - 2; j += 2 ) 2347 { 2348 double re = dataA[j*stepA]*dataB[j*stepB] + 2349 dataA[(j+1)*stepA]*dataB[(j+1)*stepB]; 2350 double im = dataA[(j+1)*stepA]*dataB[j*stepB] - 2351 dataA[j*stepA]*dataB[(j+1)*stepB]; 2352 dataC[j*stepC] = re; dataC[(j+1)*stepC] = im; 2353 } 2354 if( k == 1 ) 2355 dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1; 2356 } 2357 } 2358 2359 for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC ) 2360 { 2361 if( is_1d && cn == 1 ) 2362 { 2363 dataC[0] = dataA[0]*dataB[0]; 2364 if( cols % 2 == 0 ) 2365 dataC[j1] = dataA[j1]*dataB[j1]; 2366 } 2367 2368 if( !(flags & CV_DXT_MUL_CONJ) ) 2369 for( j = j0; j < j1; j += 2 ) 2370 { 2371 double re = dataA[j]*dataB[j] - dataA[j+1]*dataB[j+1]; 2372 double im = dataA[j+1]*dataB[j] + dataA[j]*dataB[j+1]; 2373 dataC[j] = re; dataC[j+1] = im; 2374 } 2375 else 2376 for( j = j0; j < j1; j += 2 ) 2377 { 2378 double re = dataA[j]*dataB[j] + dataA[j+1]*dataB[j+1]; 2379 double im = dataA[j+1]*dataB[j] - dataA[j]*dataB[j+1]; 2380 dataC[j] = re; dataC[j+1] = im; 2381 } 2382 } 2383 } 2384 else 2385 { 2386 CV_ERROR( CV_StsUnsupportedFormat, "Only 32f and 64f types are supported" ); 2387 } 2388 2389 __END__; 2390} 2391 2392 2393/****************************************************************************************\ 2394 Discrete Cosine Transform 2395\****************************************************************************************/ 2396 2397/* DCT is calculated using DFT, as described here: 2398 http://www.ece.utexas.edu/~bevans/courses/ee381k/lectures/09_DCT/lecture9/: 2399*/ 2400#define ICV_DCT_FWD( flavor, datatype ) \ 2401static CvStatus CV_STDCALL \ 2402icvDCT_fwd_##flavor( const datatype* src, int src_step, datatype* dft_src,\ 2403 datatype* dft_dst, datatype* dst, int dst_step, \ 2404 int n, int nf, int* factors, const int* itab, \ 2405 const CvComplex##flavor* dft_wave, \ 2406 const CvComplex##flavor* dct_wave, \ 2407 const void* spec, CvComplex##flavor* buf ) \ 2408{ \ 2409 int j, n2 = n >> 1; \ 2410 \ 2411 src_step /= sizeof(src[0]); \ 2412 dst_step /= sizeof(dst[0]); \ 2413 datatype* dst1 = dst + (n-1)*dst_step; \ 2414 \ 2415 if( n == 1 ) \ 2416 { \ 2417 dst[0] = src[0]; \ 2418 return CV_OK; \ 2419 } \ 2420 \ 2421 for( j = 0; j < n2; j++, src += src_step*2 ) \ 2422 { \ 2423 dft_src[j] = src[0]; \ 2424 dft_src[n-j-1] = src[src_step]; \ 2425 } \ 2426 \ 2427 icvRealDFT_##flavor( dft_src, dft_dst, n, nf, factors, \ 2428 itab, dft_wave, n, spec, buf, 0, 1.0 ); \ 2429 src = dft_dst; \ 2430 \ 2431 dst[0] = (datatype)(src[0]*dct_wave->re*icv_sin_45); \ 2432 dst += dst_step; \ 2433 for( j = 1, dct_wave++; j < n2; j++, dct_wave++, \ 2434 dst += dst_step, dst1 -= dst_step ) \ 2435 { \ 2436 double t0 = dct_wave->re*src[j*2-1] - dct_wave->im*src[j*2]; \ 2437 double t1 = -dct_wave->im*src[j*2-1] - dct_wave->re*src[j*2]; \ 2438 dst[0] = (datatype)t0; \ 2439 dst1[0] = (datatype)t1; \ 2440 } \ 2441 \ 2442 dst[0] = (datatype)(src[n-1]*dct_wave->re); \ 2443 return CV_OK; \ 2444} 2445 2446 2447#define ICV_DCT_INV( flavor, datatype ) \ 2448static CvStatus CV_STDCALL \ 2449icvDCT_inv_##flavor( const datatype* src, int src_step, datatype* dft_src,\ 2450 datatype* dft_dst, datatype* dst, int dst_step, \ 2451 int n, int nf, int* factors, const int* itab, \ 2452 const CvComplex##flavor* dft_wave, \ 2453 const CvComplex##flavor* dct_wave, \ 2454 const void* spec, CvComplex##flavor* buf ) \ 2455{ \ 2456 int j, n2 = n >> 1; \ 2457 \ 2458 src_step /= sizeof(src[0]); \ 2459 dst_step /= sizeof(dst[0]); \ 2460 const datatype* src1 = src + (n-1)*src_step; \ 2461 \ 2462 if( n == 1 ) \ 2463 { \ 2464 dst[0] = src[0]; \ 2465 return CV_OK; \ 2466 } \ 2467 \ 2468 dft_src[0] = (datatype)(src[0]*2*dct_wave->re*icv_sin_45); \ 2469 src += src_step; \ 2470 for( j = 1, dct_wave++; j < n2; j++, dct_wave++, \ 2471 src += src_step, src1 -= src_step ) \ 2472 { \ 2473 double t0 = dct_wave->re*src[0] - dct_wave->im*src1[0]; \ 2474 double t1 = -dct_wave->im*src[0] - dct_wave->re*src1[0]; \ 2475 dft_src[j*2-1] = (datatype)t0; \ 2476 dft_src[j*2] = (datatype)t1; \ 2477 } \ 2478 \ 2479 dft_src[n-1] = (datatype)(src[0]*2*dct_wave->re); \ 2480 icvCCSIDFT_##flavor( dft_src, dft_dst, n, nf, factors, itab, \ 2481 dft_wave, n, spec, buf, CV_DXT_INVERSE, 1.0 ); \ 2482 \ 2483 for( j = 0; j < n2; j++, dst += dst_step*2 ) \ 2484 { \ 2485 dst[0] = dft_dst[j]; \ 2486 dst[dst_step] = dft_dst[n-j-1]; \ 2487 } \ 2488 return CV_OK; \ 2489} 2490 2491 2492ICV_DCT_FWD( 64f, double ) 2493ICV_DCT_INV( 64f, double ) 2494ICV_DCT_FWD( 32f, float ) 2495ICV_DCT_INV( 32f, float ) 2496 2497static void 2498icvDCTInit( int n, int elem_size, void* _wave, int inv ) 2499{ 2500 static const double icvDctScale[] = 2501 { 2502 0.707106781186547570, 0.500000000000000000, 0.353553390593273790, 2503 0.250000000000000000, 0.176776695296636890, 0.125000000000000000, 2504 0.088388347648318447, 0.062500000000000000, 0.044194173824159223, 2505 0.031250000000000000, 0.022097086912079612, 0.015625000000000000, 2506 0.011048543456039806, 0.007812500000000000, 0.005524271728019903, 2507 0.003906250000000000, 0.002762135864009952, 0.001953125000000000, 2508 0.001381067932004976, 0.000976562500000000, 0.000690533966002488, 2509 0.000488281250000000, 0.000345266983001244, 0.000244140625000000, 2510 0.000172633491500622, 0.000122070312500000, 0.000086316745750311, 2511 0.000061035156250000, 0.000043158372875155, 0.000030517578125000 2512 }; 2513 2514 int i; 2515 CvComplex64f w, w1; 2516 double t, scale; 2517 2518 if( n == 1 ) 2519 return; 2520 2521 assert( (n&1) == 0 ); 2522 2523 if( (n & (n - 1)) == 0 ) 2524 { 2525 int m = icvlog2(n); 2526 scale = (!inv ? 2 : 1)*icvDctScale[m]; 2527 w1.re = icvDxtTab[m+2][0]; 2528 w1.im = -icvDxtTab[m+2][1]; 2529 } 2530 else 2531 { 2532 t = 1./(2*n); 2533 scale = (!inv ? 2 : 1)*sqrt(t); 2534 w1.im = sin(-CV_PI*t); 2535 w1.re = sqrt(1. - w1.im*w1.im); 2536 } 2537 n >>= 1; 2538 2539 if( elem_size == sizeof(CvComplex64f) ) 2540 { 2541 CvComplex64f* wave = (CvComplex64f*)_wave; 2542 2543 w.re = scale; 2544 w.im = 0.; 2545 2546 for( i = 0; i <= n; i++ ) 2547 { 2548 wave[i] = w; 2549 t = w.re*w1.re - w.im*w1.im; 2550 w.im = w.re*w1.im + w.im*w1.re; 2551 w.re = t; 2552 } 2553 } 2554 else 2555 { 2556 CvComplex32f* wave = (CvComplex32f*)_wave; 2557 assert( elem_size == sizeof(CvComplex32f) ); 2558 2559 w.re = (float)scale; 2560 w.im = 0.f; 2561 2562 for( i = 0; i <= n; i++ ) 2563 { 2564 wave[i].re = (float)w.re; 2565 wave[i].im = (float)w.im; 2566 t = w.re*w1.re - w.im*w1.im; 2567 w.im = w.re*w1.im + w.im*w1.re; 2568 w.re = t; 2569 } 2570 } 2571} 2572 2573 2574typedef CvStatus (CV_STDCALL * CvDCTFunc)( 2575 const void* src, int src_step, void* dft_src, 2576 void* dft_dst, void* dst, int dst_step, int n, 2577 int nf, int* factors, const int* itab, const void* dft_wave, 2578 const void* dct_wave, const void* spec, void* buf ); 2579 2580CV_IMPL void 2581cvDCT( const CvArr* srcarr, CvArr* dstarr, int flags ) 2582{ 2583 static CvDCTFunc dct_tbl[4]; 2584 static int inittab = 0; 2585 2586 void* buffer = 0; 2587 int local_alloc = 1; 2588 int inv = (flags & CV_DXT_INVERSE) != 0, depth = -1; 2589 void *spec_dft = 0, *spec = 0; 2590 2591 CV_FUNCNAME( "cvDCT" ); 2592 2593 __BEGIN__; 2594 2595 double scale = 1.; 2596 int prev_len = 0, buf_size = 0, nf = 0, stage, end_stage; 2597 CvMat *src = (CvMat*)srcarr, *dst = (CvMat*)dstarr; 2598 uchar *src_dft_buf = 0, *dst_dft_buf = 0; 2599 uchar *dft_wave = 0, *dct_wave = 0; 2600 int* itab = 0; 2601 uchar* ptr = 0; 2602 CvMat srcstub, dststub; 2603 int complex_elem_size, elem_size; 2604 int factors[34], inplace_transform; 2605 int i, len, count; 2606 CvDCTFunc dct_func; 2607 2608 if( !inittab ) 2609 { 2610 dct_tbl[0] = (CvDCTFunc)icvDCT_fwd_32f; 2611 dct_tbl[1] = (CvDCTFunc)icvDCT_inv_32f; 2612 dct_tbl[2] = (CvDCTFunc)icvDCT_fwd_64f; 2613 dct_tbl[3] = (CvDCTFunc)icvDCT_inv_64f; 2614 inittab = 1; 2615 } 2616 2617 if( !CV_IS_MAT( src )) 2618 { 2619 int coi = 0; 2620 CV_CALL( src = cvGetMat( src, &srcstub, &coi )); 2621 2622 if( coi != 0 ) 2623 CV_ERROR( CV_BadCOI, "" ); 2624 } 2625 2626 if( !CV_IS_MAT( dst )) 2627 { 2628 int coi = 0; 2629 CV_CALL( dst = cvGetMat( dst, &dststub, &coi )); 2630 2631 if( coi != 0 ) 2632 CV_ERROR( CV_BadCOI, "" ); 2633 } 2634 2635 depth = CV_MAT_DEPTH(src->type); 2636 elem_size = CV_ELEM_SIZE1(depth); 2637 complex_elem_size = elem_size*2; 2638 2639 // check types and sizes 2640 if( !CV_ARE_TYPES_EQ(src, dst) ) 2641 CV_ERROR( CV_StsUnmatchedFormats, "" ); 2642 2643 if( depth < CV_32F || CV_MAT_CN(src->type) != 1 ) 2644 CV_ERROR( CV_StsUnsupportedFormat, 2645 "Only 32fC1 and 64fC1 formats are supported" ); 2646 2647 dct_func = dct_tbl[inv + (depth == CV_64F)*2]; 2648 2649 if( (flags & CV_DXT_ROWS) || src->rows == 1 || 2650 (src->cols == 1 && CV_IS_MAT_CONT(src->type & dst->type))) 2651 { 2652 stage = end_stage = 0; 2653 } 2654 else 2655 { 2656 stage = src->cols == 1; 2657 end_stage = 1; 2658 } 2659 2660 for( ; stage <= end_stage; stage++ ) 2661 { 2662 uchar *sptr = src->data.ptr, *dptr = dst->data.ptr; 2663 int sstep0, sstep1, dstep0, dstep1; 2664 2665 if( stage == 0 ) 2666 { 2667 len = src->cols; 2668 count = src->rows; 2669 if( len == 1 && !(flags & CV_DXT_ROWS) ) 2670 { 2671 len = src->rows; 2672 count = 1; 2673 } 2674 sstep0 = src->step; 2675 dstep0 = dst->step; 2676 sstep1 = dstep1 = elem_size; 2677 } 2678 else 2679 { 2680 len = dst->rows; 2681 count = dst->cols; 2682 sstep1 = src->step; 2683 dstep1 = dst->step; 2684 sstep0 = dstep0 = elem_size; 2685 } 2686 2687 if( len != prev_len ) 2688 { 2689 int sz; 2690 2691 if( len > 1 && (len & 1) ) 2692 CV_ERROR( CV_StsNotImplemented, "Odd-size DCT\'s are not implemented" ); 2693 2694 sz = len*elem_size; 2695 sz += (len/2 + 1)*complex_elem_size; 2696 2697 spec = 0; 2698 inplace_transform = 1; 2699 if( len*count >= 64 && icvDFTInitAlloc_R_32f_p ) 2700 { 2701 int ipp_sz = 0; 2702 if( depth == CV_32F ) 2703 { 2704 if( spec_dft ) 2705 IPPI_CALL( icvDFTFree_R_32f_p( spec_dft )); 2706 IPPI_CALL( icvDFTInitAlloc_R_32f_p( &spec_dft, len, 8, cvAlgHintNone )); 2707 IPPI_CALL( icvDFTGetBufSize_R_32f_p( spec_dft, &ipp_sz )); 2708 } 2709 else 2710 { 2711 if( spec_dft ) 2712 IPPI_CALL( icvDFTFree_R_64f_p( spec_dft )); 2713 IPPI_CALL( icvDFTInitAlloc_R_64f_p( &spec_dft, len, 8, cvAlgHintNone )); 2714 IPPI_CALL( icvDFTGetBufSize_R_64f_p( spec_dft, &ipp_sz )); 2715 } 2716 spec = spec_dft; 2717 sz += ipp_sz; 2718 } 2719 else 2720 { 2721 sz += len*(complex_elem_size + sizeof(int)) + complex_elem_size; 2722 2723 nf = icvDFTFactorize( len, factors ); 2724 inplace_transform = factors[0] == factors[nf-1]; 2725 2726 i = nf > 1 && (factors[0] & 1) == 0; 2727 if( (factors[i] & 1) != 0 && factors[i] > 5 ) 2728 sz += (factors[i]+1)*complex_elem_size; 2729 2730 if( !inplace_transform ) 2731 sz += len*elem_size; 2732 } 2733 2734 if( sz > buf_size ) 2735 { 2736 if( !local_alloc && buffer ) 2737 cvFree( &buffer ); 2738 if( sz <= CV_MAX_LOCAL_DFT_SIZE ) 2739 { 2740 buf_size = sz = CV_MAX_LOCAL_DFT_SIZE; 2741 buffer = cvStackAlloc(sz + 32); 2742 local_alloc = 1; 2743 } 2744 else 2745 { 2746 CV_CALL( buffer = cvAlloc(sz + 32) ); 2747 buf_size = sz; 2748 local_alloc = 0; 2749 } 2750 } 2751 2752 ptr = (uchar*)buffer; 2753 if( !spec ) 2754 { 2755 dft_wave = ptr; 2756 ptr += len*complex_elem_size; 2757 itab = (int*)ptr; 2758 ptr = (uchar*)cvAlignPtr( ptr + len*sizeof(int), 16 ); 2759 icvDFTInit( len, nf, factors, itab, complex_elem_size, dft_wave, inv ); 2760 } 2761 2762 dct_wave = ptr; 2763 ptr += (len/2 + 1)*complex_elem_size; 2764 src_dft_buf = dst_dft_buf = ptr; 2765 ptr += len*elem_size; 2766 if( !inplace_transform ) 2767 { 2768 dst_dft_buf = ptr; 2769 ptr += len*elem_size; 2770 } 2771 icvDCTInit( len, complex_elem_size, dct_wave, inv ); 2772 if( !inv ) 2773 scale += scale; 2774 prev_len = len; 2775 } 2776 // otherwise reuse the tables calculated on the previous stage 2777 2778 for( i = 0; i < count; i++ ) 2779 { 2780 dct_func( sptr + i*sstep0, sstep1, src_dft_buf, dst_dft_buf, 2781 dptr + i*dstep0, dstep1, len, nf, factors, 2782 itab, dft_wave, dct_wave, spec, ptr ); 2783 } 2784 src = dst; 2785 } 2786 2787 __END__; 2788 2789 if( spec_dft ) 2790 { 2791 if( depth == CV_32F ) 2792 icvDFTFree_R_32f_p( spec_dft ); 2793 else 2794 icvDFTFree_R_64f_p( spec_dft ); 2795 } 2796 2797 if( buffer && !local_alloc ) 2798 cvFree( &buffer ); 2799} 2800 2801 2802static const int icvOptimalDFTSize[] = { 28031, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 280450, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, 128, 135, 144, 150, 160, 2805162, 180, 192, 200, 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375, 2806384, 400, 405, 432, 450, 480, 486, 500, 512, 540, 576, 600, 625, 640, 648, 675, 720, 2807729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 1080, 1125, 1152, 1200, 28081215, 1250, 1280, 1296, 1350, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875, 28091920, 1944, 2000, 2025, 2048, 2160, 2187, 2250, 2304, 2400, 2430, 2500, 2560, 2592, 28102700, 2880, 2916, 3000, 3072, 3125, 3200, 3240, 3375, 3456, 3600, 3645, 3750, 3840, 28113888, 4000, 4050, 4096, 4320, 4374, 4500, 4608, 4800, 4860, 5000, 5120, 5184, 5400, 28125625, 5760, 5832, 6000, 6075, 6144, 6250, 6400, 6480, 6561, 6750, 6912, 7200, 7290, 28137500, 7680, 7776, 8000, 8100, 8192, 8640, 8748, 9000, 9216, 9375, 9600, 9720, 10000, 281410125, 10240, 10368, 10800, 10935, 11250, 11520, 11664, 12000, 12150, 12288, 12500, 281512800, 12960, 13122, 13500, 13824, 14400, 14580, 15000, 15360, 15552, 15625, 16000, 281616200, 16384, 16875, 17280, 17496, 18000, 18225, 18432, 18750, 19200, 19440, 19683, 281720000, 20250, 20480, 20736, 21600, 21870, 22500, 23040, 23328, 24000, 24300, 24576, 281825000, 25600, 25920, 26244, 27000, 27648, 28125, 28800, 29160, 30000, 30375, 30720, 281931104, 31250, 32000, 32400, 32768, 32805, 33750, 34560, 34992, 36000, 36450, 36864, 282037500, 38400, 38880, 39366, 40000, 40500, 40960, 41472, 43200, 43740, 45000, 46080, 282146656, 46875, 48000, 48600, 49152, 50000, 50625, 51200, 51840, 52488, 54000, 54675, 282255296, 56250, 57600, 58320, 59049, 60000, 60750, 61440, 62208, 62500, 64000, 64800, 282365536, 65610, 67500, 69120, 69984, 72000, 72900, 73728, 75000, 76800, 77760, 78125, 282478732, 80000, 81000, 81920, 82944, 84375, 86400, 87480, 90000, 91125, 92160, 93312, 282593750, 96000, 97200, 98304, 98415, 100000, 101250, 102400, 103680, 104976, 108000, 2826109350, 110592, 112500, 115200, 116640, 118098, 120000, 121500, 122880, 124416, 125000, 2827128000, 129600, 131072, 131220, 135000, 138240, 139968, 140625, 144000, 145800, 147456, 2828150000, 151875, 153600, 155520, 156250, 157464, 160000, 162000, 163840, 164025, 165888, 2829168750, 172800, 174960, 177147, 180000, 182250, 184320, 186624, 187500, 192000, 194400, 2830196608, 196830, 200000, 202500, 204800, 207360, 209952, 216000, 218700, 221184, 225000, 2831230400, 233280, 234375, 236196, 240000, 243000, 245760, 248832, 250000, 253125, 256000, 2832259200, 262144, 262440, 270000, 273375, 276480, 279936, 281250, 288000, 291600, 294912, 2833295245, 300000, 303750, 307200, 311040, 312500, 314928, 320000, 324000, 327680, 328050, 2834331776, 337500, 345600, 349920, 354294, 360000, 364500, 368640, 373248, 375000, 384000, 2835388800, 390625, 393216, 393660, 400000, 405000, 409600, 414720, 419904, 421875, 432000, 2836437400, 442368, 450000, 455625, 460800, 466560, 468750, 472392, 480000, 486000, 491520, 2837492075, 497664, 500000, 506250, 512000, 518400, 524288, 524880, 531441, 540000, 546750, 2838552960, 559872, 562500, 576000, 583200, 589824, 590490, 600000, 607500, 614400, 622080, 2839625000, 629856, 640000, 648000, 655360, 656100, 663552, 675000, 691200, 699840, 703125, 2840708588, 720000, 729000, 737280, 746496, 750000, 759375, 768000, 777600, 781250, 786432, 2841787320, 800000, 810000, 819200, 820125, 829440, 839808, 843750, 864000, 874800, 884736, 2842885735, 900000, 911250, 921600, 933120, 937500, 944784, 960000, 972000, 983040, 984150, 2843995328, 1000000, 1012500, 1024000, 1036800, 1048576, 1049760, 1062882, 1080000, 1093500, 28441105920, 1119744, 1125000, 1152000, 1166400, 1171875, 1179648, 1180980, 1200000, 28451215000, 1228800, 1244160, 1250000, 1259712, 1265625, 1280000, 1296000, 1310720, 28461312200, 1327104, 1350000, 1366875, 1382400, 1399680, 1406250, 1417176, 1440000, 28471458000, 1474560, 1476225, 1492992, 1500000, 1518750, 1536000, 1555200, 1562500, 28481572864, 1574640, 1594323, 1600000, 1620000, 1638400, 1640250, 1658880, 1679616, 28491687500, 1728000, 1749600, 1769472, 1771470, 1800000, 1822500, 1843200, 1866240, 28501875000, 1889568, 1920000, 1944000, 1953125, 1966080, 1968300, 1990656, 2000000, 28512025000, 2048000, 2073600, 2097152, 2099520, 2109375, 2125764, 2160000, 2187000, 28522211840, 2239488, 2250000, 2278125, 2304000, 2332800, 2343750, 2359296, 2361960, 28532400000, 2430000, 2457600, 2460375, 2488320, 2500000, 2519424, 2531250, 2560000, 28542592000, 2621440, 2624400, 2654208, 2657205, 2700000, 2733750, 2764800, 2799360, 28552812500, 2834352, 2880000, 2916000, 2949120, 2952450, 2985984, 3000000, 3037500, 28563072000, 3110400, 3125000, 3145728, 3149280, 3188646, 3200000, 3240000, 3276800, 28573280500, 3317760, 3359232, 3375000, 3456000, 3499200, 3515625, 3538944, 3542940, 28583600000, 3645000, 3686400, 3732480, 3750000, 3779136, 3796875, 3840000, 3888000, 28593906250, 3932160, 3936600, 3981312, 4000000, 4050000, 4096000, 4100625, 4147200, 28604194304, 4199040, 4218750, 4251528, 4320000, 4374000, 4423680, 4428675, 4478976, 28614500000, 4556250, 4608000, 4665600, 4687500, 4718592, 4723920, 4782969, 4800000, 28624860000, 4915200, 4920750, 4976640, 5000000, 5038848, 5062500, 5120000, 5184000, 28635242880, 5248800, 5308416, 5314410, 5400000, 5467500, 5529600, 5598720, 5625000, 28645668704, 5760000, 5832000, 5859375, 5898240, 5904900, 5971968, 6000000, 6075000, 28656144000, 6220800, 6250000, 6291456, 6298560, 6328125, 6377292, 6400000, 6480000, 28666553600, 6561000, 6635520, 6718464, 6750000, 6834375, 6912000, 6998400, 7031250, 28677077888, 7085880, 7200000, 7290000, 7372800, 7381125, 7464960, 7500000, 7558272, 28687593750, 7680000, 7776000, 7812500, 7864320, 7873200, 7962624, 7971615, 8000000, 28698100000, 8192000, 8201250, 8294400, 8388608, 8398080, 8437500, 8503056, 8640000, 28708748000, 8847360, 8857350, 8957952, 9000000, 9112500, 9216000, 9331200, 9375000, 28719437184, 9447840, 9565938, 9600000, 9720000, 9765625, 9830400, 9841500, 9953280, 287210000000, 10077696, 10125000, 10240000, 10368000, 10485760, 10497600, 10546875, 10616832, 287310628820, 10800000, 10935000, 11059200, 11197440, 11250000, 11337408, 11390625, 11520000, 287411664000, 11718750, 11796480, 11809800, 11943936, 12000000, 12150000, 12288000, 12301875, 287512441600, 12500000, 12582912, 12597120, 12656250, 12754584, 12800000, 12960000, 13107200, 287613122000, 13271040, 13286025, 13436928, 13500000, 13668750, 13824000, 13996800, 14062500, 287714155776, 14171760, 14400000, 14580000, 14745600, 14762250, 14929920, 15000000, 15116544, 287815187500, 15360000, 15552000, 15625000, 15728640, 15746400, 15925248, 15943230, 16000000, 287916200000, 16384000, 16402500, 16588800, 16777216, 16796160, 16875000, 17006112, 17280000, 288017496000, 17578125, 17694720, 17714700, 17915904, 18000000, 18225000, 18432000, 18662400, 288118750000, 18874368, 18895680, 18984375, 19131876, 19200000, 19440000, 19531250, 19660800, 288219683000, 19906560, 20000000, 20155392, 20250000, 20480000, 20503125, 20736000, 20971520, 288320995200, 21093750, 21233664, 21257640, 21600000, 21870000, 22118400, 22143375, 22394880, 288422500000, 22674816, 22781250, 23040000, 23328000, 23437500, 23592960, 23619600, 23887872, 288523914845, 24000000, 24300000, 24576000, 24603750, 24883200, 25000000, 25165824, 25194240, 288625312500, 25509168, 25600000, 25920000, 26214400, 26244000, 26542080, 26572050, 26873856, 288727000000, 27337500, 27648000, 27993600, 28125000, 28311552, 28343520, 28800000, 29160000, 288829296875, 29491200, 29524500, 29859840, 30000000, 30233088, 30375000, 30720000, 31104000, 288931250000, 31457280, 31492800, 31640625, 31850496, 31886460, 32000000, 32400000, 32768000, 289032805000, 33177600, 33554432, 33592320, 33750000, 34012224, 34171875, 34560000, 34992000, 289135156250, 35389440, 35429400, 35831808, 36000000, 36450000, 36864000, 36905625, 37324800, 289237500000, 37748736, 37791360, 37968750, 38263752, 38400000, 38880000, 39062500, 39321600, 289339366000, 39813120, 39858075, 40000000, 40310784, 40500000, 40960000, 41006250, 41472000, 289441943040, 41990400, 42187500, 42467328, 42515280, 43200000, 43740000, 44236800, 44286750, 289544789760, 45000000, 45349632, 45562500, 46080000, 46656000, 46875000, 47185920, 47239200, 289647775744, 47829690, 48000000, 48600000, 48828125, 49152000, 49207500, 49766400, 50000000, 289750331648, 50388480, 50625000, 51018336, 51200000, 51840000, 52428800, 52488000, 52734375, 289853084160, 53144100, 53747712, 54000000, 54675000, 55296000, 55987200, 56250000, 56623104, 289956687040, 56953125, 57600000, 58320000, 58593750, 58982400, 59049000, 59719680, 60000000, 290060466176, 60750000, 61440000, 61509375, 62208000, 62500000, 62914560, 62985600, 63281250, 290163700992, 63772920, 64000000, 64800000, 65536000, 65610000, 66355200, 66430125, 67108864, 290267184640, 67500000, 68024448, 68343750, 69120000, 69984000, 70312500, 70778880, 70858800, 290371663616, 72000000, 72900000, 73728000, 73811250, 74649600, 75000000, 75497472, 75582720, 290475937500, 76527504, 76800000, 77760000, 78125000, 78643200, 78732000, 79626240, 79716150, 290580000000, 80621568, 81000000, 81920000, 82012500, 82944000, 83886080, 83980800, 84375000, 290684934656, 85030560, 86400000, 87480000, 87890625, 88473600, 88573500, 89579520, 90000000, 290790699264, 91125000, 92160000, 93312000, 93750000, 94371840, 94478400, 94921875, 95551488, 290895659380, 96000000, 97200000, 97656250, 98304000, 98415000, 99532800, 100000000, 2909100663296, 100776960, 101250000, 102036672, 102400000, 102515625, 103680000, 104857600, 2910104976000, 105468750, 106168320, 106288200, 107495424, 108000000, 109350000, 110592000, 2911110716875, 111974400, 112500000, 113246208, 113374080, 113906250, 115200000, 116640000, 2912117187500, 117964800, 118098000, 119439360, 119574225, 120000000, 120932352, 121500000, 2913122880000, 123018750, 124416000, 125000000, 125829120, 125971200, 126562500, 127401984, 2914127545840, 128000000, 129600000, 131072000, 131220000, 132710400, 132860250, 134217728, 2915134369280, 135000000, 136048896, 136687500, 138240000, 139968000, 140625000, 141557760, 2916141717600, 143327232, 144000000, 145800000, 146484375, 147456000, 147622500, 149299200, 2917150000000, 150994944, 151165440, 151875000, 153055008, 153600000, 155520000, 156250000, 2918157286400, 157464000, 158203125, 159252480, 159432300, 160000000, 161243136, 162000000, 2919163840000, 164025000, 165888000, 167772160, 167961600, 168750000, 169869312, 170061120, 2920170859375, 172800000, 174960000, 175781250, 176947200, 177147000, 179159040, 180000000, 2921181398528, 182250000, 184320000, 184528125, 186624000, 187500000, 188743680, 188956800, 2922189843750, 191102976, 191318760, 192000000, 194400000, 195312500, 196608000, 196830000, 2923199065600, 199290375, 200000000, 201326592, 201553920, 202500000, 204073344, 204800000, 2924205031250, 207360000, 209715200, 209952000, 210937500, 212336640, 212576400, 214990848, 2925216000000, 218700000, 221184000, 221433750, 223948800, 225000000, 226492416, 226748160, 2926227812500, 230400000, 233280000, 234375000, 235929600, 236196000, 238878720, 239148450, 2927240000000, 241864704, 243000000, 244140625, 245760000, 246037500, 248832000, 250000000, 2928251658240, 251942400, 253125000, 254803968, 255091680, 256000000, 259200000, 262144000, 2929262440000, 263671875, 265420800, 265720500, 268435456, 268738560, 270000000, 272097792, 2930273375000, 276480000, 279936000, 281250000, 283115520, 283435200, 284765625, 286654464, 2931288000000, 291600000, 292968750, 294912000, 295245000, 298598400, 300000000, 301989888, 2932302330880, 303750000, 306110016, 307200000, 307546875, 311040000, 312500000, 314572800, 2933314928000, 316406250, 318504960, 318864600, 320000000, 322486272, 324000000, 327680000, 2934328050000, 331776000, 332150625, 335544320, 335923200, 337500000, 339738624, 340122240, 2935341718750, 345600000, 349920000, 351562500, 353894400, 354294000, 358318080, 360000000, 2936362797056, 364500000, 368640000, 369056250, 373248000, 375000000, 377487360, 377913600, 2937379687500, 382205952, 382637520, 384000000, 388800000, 390625000, 393216000, 393660000, 2938398131200, 398580750, 400000000, 402653184, 403107840, 405000000, 408146688, 409600000, 2939410062500, 414720000, 419430400, 419904000, 421875000, 424673280, 425152800, 429981696, 2940432000000, 437400000, 439453125, 442368000, 442867500, 447897600, 450000000, 452984832, 2941453496320, 455625000, 460800000, 466560000, 468750000, 471859200, 472392000, 474609375, 2942477757440, 478296900, 480000000, 483729408, 486000000, 488281250, 491520000, 492075000, 2943497664000, 500000000, 503316480, 503884800, 506250000, 509607936, 510183360, 512000000, 2944512578125, 518400000, 524288000, 524880000, 527343750, 530841600, 531441000, 536870912, 2945537477120, 540000000, 544195584, 546750000, 552960000, 553584375, 559872000, 562500000, 2946566231040, 566870400, 569531250, 573308928, 576000000, 583200000, 585937500, 589824000, 2947590490000, 597196800, 597871125, 600000000, 603979776, 604661760, 607500000, 612220032, 2948614400000, 615093750, 622080000, 625000000, 629145600, 629856000, 632812500, 637009920, 2949637729200, 640000000, 644972544, 648000000, 655360000, 656100000, 663552000, 664301250, 2950671088640, 671846400, 675000000, 679477248, 680244480, 683437500, 691200000, 699840000, 2951703125000, 707788800, 708588000, 716636160, 720000000, 725594112, 729000000, 732421875, 2952737280000, 738112500, 746496000, 750000000, 754974720, 755827200, 759375000, 764411904, 2953765275040, 768000000, 777600000, 781250000, 786432000, 787320000, 791015625, 796262400, 2954797161500, 800000000, 805306368, 806215680, 810000000, 816293376, 819200000, 820125000, 2955829440000, 838860800, 839808000, 843750000, 849346560, 850305600, 854296875, 859963392, 2956864000000, 874800000, 878906250, 884736000, 885735000, 895795200, 900000000, 905969664, 2957906992640, 911250000, 921600000, 922640625, 933120000, 937500000, 943718400, 944784000, 2958949218750, 955514880, 956593800, 960000000, 967458816, 972000000, 976562500, 983040000, 2959984150000, 995328000, 996451875, 1000000000, 1006632960, 1007769600, 1012500000, 29601019215872, 1020366720, 1024000000, 1025156250, 1036800000, 1048576000, 1049760000, 29611054687500, 1061683200, 1062882000, 1073741824, 1074954240, 1080000000, 1088391168, 29621093500000, 1105920000, 1107168750, 1119744000, 1125000000, 1132462080, 1133740800, 29631139062500, 1146617856, 1152000000, 1166400000, 1171875000, 1179648000, 1180980000, 29641194393600, 1195742250, 1200000000, 1207959552, 1209323520, 1215000000, 1220703125, 29651224440064, 1228800000, 1230187500, 1244160000, 1250000000, 1258291200, 1259712000, 29661265625000, 1274019840, 1275458400, 1280000000, 1289945088, 1296000000, 1310720000, 29671312200000, 1318359375, 1327104000, 1328602500, 1342177280, 1343692800, 1350000000, 29681358954496, 1360488960, 1366875000, 1382400000, 1399680000, 1406250000, 1415577600, 29691417176000, 1423828125, 1433272320, 1440000000, 1451188224, 1458000000, 1464843750, 29701474560000, 1476225000, 1492992000, 1500000000, 1509949440, 1511654400, 1518750000, 29711528823808, 1530550080, 1536000000, 1537734375, 1555200000, 1562500000, 1572864000, 29721574640000, 1582031250, 1592524800, 1594323000, 1600000000, 1610612736, 1612431360, 29731620000000, 1632586752, 1638400000, 1640250000, 1658880000, 1660753125, 1677721600, 29741679616000, 1687500000, 1698693120, 1700611200, 1708593750, 1719926784, 1728000000, 29751749600000, 1757812500, 1769472000, 1771470000, 1791590400, 1800000000, 1811939328, 29761813985280, 1822500000, 1843200000, 1845281250, 1866240000, 1875000000, 1887436800, 29771889568000, 1898437500, 1911029760, 1913187600, 1920000000, 1934917632, 1944000000, 29781953125000, 1966080000, 1968300000, 1990656000, 1992903750, 2000000000, 2013265920, 29792015539200, 2025000000, 2038431744, 2040733440, 2048000000, 2050312500, 2073600000, 29802097152000, 2099520000, 2109375000, 2123366400, 2125764000 2981}; 2982 2983 2984CV_IMPL int 2985cvGetOptimalDFTSize( int size0 ) 2986{ 2987 int a = 0, b = sizeof(icvOptimalDFTSize)/sizeof(icvOptimalDFTSize[0]) - 1; 2988 if( (unsigned)size0 >= (unsigned)icvOptimalDFTSize[b] ) 2989 return -1; 2990 2991 while( a < b ) 2992 { 2993 int c = (a + b) >> 1; 2994 if( size0 <= icvOptimalDFTSize[c] ) 2995 b = c; 2996 else 2997 a = c+1; 2998 } 2999 3000 return icvOptimalDFTSize[b]; 3001} 3002 3003/* End of file. */ 3004