1 /*M///////////////////////////////////////////////////////////////////////////////////////
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install,
7 // copy or use the software.
10 // Intel License Agreement
11 // For Open Source Computer Vision Library
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Third party copyrights are property of their respective owners.
16 // Redistribution and use in source and binary forms, with or without modification,
17 // are permitted provided that the following conditions are met:
19 // * Redistribution's of source code must retain the above copyright notice,
20 // this list of conditions and the following disclaimer.
22 // * Redistribution's in binary form must reproduce the above copyright notice,
23 // this list of conditions and the following disclaimer in the documentation
24 // and/or other materials provided with the distribution.
26 // * The name of Intel Corporation may not be used to endorse or promote products
27 // derived from this software without specific prior written permission.
29 // This software is provided by the copyright holders and contributors "as is" and
30 // any express or implied warranties, including, but not limited to, the implied
31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
32 // In no event shall the Intel Corporation or contributors be liable for any direct,
33 // indirect, incidental, special, exemplary, or consequential damages
34 // (including, but not limited to, procurement of substitute goods or services;
35 // loss of use, data, or profits; or business interruption) however caused
36 // and on any theory of liability, whether in contract, strict liability,
37 // or tort (including negligence or otherwise) arising in any way out of
38 // the use of this software, even if advised of the possibility of such damage.
44 // On Win64 (IA64) optimized versions of DFT and DCT fail the tests
45 #if defined WIN64 && !defined EM64T
46 #pragma optimize("", off)
49 icvDFTInitAlloc_C_32fc_t icvDFTInitAlloc_C_32fc_p = 0;
50 icvDFTFree_C_32fc_t icvDFTFree_C_32fc_p = 0;
51 icvDFTGetBufSize_C_32fc_t icvDFTGetBufSize_C_32fc_p = 0;
52 icvDFTFwd_CToC_32fc_t icvDFTFwd_CToC_32fc_p = 0;
53 icvDFTInv_CToC_32fc_t icvDFTInv_CToC_32fc_p = 0;
55 icvDFTInitAlloc_C_64fc_t icvDFTInitAlloc_C_64fc_p = 0;
56 icvDFTFree_C_64fc_t icvDFTFree_C_64fc_p = 0;
57 icvDFTGetBufSize_C_64fc_t icvDFTGetBufSize_C_64fc_p = 0;
58 icvDFTFwd_CToC_64fc_t icvDFTFwd_CToC_64fc_p = 0;
59 icvDFTInv_CToC_64fc_t icvDFTInv_CToC_64fc_p = 0;
61 icvDFTInitAlloc_R_32f_t icvDFTInitAlloc_R_32f_p = 0;
62 icvDFTFree_R_32f_t icvDFTFree_R_32f_p = 0;
63 icvDFTGetBufSize_R_32f_t icvDFTGetBufSize_R_32f_p = 0;
64 icvDFTFwd_RToPack_32f_t icvDFTFwd_RToPack_32f_p = 0;
65 icvDFTInv_PackToR_32f_t icvDFTInv_PackToR_32f_p = 0;
67 icvDFTInitAlloc_R_64f_t icvDFTInitAlloc_R_64f_p = 0;
68 icvDFTFree_R_64f_t icvDFTFree_R_64f_p = 0;
69 icvDFTGetBufSize_R_64f_t icvDFTGetBufSize_R_64f_p = 0;
70 icvDFTFwd_RToPack_64f_t icvDFTFwd_RToPack_64f_p = 0;
71 icvDFTInv_PackToR_64f_t icvDFTInv_PackToR_64f_p = 0;
73 /*icvDCTFwdInitAlloc_32f_t icvDCTFwdInitAlloc_32f_p = 0;
74 icvDCTFwdFree_32f_t icvDCTFwdFree_32f_p = 0;
75 icvDCTFwdGetBufSize_32f_t icvDCTFwdGetBufSize_32f_p = 0;
76 icvDCTFwd_32f_t icvDCTFwd_32f_p = 0;
78 icvDCTInvInitAlloc_32f_t icvDCTInvInitAlloc_32f_p = 0;
79 icvDCTInvFree_32f_t icvDCTInvFree_32f_p = 0;
80 icvDCTInvGetBufSize_32f_t icvDCTInvGetBufSize_32f_p = 0;
81 icvDCTInv_32f_t icvDCTInv_32f_p = 0;
83 icvDCTFwdInitAlloc_64f_t icvDCTFwdInitAlloc_64f_p = 0;
84 icvDCTFwdFree_64f_t icvDCTFwdFree_64f_p = 0;
85 icvDCTFwdGetBufSize_64f_t icvDCTFwdGetBufSize_64f_p = 0;
86 icvDCTFwd_64f_t icvDCTFwd_64f_p = 0;
88 icvDCTInvInitAlloc_64f_t icvDCTInvInitAlloc_64f_p = 0;
89 icvDCTInvFree_64f_t icvDCTInvFree_64f_p = 0;
90 icvDCTInvGetBufSize_64f_t icvDCTInvGetBufSize_64f_p = 0;
91 icvDCTInv_64f_t icvDCTInv_64f_p = 0;*/
93 /****************************************************************************************\
94 Discrete Fourier Transform
95 \****************************************************************************************/
97 #define CV_MAX_LOCAL_DFT_SIZE (1 << 15)
99 static const uchar log2tab[] = { 0, 0, 1, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0 };
100 static int icvlog2( int n )
103 int f = (n >= (1 << 16))*16;
106 f = (n >= (1 << 8))*8;
109 f = (n >= (1 << 4))*4;
111 return m + f + log2tab[n];
114 static unsigned char icvRevTable[] =
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
134 static const double icvDxtTab[][2] =
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 }
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)))
177 icvDFTFactorize( int n, int* factors )
187 f = (((n - 1)^n)+1) >> 1;
191 n = f == n ? 1 : n/f;
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 );
221 icvDFTInit( int n0, int nf, int* factors, int* itab, int elem_size, void* _wave, int inv_itab )
223 int digits[34], radix[34];
224 int n = factors[0], m = 0;
237 for( i = 1; i < n0-1; i++ )
247 if( elem_size == sizeof(CvComplex64f) )
248 ((CvComplex64f*)_wave)[0] = CvComplex64f(1.,0.);
250 ((CvComplex32f*)_wave)[0] = CvComplex32f(1.f,0.f);
258 // radix[] is initialized from index 'nf' down to zero
262 for( i = 0; i < nf; i++ )
265 radix[nf-i-1] = radix[nf-i]*factors[nf-i-1];
268 if( inv_itab && factors[0] != factors[nf-1] )
273 int a = radix[1], na2 = n*a>>1, na4 = na2 >> 1;
284 for( i = 0; i <= n - 4; i += 4 )
286 j = (icvRevTable[i>>2]>>shift)*a;
290 itab[i+3] = j + na2 + na4;
296 for( i = 0; i < n; i += 4 )
299 j = icvBitRev(i4,shift)*a;
303 itab[i+3] = j + na2 + na4;
309 assert (nf >= 2); // because we read radix[2] which is uninitialized otherwise
310 for( i = n, j = radix[2]; i < n0; )
312 for( k = 0; k < n; k++ )
313 itab[i+k] = itab[k] + j;
317 for( k = 1; ++digits[k] >= factors[k]; k++ )
320 j += radix[k+2] - radix[k];
326 for( i = 0, j = 0;; )
332 for( k = 0; ++digits[k] >= factors[k]; k++ )
335 j += radix[k+2] - radix[k];
343 for( i = n0 & 1; i < n0; i += 2 )
353 if( (n0 & (n0-1)) == 0 )
355 w.re = w1.re = icvDxtTab[m][0];
356 w.im = w1.im = -icvDxtTab[m][1];
361 w.im = w1.im = sin(t);
362 w.re = w1.re = sqrt(1. - w1.im*w1.im);
366 if( elem_size == sizeof(CvComplex64f) )
368 CvComplex64f* wave = (CvComplex64f*)_wave;
379 for( i = 1; i < n; i++ )
382 wave[n0-i].re = w.re;
383 wave[n0-i].im = -w.im;
385 t = w.re*w1.re - w.im*w1.im;
386 w.im = w.re*w1.im + w.im*w1.re;
392 CvComplex32f* wave = (CvComplex32f*)_wave;
393 assert( elem_size == sizeof(CvComplex32f) );
404 for( i = 1; i < n; i++ )
406 wave[i].re = (float)w.re;
407 wave[i].im = (float)w.im;
408 wave[n0-i].re = (float)w.re;
409 wave[n0-i].im = (float)-w.im;
411 t = w.re*w1.re - w.im*w1.im;
412 w.im = w.re*w1.im + w.im*w1.re;
419 static const double icv_sin_120 = 0.86602540378443864676372317075294;
420 static const double icv_sin_45 = 0.70710678118654752440084436210485;
421 static const double icv_fft5_2 = 0.559016994374947424102293417182819;
422 static const double icv_fft5_3 = -0.951056516295153572116439333379382;
423 static const double icv_fft5_4 = -1.538841768587626701285145288018455;
424 static const double icv_fft5_5 = 0.363271264002680442947733378740309;
426 #define ICV_DFT_NO_PERMUTE 2
427 #define ICV_DFT_COMPLEX_INPUT_OR_OUTPUT 4
429 // mixed-radix complex discrete Fourier transform: double-precision version
430 static CvStatus CV_STDCALL
431 icvDFT_64fc( const CvComplex64f* src, CvComplex64f* dst, int n,
432 int nf, int* factors, const int* itab,
433 const CvComplex64f* wave, int tab_size,
434 const void* spec, CvComplex64f* buf,
435 int flags, double scale )
437 int n0 = n, f_idx, nx;
438 int inv = flags & CV_DXT_INVERSE;
439 int dw0 = tab_size, dw;
446 assert( icvDFTFwd_CToC_64fc_p != 0 && icvDFTInv_CToC_64fc_p != 0 );
448 icvDFTFwd_CToC_64fc_p( src, dst, spec, buf ):
449 icvDFTInv_CToC_64fc_p( src, dst, spec, buf );
452 tab_step = tab_size == n ? 1 : tab_size == n*2 ? 2 : tab_size/n;
457 assert( (flags & ICV_DFT_NO_PERMUTE) == 0 );
460 for( i = 0; i <= n - 2; i += 2, itab += 2*tab_step )
462 int k0 = itab[0], k1 = itab[tab_step];
463 assert( (unsigned)k0 < (unsigned)n && (unsigned)k1 < (unsigned)n );
464 dst[i] = src[k0]; dst[i+1] = src[k1];
472 for( i = 0; i <= n - 2; i += 2, itab += 2*tab_step )
474 int k0 = itab[0], k1 = itab[tab_step];
475 assert( (unsigned)k0 < (unsigned)n && (unsigned)k1 < (unsigned)n );
476 t.re = src[k0].re; t.im = -src[k0].im;
478 t.re = src[k1].re; t.im = -src[k1].im;
484 t.re = src[n-1].re; t.im = -src[n-1].im;
491 if( (flags & ICV_DFT_NO_PERMUTE) == 0 )
493 if( factors[0] != factors[nf-1] )
494 return CV_INPLACE_NOT_SUPPORTED_ERR;
500 CvComplex64f* dsth = dst + n2;
502 for( i = 0; i < n2; i += 2, itab += tab_step*2 )
505 assert( (unsigned)j < (unsigned)n2 );
507 CV_SWAP(dst[i+1], dsth[j], t);
510 CV_SWAP(dst[i], dst[j], t);
511 CV_SWAP(dsth[i+1], dsth[j+1], t);
519 for( i = 0; i < n; i++, itab += tab_step )
522 assert( (unsigned)j < (unsigned)n );
524 CV_SWAP(dst[i], dst[j], t);
531 for( i = 0; i <= n - 2; i += 2 )
533 double t0 = -dst[i].im;
534 double t1 = -dst[i+1].im;
535 dst[i].im = t0; dst[i+1].im = t1;
539 dst[n-1].im = -dst[n-1].im;
544 // 1. power-2 transforms
545 if( (factors[0] & 1) == 0 )
548 for( ; n*4 <= factors[0]; )
554 for( i = 0; i < n0; i += n )
558 double r0, i0, r1, i1, r2, i2, r3, i3, r4, i4;
563 r2 = v0[0].re; i2 = v0[0].im;
564 r1 = v0[nx].re; i1 = v0[nx].im;
566 r0 = r1 + r2; i0 = i1 + i2;
569 i3 = v1[nx].re; r3 = v1[nx].im;
570 i4 = v1[0].re; r4 = v1[0].im;
572 r1 = i4 + i3; i1 = r4 + r3;
573 r3 = r4 - r3; i3 = i3 - i4;
575 v0[0].re = r0 + r1; v0[0].im = i0 + i1;
576 v1[0].re = r0 - r1; v1[0].im = i0 - i1;
577 v0[nx].re = r2 + r3; v0[nx].im = i2 + i3;
578 v1[nx].re = r2 - r3; v1[nx].im = i2 - i3;
580 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 )
585 r2 = v0[nx].re*wave[dw*2].re - v0[nx].im*wave[dw*2].im;
586 i2 = v0[nx].re*wave[dw*2].im + v0[nx].im*wave[dw*2].re;
587 r0 = v1[0].re*wave[dw].im + v1[0].im*wave[dw].re;
588 i0 = v1[0].re*wave[dw].re - v1[0].im*wave[dw].im;
589 r3 = v1[nx].re*wave[dw*3].im + v1[nx].im*wave[dw*3].re;
590 i3 = v1[nx].re*wave[dw*3].re - v1[nx].im*wave[dw*3].im;
592 r1 = i0 + i3; i1 = r0 + r3;
593 r3 = r0 - r3; i3 = i3 - i0;
594 r4 = v0[0].re; i4 = v0[0].im;
596 r0 = r4 + r2; i0 = i4 + i2;
597 r2 = r4 - r2; i2 = i4 - i2;
599 v0[0].re = r0 + r1; v0[0].im = i0 + i1;
600 v1[0].re = r0 - r1; v1[0].im = i0 - i1;
601 v0[nx].re = r2 + r3; v0[nx].im = i2 + i3;
602 v1[nx].re = r2 - r3; v1[nx].im = i2 - i3;
607 for( ; n < factors[0]; )
609 // do the remaining radix-2 transform
614 for( i = 0; i < n0; i += n )
616 CvComplex64f* v = dst + i;
617 double r0 = v[0].re + v[nx].re;
618 double i0 = v[0].im + v[nx].im;
619 double r1 = v[0].re - v[nx].re;
620 double i1 = v[0].im - v[nx].im;
621 v[0].re = r0; v[0].im = i0;
622 v[nx].re = r1; v[nx].im = i1;
624 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 )
627 r1 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im;
628 i1 = v[nx].im*wave[dw].re + v[nx].re*wave[dw].im;
629 r0 = v[0].re; i0 = v[0].im;
631 v[0].re = r0 + r1; v[0].im = i0 + i1;
632 v[nx].re = r0 - r1; v[nx].im = i0 - i1;
638 // 2. all the other transforms
639 for( f_idx = (factors[0]&1) ? 0 : 1; f_idx < nf; f_idx++ )
641 int factor = factors[f_idx];
649 for( i = 0; i < n0; i += n )
651 CvComplex64f* v = dst + i;
653 double r1 = v[nx].re + v[nx*2].re;
654 double i1 = v[nx].im + v[nx*2].im;
657 double r2 = icv_sin_120*(v[nx].im - v[nx*2].im);
658 double i2 = icv_sin_120*(v[nx*2].re - v[nx].re);
659 v[0].re = r0 + r1; v[0].im = i0 + i1;
660 r0 -= 0.5*r1; i0 -= 0.5*i1;
661 v[nx].re = r0 + r2; v[nx].im = i0 + i2;
662 v[nx*2].re = r0 - r2; v[nx*2].im = i0 - i2;
664 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 )
667 r0 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im;
668 i0 = v[nx].re*wave[dw].im + v[nx].im*wave[dw].re;
669 i2 = v[nx*2].re*wave[dw*2].re - v[nx*2].im*wave[dw*2].im;
670 r2 = v[nx*2].re*wave[dw*2].im + v[nx*2].im*wave[dw*2].re;
671 r1 = r0 + i2; i1 = i0 + r2;
673 r2 = icv_sin_120*(i0 - r2); i2 = icv_sin_120*(i2 - r0);
674 r0 = v[0].re; i0 = v[0].im;
675 v[0].re = r0 + r1; v[0].im = i0 + i1;
676 r0 -= 0.5*r1; i0 -= 0.5*i1;
677 v[nx].re = r0 + r2; v[nx].im = i0 + i2;
678 v[nx*2].re = r0 - r2; v[nx*2].im = i0 - i2;
682 else if( factor == 5 )
685 for( i = 0; i < n0; i += n )
687 for( j = 0, dw = 0; j < nx; j++, dw += dw0 )
689 CvComplex64f* v0 = dst + i + j;
690 CvComplex64f* v1 = v0 + nx*2;
691 CvComplex64f* v2 = v1 + nx*2;
693 double r0, i0, r1, i1, r2, i2, r3, i3, r4, i4, r5, i5;
695 r3 = v0[nx].re*wave[dw].re - v0[nx].im*wave[dw].im;
696 i3 = v0[nx].re*wave[dw].im + v0[nx].im*wave[dw].re;
697 r2 = v2[0].re*wave[dw*4].re - v2[0].im*wave[dw*4].im;
698 i2 = v2[0].re*wave[dw*4].im + v2[0].im*wave[dw*4].re;
700 r1 = r3 + r2; i1 = i3 + i2;
703 r4 = v1[nx].re*wave[dw*3].re - v1[nx].im*wave[dw*3].im;
704 i4 = v1[nx].re*wave[dw*3].im + v1[nx].im*wave[dw*3].re;
705 r0 = v1[0].re*wave[dw*2].re - v1[0].im*wave[dw*2].im;
706 i0 = v1[0].re*wave[dw*2].im + v1[0].im*wave[dw*2].re;
708 r2 = r4 + r0; i2 = i4 + i0;
711 r0 = v0[0].re; i0 = v0[0].im;
712 r5 = r1 + r2; i5 = i1 + i2;
714 v0[0].re = r0 + r5; v0[0].im = i0 + i5;
716 r0 -= 0.25*r5; i0 -= 0.25*i5;
717 r1 = icv_fft5_2*(r1 - r2); i1 = icv_fft5_2*(i1 - i2);
718 r2 = -icv_fft5_3*(i3 + i4); i2 = icv_fft5_3*(r3 + r4);
720 i3 *= -icv_fft5_5; r3 *= icv_fft5_5;
721 i4 *= -icv_fft5_4; r4 *= icv_fft5_4;
723 r5 = r2 + i3; i5 = i2 + r3;
726 r3 = r0 + r1; i3 = i0 + i1;
729 v0[nx].re = r3 + r2; v0[nx].im = i3 + i2;
730 v2[0].re = r3 - r2; v2[0].im = i3 - i2;
732 v1[0].re = r0 + r5; v1[0].im = i0 + i5;
733 v1[nx].re = r0 - r5; v1[nx].im = i0 - i5;
739 // radix-"factor" - an odd number
740 int p, q, factor2 = (factor - 1)/2;
741 int d, dd, dw_f = tab_size/factor;
742 CvComplex64f* a = buf;
743 CvComplex64f* b = buf + factor2;
745 for( i = 0; i < n0; i += n )
747 for( j = 0, dw = 0; j < nx; j++, dw += dw0 )
749 CvComplex64f* v = dst + i + j;
750 CvComplex64f v_0 = v[0];
751 CvComplex64f vn_0 = v_0;
755 for( p = 1, k = nx; p <= factor2; p++, k += nx )
757 double r0 = v[k].re + v[n-k].re;
758 double i0 = v[k].im - v[n-k].im;
759 double r1 = v[k].re - v[n-k].re;
760 double i1 = v[k].im + v[n-k].im;
762 vn_0.re += r0; vn_0.im += i1;
763 a[p-1].re = r0; a[p-1].im = i0;
764 b[p-1].re = r1; b[p-1].im = i1;
769 const CvComplex64f* wave_ = wave + dw*factor;
772 for( p = 1, k = nx; p <= factor2; p++, k += nx, d += dw )
774 double r2 = v[k].re*wave[d].re - v[k].im*wave[d].im;
775 double i2 = v[k].re*wave[d].im + v[k].im*wave[d].re;
777 double r1 = v[n-k].re*wave_[-d].re - v[n-k].im*wave_[-d].im;
778 double i1 = v[n-k].re*wave_[-d].im + v[n-k].im*wave_[-d].re;
785 vn_0.re += r0; vn_0.im += i1;
786 a[p-1].re = r0; a[p-1].im = i0;
787 b[p-1].re = r1; b[p-1].im = i1;
793 for( p = 1, k = nx; p <= factor2; p++, k += nx )
795 CvComplex64f s0 = v_0, s1 = v_0;
798 for( q = 0; q < factor2; q++ )
800 double r0 = wave[d].re * a[q].re;
801 double i0 = wave[d].im * a[q].im;
802 double r1 = wave[d].re * b[q].im;
803 double i1 = wave[d].im * b[q].re;
805 s1.re += r0 + i0; s0.re += r0 - i0;
806 s1.im += r1 - i1; s0.im += r1 + i1;
809 d -= -(d >= tab_size) & tab_size;
820 if( fabs(scale - 1.) > DBL_EPSILON )
822 double re_scale = scale, im_scale = scale;
824 im_scale = -im_scale;
826 for( i = 0; i < n0; i++ )
828 double t0 = dst[i].re*re_scale;
829 double t1 = dst[i].im*im_scale;
836 for( i = 0; i <= n0 - 2; i += 2 )
838 double t0 = -dst[i].im;
839 double t1 = -dst[i+1].im;
845 dst[n0-1].im = -dst[n0-1].im;
852 // mixed-radix complex discrete Fourier transform: single-precision version
853 static CvStatus CV_STDCALL
854 icvDFT_32fc( const CvComplex32f* src, CvComplex32f* dst, int n,
855 int nf, int* factors, const int* itab,
856 const CvComplex32f* wave, int tab_size,
857 const void* spec, CvComplex32f* buf,
858 int flags, double scale )
860 int n0 = n, f_idx, nx;
861 int inv = flags & CV_DXT_INVERSE;
862 int dw0 = tab_size, dw;
865 int tab_step = tab_size == n ? 1 : tab_size == n*2 ? 2 : tab_size/n;
869 assert( icvDFTFwd_CToC_32fc_p != 0 && icvDFTInv_CToC_32fc_p != 0 );
871 icvDFTFwd_CToC_32fc_p( src, dst, spec, buf ):
872 icvDFTInv_CToC_32fc_p( src, dst, spec, buf );
878 assert( (flags & ICV_DFT_NO_PERMUTE) == 0 );
881 for( i = 0; i <= n - 2; i += 2, itab += 2*tab_step )
883 int k0 = itab[0], k1 = itab[tab_step];
884 assert( (unsigned)k0 < (unsigned)n && (unsigned)k1 < (unsigned)n );
885 dst[i] = src[k0]; dst[i+1] = src[k1];
893 for( i = 0; i <= n - 2; i += 2, itab += 2*tab_step )
895 int k0 = itab[0], k1 = itab[tab_step];
896 assert( (unsigned)k0 < (unsigned)n && (unsigned)k1 < (unsigned)n );
897 t.re = src[k0].re; t.im = -src[k0].im;
899 t.re = src[k1].re; t.im = -src[k1].im;
905 t.re = src[n-1].re; t.im = -src[n-1].im;
912 if( (flags & ICV_DFT_NO_PERMUTE) == 0 )
914 if( factors[0] != factors[nf-1] )
915 return CV_INPLACE_NOT_SUPPORTED_ERR;
921 CvComplex32f* dsth = dst + n2;
923 for( i = 0; i < n2; i += 2, itab += tab_step*2 )
926 assert( (unsigned)j < (unsigned)n2 );
928 CV_SWAP(dst[i+1], dsth[j], t);
931 CV_SWAP(dst[i], dst[j], t);
932 CV_SWAP(dsth[i+1], dsth[j+1], t);
945 assert( (unsigned)j < (unsigned)n );
947 CV_SWAP(dst[i], dst[j], t);
955 // conjugate the vector - i.e. invert sign of the imaginary part
956 int* idst = (int*)dst;
957 for( i = 0; i <= n - 2; i += 2 )
959 int t0 = idst[i*2+1] ^ 0x80000000;
960 int t1 = idst[i*2+3] ^ 0x80000000;
961 idst[i*2+1] = t0; idst[i*2+3] = t1;
965 idst[2*i+1] ^= 0x80000000;
970 // 1. power-2 transforms
971 if( (factors[0] & 1) == 0 )
974 for( ; n*4 <= factors[0]; )
980 for( i = 0; i < n0; i += n )
984 double r0, i0, r1, i1, r2, i2, r3, i3, r4, i4;
989 r2 = v0[0].re; i2 = v0[0].im;
990 r1 = v0[nx].re; i1 = v0[nx].im;
992 r0 = r1 + r2; i0 = i1 + i2;
995 i3 = v1[nx].re; r3 = v1[nx].im;
996 i4 = v1[0].re; r4 = v1[0].im;
998 r1 = i4 + i3; i1 = r4 + r3;
999 r3 = r4 - r3; i3 = i3 - i4;
1001 v0[0].re = (float)(r0 + r1); v0[0].im = (float)(i0 + i1);
1002 v1[0].re = (float)(r0 - r1); v1[0].im = (float)(i0 - i1);
1003 v0[nx].re = (float)(r2 + r3); v0[nx].im = (float)(i2 + i3);
1004 v1[nx].re = (float)(r2 - r3); v1[nx].im = (float)(i2 - i3);
1006 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 )
1011 r2 = v0[nx].re*wave[dw*2].re - v0[nx].im*wave[dw*2].im;
1012 i2 = v0[nx].re*wave[dw*2].im + v0[nx].im*wave[dw*2].re;
1013 r0 = v1[0].re*wave[dw].im + v1[0].im*wave[dw].re;
1014 i0 = v1[0].re*wave[dw].re - v1[0].im*wave[dw].im;
1015 r3 = v1[nx].re*wave[dw*3].im + v1[nx].im*wave[dw*3].re;
1016 i3 = v1[nx].re*wave[dw*3].re - v1[nx].im*wave[dw*3].im;
1018 r1 = i0 + i3; i1 = r0 + r3;
1019 r3 = r0 - r3; i3 = i3 - i0;
1020 r4 = v0[0].re; i4 = v0[0].im;
1022 r0 = r4 + r2; i0 = i4 + i2;
1023 r2 = r4 - r2; i2 = i4 - i2;
1025 v0[0].re = (float)(r0 + r1); v0[0].im = (float)(i0 + i1);
1026 v1[0].re = (float)(r0 - r1); v1[0].im = (float)(i0 - i1);
1027 v0[nx].re = (float)(r2 + r3); v0[nx].im = (float)(i2 + i3);
1028 v1[nx].re = (float)(r2 - r3); v1[nx].im = (float)(i2 - i3);
1033 for( ; n < factors[0]; )
1035 // do the remaining radix-2 transform
1040 for( i = 0; i < n0; i += n )
1042 CvComplex32f* v = dst + i;
1043 double r0 = v[0].re + v[nx].re;
1044 double i0 = v[0].im + v[nx].im;
1045 double r1 = v[0].re - v[nx].re;
1046 double i1 = v[0].im - v[nx].im;
1047 v[0].re = (float)r0; v[0].im = (float)i0;
1048 v[nx].re = (float)r1; v[nx].im = (float)i1;
1050 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 )
1053 r1 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im;
1054 i1 = v[nx].im*wave[dw].re + v[nx].re*wave[dw].im;
1055 r0 = v[0].re; i0 = v[0].im;
1057 v[0].re = (float)(r0 + r1); v[0].im = (float)(i0 + i1);
1058 v[nx].re = (float)(r0 - r1); v[nx].im = (float)(i0 - i1);
1064 // 2. all the other transforms
1065 for( f_idx = (factors[0]&1) ? 0 : 1; f_idx < nf; f_idx++ )
1067 int factor = factors[f_idx];
1075 for( i = 0; i < n0; i += n )
1077 CvComplex32f* v = dst + i;
1079 double r1 = v[nx].re + v[nx*2].re;
1080 double i1 = v[nx].im + v[nx*2].im;
1081 double r0 = v[0].re;
1082 double i0 = v[0].im;
1083 double r2 = icv_sin_120*(v[nx].im - v[nx*2].im);
1084 double i2 = icv_sin_120*(v[nx*2].re - v[nx].re);
1085 v[0].re = (float)(r0 + r1); v[0].im = (float)(i0 + i1);
1086 r0 -= 0.5*r1; i0 -= 0.5*i1;
1087 v[nx].re = (float)(r0 + r2); v[nx].im = (float)(i0 + i2);
1088 v[nx*2].re = (float)(r0 - r2); v[nx*2].im = (float)(i0 - i2);
1090 for( j = 1, dw = dw0; j < nx; j++, dw += dw0 )
1093 r0 = v[nx].re*wave[dw].re - v[nx].im*wave[dw].im;
1094 i0 = v[nx].re*wave[dw].im + v[nx].im*wave[dw].re;
1095 i2 = v[nx*2].re*wave[dw*2].re - v[nx*2].im*wave[dw*2].im;
1096 r2 = v[nx*2].re*wave[dw*2].im + v[nx*2].im*wave[dw*2].re;
1097 r1 = r0 + i2; i1 = i0 + r2;
1099 r2 = icv_sin_120*(i0 - r2); i2 = icv_sin_120*(i2 - r0);
1100 r0 = v[0].re; i0 = v[0].im;
1101 v[0].re = (float)(r0 + r1); v[0].im = (float)(i0 + i1);
1102 r0 -= 0.5*r1; i0 -= 0.5*i1;
1103 v[nx].re = (float)(r0 + r2); v[nx].im = (float)(i0 + i2);
1104 v[nx*2].re = (float)(r0 - r2); v[nx*2].im = (float)(i0 - i2);
1108 else if( factor == 5 )
1111 for( i = 0; i < n0; i += n )
1113 for( j = 0, dw = 0; j < nx; j++, dw += dw0 )
1115 CvComplex32f* v0 = dst + i + j;
1116 CvComplex32f* v1 = v0 + nx*2;
1117 CvComplex32f* v2 = v1 + nx*2;
1119 double r0, i0, r1, i1, r2, i2, r3, i3, r4, i4, r5, i5;
1121 r3 = v0[nx].re*wave[dw].re - v0[nx].im*wave[dw].im;
1122 i3 = v0[nx].re*wave[dw].im + v0[nx].im*wave[dw].re;
1123 r2 = v2[0].re*wave[dw*4].re - v2[0].im*wave[dw*4].im;
1124 i2 = v2[0].re*wave[dw*4].im + v2[0].im*wave[dw*4].re;
1126 r1 = r3 + r2; i1 = i3 + i2;
1129 r4 = v1[nx].re*wave[dw*3].re - v1[nx].im*wave[dw*3].im;
1130 i4 = v1[nx].re*wave[dw*3].im + v1[nx].im*wave[dw*3].re;
1131 r0 = v1[0].re*wave[dw*2].re - v1[0].im*wave[dw*2].im;
1132 i0 = v1[0].re*wave[dw*2].im + v1[0].im*wave[dw*2].re;
1134 r2 = r4 + r0; i2 = i4 + i0;
1137 r0 = v0[0].re; i0 = v0[0].im;
1138 r5 = r1 + r2; i5 = i1 + i2;
1140 v0[0].re = (float)(r0 + r5); v0[0].im = (float)(i0 + i5);
1142 r0 -= 0.25*r5; i0 -= 0.25*i5;
1143 r1 = icv_fft5_2*(r1 - r2); i1 = icv_fft5_2*(i1 - i2);
1144 r2 = -icv_fft5_3*(i3 + i4); i2 = icv_fft5_3*(r3 + r4);
1146 i3 *= -icv_fft5_5; r3 *= icv_fft5_5;
1147 i4 *= -icv_fft5_4; r4 *= icv_fft5_4;
1149 r5 = r2 + i3; i5 = i2 + r3;
1152 r3 = r0 + r1; i3 = i0 + i1;
1155 v0[nx].re = (float)(r3 + r2); v0[nx].im = (float)(i3 + i2);
1156 v2[0].re = (float)(r3 - r2); v2[0].im = (float)(i3 - i2);
1158 v1[0].re = (float)(r0 + r5); v1[0].im = (float)(i0 + i5);
1159 v1[nx].re = (float)(r0 - r5); v1[nx].im = (float)(i0 - i5);
1165 // radix-"factor" - an odd number
1166 int p, q, factor2 = (factor - 1)/2;
1167 int d, dd, dw_f = tab_size/factor;
1168 CvComplex32f* a = buf;
1169 CvComplex32f* b = buf + factor2;
1171 for( i = 0; i < n0; i += n )
1173 for( j = 0, dw = 0; j < nx; j++, dw += dw0 )
1175 CvComplex32f* v = dst + i + j;
1176 CvComplex32f v_0 = v[0];
1177 CvComplex64f vn_0(v_0);
1181 for( p = 1, k = nx; p <= factor2; p++, k += nx )
1183 double r0 = v[k].re + v[n-k].re;
1184 double i0 = v[k].im - v[n-k].im;
1185 double r1 = v[k].re - v[n-k].re;
1186 double i1 = v[k].im + v[n-k].im;
1188 vn_0.re += r0; vn_0.im += i1;
1189 a[p-1].re = (float)r0; a[p-1].im = (float)i0;
1190 b[p-1].re = (float)r1; b[p-1].im = (float)i1;
1195 const CvComplex32f* wave_ = wave + dw*factor;
1198 for( p = 1, k = nx; p <= factor2; p++, k += nx, d += dw )
1200 double r2 = v[k].re*wave[d].re - v[k].im*wave[d].im;
1201 double i2 = v[k].re*wave[d].im + v[k].im*wave[d].re;
1203 double r1 = v[n-k].re*wave_[-d].re - v[n-k].im*wave_[-d].im;
1204 double i1 = v[n-k].re*wave_[-d].im + v[n-k].im*wave_[-d].re;
1206 double r0 = r2 + r1;
1207 double i0 = i2 - i1;
1211 vn_0.re += r0; vn_0.im += i1;
1212 a[p-1].re = (float)r0; a[p-1].im = (float)i0;
1213 b[p-1].re = (float)r1; b[p-1].im = (float)i1;
1217 v[0].re = (float)vn_0.re;
1218 v[0].im = (float)vn_0.im;
1220 for( p = 1, k = nx; p <= factor2; p++, k += nx )
1222 CvComplex64f s0(v_0), s1 = s0;
1225 for( q = 0; q < factor2; q++ )
1227 double r0 = wave[d].re * a[q].re;
1228 double i0 = wave[d].im * a[q].im;
1229 double r1 = wave[d].re * b[q].im;
1230 double i1 = wave[d].im * b[q].re;
1232 s1.re += r0 + i0; s0.re += r0 - i0;
1233 s1.im += r1 - i1; s0.im += r1 + i1;
1236 d -= -(d >= tab_size) & tab_size;
1239 v[k].re = (float)s0.re;
1240 v[k].im = (float)s0.im;
1241 v[n-k].re = (float)s1.re;
1242 v[n-k].im = (float)s1.im;
1249 if( fabs(scale - 1.) > DBL_EPSILON )
1251 double re_scale = scale, im_scale = scale;
1253 im_scale = -im_scale;
1255 for( i = 0; i < n0; i++ )
1257 double t0 = dst[i].re*re_scale;
1258 double t1 = dst[i].im*im_scale;
1259 dst[i].re = (float)t0;
1260 dst[i].im = (float)t1;
1265 for( i = 0; i <= n0 - 2; i += 2 )
1267 double t0 = -dst[i].im;
1268 double t1 = -dst[i+1].im;
1269 dst[i].im = (float)t0;
1270 dst[i+1].im = (float)t1;
1274 dst[n0-1].im = -dst[n0-1].im;
1281 /* FFT of real vector
1282 output vector format:
1283 re(0), re(1), im(1), ... , re(n/2-1), im((n+1)/2-1) [, re((n+1)/2)] OR ...
1284 re(0), 0, re(1), im(1), ..., re(n/2-1), im((n+1)/2-1) [, re((n+1)/2), 0] */
1285 #define ICV_REAL_DFT( flavor, datatype ) \
1286 static CvStatus CV_STDCALL \
1287 icvRealDFT_##flavor( const datatype* src, datatype* dst, \
1288 int n, int nf, int* factors, const int* itab, \
1289 const CvComplex##flavor* wave, int tab_size, \
1290 const void* spec, CvComplex##flavor* buf, \
1291 int flags, double scale ) \
1293 int complex_output = (flags & ICV_DFT_COMPLEX_INPUT_OR_OUTPUT) != 0;\
1294 int j, n2 = n >> 1; \
1295 dst += complex_output; \
1299 icvDFTFwd_RToPack_##flavor##_p( src, dst, spec, buf ); \
1303 assert( tab_size == n ); \
1307 dst[0] = (datatype)(src[0]*scale); \
1311 double t = (src[0] + src[1])*scale; \
1312 dst[1] = (datatype)((src[0] - src[1])*scale); \
1313 dst[0] = (datatype)t; \
1317 dst -= complex_output; \
1318 CvComplex##flavor* _dst = (CvComplex##flavor*)dst; \
1319 _dst[0].re = (datatype)(src[0]*scale); \
1321 for( j = 1; j < n; j += 2 ) \
1323 double t0 = src[itab[j]]*scale; \
1324 double t1 = src[itab[j+1]]*scale; \
1325 _dst[j].re = (datatype)t0; \
1327 _dst[j+1].re = (datatype)t1; \
1330 icvDFT_##flavor##c( _dst, _dst, n, nf, factors, itab, wave, \
1331 tab_size, 0, buf, ICV_DFT_NO_PERMUTE, 1. ); \
1332 if( !complex_output ) \
1339 double h1_re, h1_im, h2_re, h2_im; \
1340 double scale2 = scale*0.5; \
1343 icvDFT_##flavor##c( (CvComplex##flavor*)src, \
1344 (CvComplex##flavor*)dst, n2, \
1345 nf - (factors[0] == 1), \
1346 factors + (factors[0] == 1), \
1347 itab, wave, tab_size, 0, buf, 0, 1. ); \
1350 t = dst[0] - dst[1]; \
1351 dst[0] = (datatype)((dst[0] + dst[1])*scale); \
1352 dst[1] = (datatype)(t*scale); \
1356 dst[n-1] = dst[1]; \
1358 for( j = 2, wave++; j < n2; j += 2, wave++ ) \
1361 h2_re = scale2*(dst[j+1] + t); \
1362 h2_im = scale2*(dst[n-j] - dst[j]); \
1365 h1_re = scale2*(dst[j] + dst[n-j]); \
1366 h1_im = scale2*(dst[j+1] - t); \
1369 t = h2_re*wave->re - h2_im*wave->im; \
1370 h2_im = h2_re*wave->im + h2_im*wave->re; \
1374 dst[j-1] = (datatype)(h1_re + h2_re); \
1375 dst[n-j-1] = (datatype)(h1_re - h2_re); \
1376 dst[j] = (datatype)(h1_im + h2_im); \
1377 dst[n-j] = (datatype)(h2_im - h1_im); \
1382 dst[n2-1] = (datatype)(t0*scale); \
1383 dst[n2] = (datatype)(-t*scale); \
1388 if( complex_output ) \
1392 if( (n & 1) == 0 ) \
1400 /* Inverse FFT of complex conjugate-symmetric vector
1401 input vector format:
1402 re[0], re[1], im[1], ... , re[n/2-1], im[n/2-1], re[n/2] OR
1403 re(0), 0, re(1), im(1), ..., re(n/2-1), im((n+1)/2-1) [, re((n+1)/2), 0] */
1404 #define ICV_CCS_IDFT( flavor, datatype ) \
1405 static CvStatus CV_STDCALL \
1406 icvCCSIDFT_##flavor( const datatype* src, datatype* dst, \
1407 int n, int nf, int* factors, const int* itab, \
1408 const CvComplex##flavor* wave, int tab_size, \
1409 const void* spec, CvComplex##flavor* buf, \
1410 int flags, double scale ) \
1412 int complex_input = (flags & ICV_DFT_COMPLEX_INPUT_OR_OUTPUT) != 0; \
1413 int j, k, n2 = (n+1) >> 1; \
1414 double save_s1 = 0.; \
1415 double t0, t1, t2, t3, t; \
1417 assert( tab_size == n ); \
1419 if( complex_input ) \
1421 assert( src != dst ); \
1423 ((datatype*)src)[1] = src[0]; \
1429 icvDFTInv_PackToR_##flavor##_p( src, dst, spec, buf ); \
1435 dst[0] = (datatype)(src[0]*scale); \
1439 t = (src[0] + src[1])*scale; \
1440 dst[1] = (datatype)((src[0] - src[1])*scale); \
1441 dst[0] = (datatype)t; \
1445 CvComplex##flavor* _src = (CvComplex##flavor*)(src-1); \
1446 CvComplex##flavor* _dst = (CvComplex##flavor*)dst; \
1448 _dst[0].re = src[0]; \
1450 for( j = 1; j < n2; j++ ) \
1452 int k0 = itab[j], k1 = itab[n-j]; \
1453 t0 = _src[j].re; t1 = _src[j].im; \
1454 _dst[k0].re = (datatype)t0; _dst[k0].im = (datatype)-t1; \
1455 _dst[k1].re = (datatype)t0; _dst[k1].im = (datatype)t1; \
1458 icvDFT_##flavor##c( _dst, _dst, n, nf, factors, itab, wave, \
1459 tab_size, 0, buf, ICV_DFT_NO_PERMUTE, 1. ); \
1460 dst[0] = (datatype)(dst[0]*scale); \
1461 for( j = 1; j < n; j += 2 ) \
1463 t0 = dst[j*2]*scale; \
1464 t1 = dst[j*2+2]*scale; \
1465 dst[j] = (datatype)t0; \
1466 dst[j+1] = (datatype)t1; \
1471 int inplace = src == dst; \
1472 const CvComplex##flavor* w = wave; \
1475 t0 = (src[0] + src[n-1]); \
1476 t1 = (src[n-1] - src[0]); \
1477 dst[0] = (datatype)t0; \
1478 dst[1] = (datatype)t1; \
1480 for( j = 2, w++; j < n2; j += 2, w++ ) \
1482 double h1_re, h1_im, h2_re, h2_im; \
1484 h1_re = (t + src[n-j-1]); \
1485 h1_im = (src[j] - src[n-j]); \
1487 h2_re = (t - src[n-j-1]); \
1488 h2_im = (src[j] + src[n-j]); \
1490 t = h2_re*w->re + h2_im*w->im; \
1491 h2_im = h2_im*w->re - h2_re*w->im; \
1495 t0 = h1_re - h2_im; \
1496 t1 = -h1_im - h2_re; \
1497 t2 = h1_re + h2_im; \
1498 t3 = h1_im - h2_re; \
1502 dst[j] = (datatype)t0; \
1503 dst[j+1] = (datatype)t1; \
1504 dst[n-j] = (datatype)t2; \
1505 dst[n-j+1]= (datatype)t3; \
1511 dst[k] = (datatype)t0; \
1512 dst[k+1] = (datatype)t1; \
1514 dst[k] = (datatype)t2; \
1515 dst[k+1]= (datatype)t3; \
1526 dst[n2] = (datatype)t0; \
1527 dst[n2+1] = (datatype)t1; \
1532 dst[k*2] = (datatype)t0; \
1533 dst[k*2+1] = (datatype)t1; \
1538 icvDFT_##flavor##c( (CvComplex##flavor*)dst, \
1539 (CvComplex##flavor*)dst, n2, \
1540 nf - (factors[0] == 1), \
1541 factors + (factors[0] == 1), itab, \
1542 wave, tab_size, 0, buf, \
1543 inplace ? 0 : ICV_DFT_NO_PERMUTE, 1. ); \
1546 for( j = 0; j < n; j += 2 ) \
1548 t0 = dst[j]*scale; \
1549 t1 = dst[j+1]*(-scale); \
1550 dst[j] = (datatype)t0; \
1551 dst[j+1] = (datatype)t1; \
1556 if( complex_input ) \
1557 ((datatype*)src)[0] = (datatype)save_s1; \
1563 ICV_REAL_DFT( 64f, double )
1564 ICV_CCS_IDFT( 64f, double )
1565 ICV_REAL_DFT( 32f, float )
1566 ICV_CCS_IDFT( 32f, float )
1570 icvCopyColumn( const uchar* _src, int src_step,
1571 uchar* _dst, int dst_step,
1572 int len, int elem_size )
1575 const int* src = (const int*)_src;
1576 int* dst = (int*)_dst;
1577 src_step /= sizeof(src[0]);
1578 dst_step /= sizeof(dst[0]);
1580 if( elem_size == sizeof(int) )
1582 for( i = 0; i < len; i++, src += src_step, dst += dst_step )
1585 else if( elem_size == sizeof(int)*2 )
1587 for( i = 0; i < len; i++, src += src_step, dst += dst_step )
1589 t0 = src[0]; t1 = src[1];
1590 dst[0] = t0; dst[1] = t1;
1593 else if( elem_size == sizeof(int)*4 )
1595 for( i = 0; i < len; i++, src += src_step, dst += dst_step )
1597 t0 = src[0]; t1 = src[1];
1598 dst[0] = t0; dst[1] = t1;
1599 t0 = src[2]; t1 = src[3];
1600 dst[2] = t0; dst[3] = t1;
1607 icvCopyFrom2Columns( const uchar* _src, int src_step,
1608 uchar* _dst0, uchar* _dst1,
1609 int len, int elem_size )
1612 const int* src = (const int*)_src;
1613 int* dst0 = (int*)_dst0;
1614 int* dst1 = (int*)_dst1;
1615 src_step /= sizeof(src[0]);
1617 if( elem_size == sizeof(int) )
1619 for( i = 0; i < len; i++, src += src_step )
1621 t0 = src[0]; t1 = src[1];
1622 dst0[i] = t0; dst1[i] = t1;
1625 else if( elem_size == sizeof(int)*2 )
1627 for( i = 0; i < len*2; i += 2, src += src_step )
1629 t0 = src[0]; t1 = src[1];
1630 dst0[i] = t0; dst0[i+1] = t1;
1631 t0 = src[2]; t1 = src[3];
1632 dst1[i] = t0; dst1[i+1] = t1;
1635 else if( elem_size == sizeof(int)*4 )
1637 for( i = 0; i < len*4; i += 4, src += src_step )
1639 t0 = src[0]; t1 = src[1];
1640 dst0[i] = t0; dst0[i+1] = t1;
1641 t0 = src[2]; t1 = src[3];
1642 dst0[i+2] = t0; dst0[i+3] = t1;
1643 t0 = src[4]; t1 = src[5];
1644 dst1[i] = t0; dst1[i+1] = t1;
1645 t0 = src[6]; t1 = src[7];
1646 dst1[i+2] = t0; dst1[i+3] = t1;
1653 icvCopyTo2Columns( const uchar* _src0, const uchar* _src1,
1654 uchar* _dst, int dst_step,
1655 int len, int elem_size )
1658 const int* src0 = (const int*)_src0;
1659 const int* src1 = (const int*)_src1;
1660 int* dst = (int*)_dst;
1661 dst_step /= sizeof(dst[0]);
1663 if( elem_size == sizeof(int) )
1665 for( i = 0; i < len; i++, dst += dst_step )
1667 t0 = src0[i]; t1 = src1[i];
1668 dst[0] = t0; dst[1] = t1;
1671 else if( elem_size == sizeof(int)*2 )
1673 for( i = 0; i < len*2; i += 2, dst += dst_step )
1675 t0 = src0[i]; t1 = src0[i+1];
1676 dst[0] = t0; dst[1] = t1;
1677 t0 = src1[i]; t1 = src1[i+1];
1678 dst[2] = t0; dst[3] = t1;
1681 else if( elem_size == sizeof(int)*4 )
1683 for( i = 0; i < len*4; i += 4, dst += dst_step )
1685 t0 = src0[i]; t1 = src0[i+1];
1686 dst[0] = t0; dst[1] = t1;
1687 t0 = src0[i+2]; t1 = src0[i+3];
1688 dst[2] = t0; dst[3] = t1;
1689 t0 = src1[i]; t1 = src1[i+1];
1690 dst[4] = t0; dst[5] = t1;
1691 t0 = src1[i+2]; t1 = src1[i+3];
1692 dst[6] = t0; dst[7] = t1;
1699 icvExpandCCS( uchar* _ptr, int len, int elem_size )
1703 memcpy( _ptr, _ptr + elem_size, elem_size );
1704 memset( _ptr + elem_size, 0, elem_size );
1705 if( (len & 1) == 0 )
1706 memset( _ptr + (len+1)*elem_size, 0, elem_size );
1708 if( elem_size == sizeof(float) )
1710 CvComplex32f* ptr = (CvComplex32f*)_ptr;
1712 for( i = 1; i < (len+1)/2; i++ )
1722 CvComplex64f* ptr = (CvComplex64f*)_ptr;
1724 for( i = 1; i < (len+1)/2; i++ )
1735 typedef CvStatus (CV_STDCALL *CvDFTFunc)(
1736 const void* src, void* dst, int n, int nf, int* factors,
1737 const int* itab, const void* wave, int tab_size,
1738 const void* spec, void* buf, int inv, double scale );
1741 cvDFT( const CvArr* srcarr, CvArr* dstarr, int flags, int nonzero_rows )
1743 static CvDFTFunc dft_tbl[6];
1744 static int inittab = 0;
1747 int local_alloc = 1;
1749 void *spec_c = 0, *spec_r = 0, *spec = 0;
1751 CV_FUNCNAME( "cvDFT" );
1755 int prev_len = 0, buf_size = 0, stage = 0;
1756 int nf = 0, inv = (flags & CV_DXT_INVERSE) != 0;
1757 int real_transform = 0;
1758 CvMat *src = (CvMat*)srcarr, *dst = (CvMat*)dstarr;
1759 CvMat srcstub, dststub, *src0;
1760 int complex_elem_size, elem_size;
1761 int factors[34], inplace_transform = 0;
1762 int ipp_norm_flag = 0;
1766 dft_tbl[0] = (CvDFTFunc)icvDFT_32fc;
1767 dft_tbl[1] = (CvDFTFunc)icvRealDFT_32f;
1768 dft_tbl[2] = (CvDFTFunc)icvCCSIDFT_32f;
1769 dft_tbl[3] = (CvDFTFunc)icvDFT_64fc;
1770 dft_tbl[4] = (CvDFTFunc)icvRealDFT_64f;
1771 dft_tbl[5] = (CvDFTFunc)icvCCSIDFT_64f;
1775 if( !CV_IS_MAT( src ))
1778 CV_CALL( src = cvGetMat( src, &srcstub, &coi ));
1781 CV_ERROR( CV_BadCOI, "" );
1784 if( !CV_IS_MAT( dst ))
1787 CV_CALL( dst = cvGetMat( dst, &dststub, &coi ));
1790 CV_ERROR( CV_BadCOI, "" );
1794 elem_size = CV_ELEM_SIZE1(src->type);
1795 complex_elem_size = elem_size*2;
1797 // check types and sizes
1798 if( !CV_ARE_DEPTHS_EQ(src, dst) )
1799 CV_ERROR( CV_StsUnmatchedFormats, "" );
1801 depth = CV_MAT_DEPTH(src->type);
1802 if( depth < CV_32F )
1803 CV_ERROR( CV_StsUnsupportedFormat,
1804 "Only 32fC1, 32fC2, 64fC1 and 64fC2 formats are supported" );
1806 if( CV_ARE_CNS_EQ(src, dst) )
1808 if( CV_MAT_CN(src->type) > 2 )
1809 CV_ERROR( CV_StsUnsupportedFormat,
1810 "Only 32fC1, 32fC2, 64fC1 and 64fC2 formats are supported" );
1812 if( !CV_ARE_SIZES_EQ(src, dst) )
1813 CV_ERROR( CV_StsUnmatchedSizes, "" );
1814 real_transform = CV_MAT_CN(src->type) == 1;
1815 if( !real_transform )
1816 elem_size = complex_elem_size;
1818 else if( !inv && CV_MAT_CN(src->type) == 1 && CV_MAT_CN(dst->type) == 2 )
1820 if( (src->cols != 1 || dst->cols != 1 ||
1821 src->rows/2+1 != dst->rows && src->rows != dst->rows) &&
1822 (src->cols/2+1 != dst->cols || src->rows != dst->rows) )
1823 CV_ERROR( CV_StsUnmatchedSizes, "" );
1826 else if( inv && CV_MAT_CN(src->type) == 2 && CV_MAT_CN(dst->type) == 1 )
1828 if( (src->cols != 1 || dst->cols != 1 ||
1829 dst->rows/2+1 != src->rows && src->rows != dst->rows) &&
1830 (dst->cols/2+1 != src->cols || src->rows != dst->rows) )
1831 CV_ERROR( CV_StsUnmatchedSizes, "" );
1835 CV_ERROR( CV_StsUnmatchedFormats,
1836 "Incorrect or unsupported combination of input & output formats" );
1838 if( src->cols == 1 && nonzero_rows > 0 )
1839 CV_ERROR( CV_StsNotImplemented,
1840 "This mode (using nonzero_rows with a single-column matrix) breaks the function logic, so it is prohibited.\n"
1841 "For fast convolution/correlation use 2-column matrix or single-row matrix instead" );
1843 // determine, which transform to do first - row-wise
1844 // (stage 0) or column-wise (stage 1) transform
1845 if( !(flags & CV_DXT_ROWS) && src->rows > 1 &&
1846 (src->cols == 1 && !CV_IS_MAT_CONT(src->type & dst->type) ||
1847 src->cols > 1 && inv && real_transform) )
1850 ipp_norm_flag = !(flags & CV_DXT_SCALE) ? 8 : (flags & CV_DXT_INVERSE) ? 2 : 1;
1858 int i, len, count, sz = 0;
1859 int use_buf = 0, odd_real = 0;
1862 if( stage == 0 ) // row-wise transform
1864 len = !inv ? src->cols : dst->cols;
1866 if( len == 1 && !(flags & CV_DXT_ROWS) )
1868 len = !inv ? src->rows : dst->rows;
1871 odd_real = real_transform && (len & 1);
1876 count = !inv ? src0->cols : dst->cols;
1877 sz = 2*len*complex_elem_size;
1881 if( len*count >= 64 && icvDFTInitAlloc_R_32f_p != 0 ) // use IPP DFT if available
1885 if( real_transform && stage == 0 )
1887 if( depth == CV_32F )
1890 IPPI_CALL( icvDFTFree_R_32f_p( spec_r ));
1891 IPPI_CALL( icvDFTInitAlloc_R_32f_p(
1892 &spec_r, len, ipp_norm_flag, cvAlgHintNone ));
1893 IPPI_CALL( icvDFTGetBufSize_R_32f_p( spec_r, &ipp_sz ));
1898 IPPI_CALL( icvDFTFree_R_64f_p( spec_r ));
1899 IPPI_CALL( icvDFTInitAlloc_R_64f_p(
1900 &spec_r, len, ipp_norm_flag, cvAlgHintNone ));
1901 IPPI_CALL( icvDFTGetBufSize_R_64f_p( spec_r, &ipp_sz ));
1907 if( depth == CV_32F )
1910 IPPI_CALL( icvDFTFree_C_32fc_p( spec_c ));
1911 IPPI_CALL( icvDFTInitAlloc_C_32fc_p(
1912 &spec_c, len, ipp_norm_flag, cvAlgHintNone ));
1913 IPPI_CALL( icvDFTGetBufSize_C_32fc_p( spec_c, &ipp_sz ));
1918 IPPI_CALL( icvDFTFree_C_64fc_p( spec_c ));
1919 IPPI_CALL( icvDFTInitAlloc_C_64fc_p(
1920 &spec_c, len, ipp_norm_flag, cvAlgHintNone ));
1921 IPPI_CALL( icvDFTGetBufSize_C_64fc_p( spec_c, &ipp_sz ));
1930 if( len != prev_len )
1931 nf = icvDFTFactorize( len, factors );
1933 inplace_transform = factors[0] == factors[nf-1];
1934 sz += len*(complex_elem_size + sizeof(int));
1935 i = nf > 1 && (factors[0] & 1) == 0;
1936 if( (factors[i] & 1) != 0 && factors[i] > 5 )
1937 sz += (factors[i]+1)*complex_elem_size;
1939 if( stage == 0 && (src->data.ptr == dst->data.ptr && !inplace_transform || odd_real) ||
1940 stage == 1 && !inplace_transform )
1943 sz += len*complex_elem_size;
1949 prev_len = 0; // because we release the buffer,
1950 // force recalculation of
1951 // twiddle factors and permutation table
1952 if( !local_alloc && buffer )
1954 if( sz <= CV_MAX_LOCAL_DFT_SIZE )
1956 buf_size = sz = CV_MAX_LOCAL_DFT_SIZE;
1957 buffer = cvStackAlloc(sz + 32);
1962 CV_CALL( buffer = cvAlloc(sz + 32) );
1968 ptr = (uchar*)buffer;
1972 ptr += len*complex_elem_size;
1974 ptr = (uchar*)cvAlignPtr( ptr + len*sizeof(int), 16 );
1976 if( len != prev_len || (!inplace_transform && inv && real_transform))
1977 icvDFTInit( len, nf, factors, itab, complex_elem_size,
1978 wave, stage == 0 && inv && real_transform );
1979 // otherwise reuse the tables calculated on the previous stage
1985 int dptr_offset = 0;
1986 int dst_full_len = len*elem_size;
1987 int _flags = inv + (CV_MAT_CN(src->type) != CV_MAT_CN(dst->type) ?
1988 ICV_DFT_COMPLEX_INPUT_OR_OUTPUT : 0);
1992 ptr += len*complex_elem_size;
1993 if( odd_real && !inv && len > 1 &&
1994 !(_flags & ICV_DFT_COMPLEX_INPUT_OR_OUTPUT))
1995 dptr_offset = elem_size;
1998 if( !inv && (_flags & ICV_DFT_COMPLEX_INPUT_OR_OUTPUT) )
1999 dst_full_len += (len & 1) ? elem_size : complex_elem_size;
2001 dft_func = dft_tbl[(!real_transform ? 0 : !inv ? 1 : 2) + (depth == CV_64F)*3];
2003 if( count > 1 && !(flags & CV_DXT_ROWS) && (!inv || !real_transform) )
2005 else if( flags & CV_DXT_SCALE )
2006 scale = 1./(len * (flags & CV_DXT_ROWS ? 1 : count));
2008 if( nonzero_rows <= 0 || nonzero_rows > count )
2009 nonzero_rows = count;
2011 for( i = 0; i < nonzero_rows; i++ )
2013 uchar* sptr = src->data.ptr + i*src->step;
2014 uchar* dptr0 = dst->data.ptr + i*dst->step;
2015 uchar* dptr = dptr0;
2020 dft_func( sptr, dptr, len, nf, factors, itab, wave, len, spec, ptr, _flags, scale );
2022 memcpy( dptr0, dptr + dptr_offset, dst_full_len );
2025 for( ; i < count; i++ )
2027 uchar* dptr0 = dst->data.ptr + i*dst->step;
2028 memset( dptr0, 0, dst_full_len );
2037 int a = 0, b = count;
2038 uchar *buf0, *buf1, *dbuf0, *dbuf1;
2039 uchar* sptr0 = src->data.ptr;
2040 uchar* dptr0 = dst->data.ptr;
2042 ptr += len*complex_elem_size;
2044 ptr += len*complex_elem_size;
2045 dbuf0 = buf0, dbuf1 = buf1;
2051 ptr += len*complex_elem_size;
2054 dft_func = dft_tbl[(depth == CV_64F)*3];
2056 if( real_transform && inv && src->cols > 1 )
2058 else if( flags & CV_DXT_SCALE )
2059 scale = 1./(len * count);
2061 if( real_transform )
2065 even = (count & 1) == 0;
2069 memset( buf0, 0, len*complex_elem_size );
2070 icvCopyColumn( sptr0, src->step, buf0, complex_elem_size, len, elem_size );
2071 sptr0 += CV_MAT_CN(dst->type)*elem_size;
2074 memset( buf1, 0, len*complex_elem_size );
2075 icvCopyColumn( sptr0 + (count-2)*elem_size, src->step,
2076 buf1, complex_elem_size, len, elem_size );
2079 else if( CV_MAT_CN(src->type) == 1 )
2081 icvCopyColumn( sptr0, src->step, buf0 + elem_size, elem_size, len, elem_size );
2082 icvExpandCCS( buf0 + elem_size, len, elem_size );
2085 icvCopyColumn( sptr0 + (count-1)*elem_size, src->step,
2086 buf1 + elem_size, elem_size, len, elem_size );
2087 icvExpandCCS( buf1 + elem_size, len, elem_size );
2093 icvCopyColumn( sptr0, src->step, buf0, complex_elem_size, len, complex_elem_size );
2094 //memcpy( buf0 + elem_size, buf0, elem_size );
2095 //icvExpandCCS( buf0 + elem_size, len, elem_size );
2098 icvCopyColumn( sptr0 + b*complex_elem_size, src->step,
2099 buf1, complex_elem_size, len, complex_elem_size );
2100 //memcpy( buf0 + elem_size, buf0, elem_size );
2101 //icvExpandCCS( buf0 + elem_size, len, elem_size );
2103 sptr0 += complex_elem_size;
2107 IPPI_CALL( dft_func( buf1, dbuf1, len, nf, factors, itab,
2108 wave, len, spec, ptr, inv, scale ));
2109 IPPI_CALL( dft_func( buf0, dbuf0, len, nf, factors, itab,
2110 wave, len, spec, ptr, inv, scale ));
2112 if( CV_MAT_CN(dst->type) == 1 )
2116 // copy the half of output vector to the first/last column.
2117 // before doing that, defgragment the vector
2118 memcpy( dbuf0 + elem_size, dbuf0, elem_size );
2119 icvCopyColumn( dbuf0 + elem_size, elem_size, dptr0,
2120 dst->step, len, elem_size );
2123 memcpy( dbuf1 + elem_size, dbuf1, elem_size );
2124 icvCopyColumn( dbuf1 + elem_size, elem_size,
2125 dptr0 + (count-1)*elem_size,
2126 dst->step, len, elem_size );
2132 // copy the real part of the complex vector to the first/last column
2133 icvCopyColumn( dbuf0, complex_elem_size, dptr0, dst->step, len, elem_size );
2135 icvCopyColumn( dbuf1, complex_elem_size, dptr0 + (count-1)*elem_size,
2136 dst->step, len, elem_size );
2143 icvCopyColumn( dbuf0, complex_elem_size, dptr0,
2144 dst->step, len, complex_elem_size );
2146 icvCopyColumn( dbuf1, complex_elem_size,
2147 dptr0 + b*complex_elem_size,
2148 dst->step, len, complex_elem_size );
2149 dptr0 += complex_elem_size;
2153 for( i = a; i < b; i += 2 )
2157 icvCopyFrom2Columns( sptr0, src->step, buf0, buf1, len, complex_elem_size );
2158 IPPI_CALL( dft_func( buf1, dbuf1, len, nf, factors, itab,
2159 wave, len, spec, ptr, inv, scale ));
2162 icvCopyColumn( sptr0, src->step, buf0, complex_elem_size, len, complex_elem_size );
2164 IPPI_CALL( dft_func( buf0, dbuf0, len, nf, factors, itab,
2165 wave, len, spec, ptr, inv, scale ));
2168 icvCopyTo2Columns( dbuf0, dbuf1, dptr0, dst->step, len, complex_elem_size );
2170 icvCopyColumn( dbuf0, complex_elem_size, dptr0, dst->step, len, complex_elem_size );
2171 sptr0 += 2*complex_elem_size;
2172 dptr0 += 2*complex_elem_size;
2183 if( buffer && !local_alloc )
2188 if( depth == CV_32F )
2189 icvDFTFree_C_32fc_p( spec_c );
2191 icvDFTFree_C_64fc_p( spec_c );
2196 if( depth == CV_32F )
2197 icvDFTFree_R_32f_p( spec_r );
2199 icvDFTFree_R_64f_p( spec_r );
2205 cvMulSpectrums( const CvArr* srcAarr, const CvArr* srcBarr,
2206 CvArr* dstarr, int flags )
2208 CV_FUNCNAME( "cvMulSpectrums" );
2212 CvMat stubA, *srcA = (CvMat*)srcAarr;
2213 CvMat stubB, *srcB = (CvMat*)srcBarr;
2214 CvMat dststub, *dst = (CvMat*)dstarr;
2215 int stepA, stepB, stepC;
2216 int type, cn, is_1d;
2217 int j, j0, j1, k, rows, cols, ncols;
2219 if( !CV_IS_MAT(srcA))
2220 CV_CALL( srcA = cvGetMat( srcA, &stubA, 0 ));
2222 if( !CV_IS_MAT(srcB))
2223 CV_CALL( srcB = cvGetMat( srcB, &stubB, 0 ));
2225 if( !CV_IS_MAT(dst))
2226 CV_CALL( dst = cvGetMat( dst, &dststub, 0 ));
2228 if( !CV_ARE_TYPES_EQ( srcA, srcB ) || !CV_ARE_TYPES_EQ( srcA, dst ))
2229 CV_ERROR( CV_StsUnmatchedFormats, "" );
2231 if( !CV_ARE_SIZES_EQ( srcA, srcB ) || !CV_ARE_SIZES_EQ( srcA, dst ))
2232 CV_ERROR( CV_StsUnmatchedSizes, "" );
2234 type = CV_MAT_TYPE( dst->type );
2235 cn = CV_MAT_CN(type);
2238 is_1d = (flags & CV_DXT_ROWS) ||
2239 (rows == 1 || cols == 1 &&
2240 CV_IS_MAT_CONT( srcA->type & srcB->type & dst->type ));
2242 if( is_1d && !(flags & CV_DXT_ROWS) )
2243 cols = cols + rows - 1, rows = 1;
2246 j1 = ncols - (cols % 2 == 0 && cn == 1);
2248 if( CV_MAT_DEPTH(type) == CV_32F )
2250 float* dataA = srcA->data.fl;
2251 float* dataB = srcB->data.fl;
2252 float* dataC = dst->data.fl;
2254 stepA = srcA->step/sizeof(dataA[0]);
2255 stepB = srcB->step/sizeof(dataB[0]);
2256 stepC = dst->step/sizeof(dataC[0]);
2258 if( !is_1d && cn == 1 )
2260 for( k = 0; k < (cols % 2 ? 1 : 2); k++ )
2263 dataA += cols - 1, dataB += cols - 1, dataC += cols - 1;
2264 dataC[0] = dataA[0]*dataB[0];
2266 dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA]*dataB[(rows-1)*stepB];
2267 if( !(flags & CV_DXT_MUL_CONJ) )
2268 for( j = 1; j <= rows - 2; j += 2 )
2270 double re = (double)dataA[j*stepA]*dataB[j*stepB] -
2271 (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
2272 double im = (double)dataA[j*stepA]*dataB[(j+1)*stepB] +
2273 (double)dataA[(j+1)*stepA]*dataB[j*stepB];
2274 dataC[j*stepC] = (float)re; dataC[(j+1)*stepC] = (float)im;
2277 for( j = 1; j <= rows - 2; j += 2 )
2279 double re = (double)dataA[j*stepA]*dataB[j*stepB] +
2280 (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
2281 double im = (double)dataA[(j+1)*stepA]*dataB[j*stepB] -
2282 (double)dataA[j*stepA]*dataB[(j+1)*stepB];
2283 dataC[j*stepC] = (float)re; dataC[(j+1)*stepC] = (float)im;
2286 dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1;
2290 for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC )
2292 if( is_1d && cn == 1 )
2294 dataC[0] = dataA[0]*dataB[0];
2296 dataC[j1] = dataA[j1]*dataB[j1];
2299 if( !(flags & CV_DXT_MUL_CONJ) )
2300 for( j = j0; j < j1; j += 2 )
2302 double re = (double)dataA[j]*dataB[j] - (double)dataA[j+1]*dataB[j+1];
2303 double im = (double)dataA[j+1]*dataB[j] + (double)dataA[j]*dataB[j+1];
2304 dataC[j] = (float)re; dataC[j+1] = (float)im;
2307 for( j = j0; j < j1; j += 2 )
2309 double re = (double)dataA[j]*dataB[j] + (double)dataA[j+1]*dataB[j+1];
2310 double im = (double)dataA[j+1]*dataB[j] - (double)dataA[j]*dataB[j+1];
2311 dataC[j] = (float)re; dataC[j+1] = (float)im;
2315 else if( CV_MAT_DEPTH(type) == CV_64F )
2317 double* dataA = srcA->data.db;
2318 double* dataB = srcB->data.db;
2319 double* dataC = dst->data.db;
2321 stepA = srcA->step/sizeof(dataA[0]);
2322 stepB = srcB->step/sizeof(dataB[0]);
2323 stepC = dst->step/sizeof(dataC[0]);
2325 if( !is_1d && cn == 1 )
2327 for( k = 0; k < (cols % 2 ? 1 : 2); k++ )
2330 dataA += cols - 1, dataB += cols - 1, dataC += cols - 1;
2331 dataC[0] = dataA[0]*dataB[0];
2333 dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA]*dataB[(rows-1)*stepB];
2334 if( !(flags & CV_DXT_MUL_CONJ) )
2335 for( j = 1; j <= rows - 2; j += 2 )
2337 double re = dataA[j*stepA]*dataB[j*stepB] -
2338 dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
2339 double im = dataA[j*stepA]*dataB[(j+1)*stepB] +
2340 dataA[(j+1)*stepA]*dataB[j*stepB];
2341 dataC[j*stepC] = re; dataC[(j+1)*stepC] = im;
2344 for( j = 1; j <= rows - 2; j += 2 )
2346 double re = dataA[j*stepA]*dataB[j*stepB] +
2347 dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
2348 double im = dataA[(j+1)*stepA]*dataB[j*stepB] -
2349 dataA[j*stepA]*dataB[(j+1)*stepB];
2350 dataC[j*stepC] = re; dataC[(j+1)*stepC] = im;
2353 dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1;
2357 for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC )
2359 if( is_1d && cn == 1 )
2361 dataC[0] = dataA[0]*dataB[0];
2363 dataC[j1] = dataA[j1]*dataB[j1];
2366 if( !(flags & CV_DXT_MUL_CONJ) )
2367 for( j = j0; j < j1; j += 2 )
2369 double re = dataA[j]*dataB[j] - dataA[j+1]*dataB[j+1];
2370 double im = dataA[j+1]*dataB[j] + dataA[j]*dataB[j+1];
2371 dataC[j] = re; dataC[j+1] = im;
2374 for( j = j0; j < j1; j += 2 )
2376 double re = dataA[j]*dataB[j] + dataA[j+1]*dataB[j+1];
2377 double im = dataA[j+1]*dataB[j] - dataA[j]*dataB[j+1];
2378 dataC[j] = re; dataC[j+1] = im;
2384 CV_ERROR( CV_StsUnsupportedFormat, "Only 32f and 64f types are supported" );
2391 /****************************************************************************************\
2392 Discrete Cosine Transform
2393 \****************************************************************************************/
2395 /* DCT is calculated using DFT, as described here:
2396 http://www.ece.utexas.edu/~bevans/courses/ee381k/lectures/09_DCT/lecture9/:
2398 #define ICV_DCT_FWD( flavor, datatype ) \
2399 static CvStatus CV_STDCALL \
2400 icvDCT_fwd_##flavor( const datatype* src, int src_step, datatype* dft_src,\
2401 datatype* dft_dst, datatype* dst, int dst_step, \
2402 int n, int nf, int* factors, const int* itab, \
2403 const CvComplex##flavor* dft_wave, \
2404 const CvComplex##flavor* dct_wave, \
2405 const void* spec, CvComplex##flavor* buf ) \
2407 int j, n2 = n >> 1; \
2409 src_step /= sizeof(src[0]); \
2410 dst_step /= sizeof(dst[0]); \
2411 datatype* dst1 = dst + (n-1)*dst_step; \
2419 for( j = 0; j < n2; j++, src += src_step*2 ) \
2421 dft_src[j] = src[0]; \
2422 dft_src[n-j-1] = src[src_step]; \
2425 icvRealDFT_##flavor( dft_src, dft_dst, n, nf, factors, \
2426 itab, dft_wave, n, spec, buf, 0, 1.0 ); \
2429 dst[0] = (datatype)(src[0]*dct_wave->re*icv_sin_45); \
2431 for( j = 1, dct_wave++; j < n2; j++, dct_wave++, \
2432 dst += dst_step, dst1 -= dst_step ) \
2434 double t0 = dct_wave->re*src[j*2-1] - dct_wave->im*src[j*2]; \
2435 double t1 = -dct_wave->im*src[j*2-1] - dct_wave->re*src[j*2]; \
2436 dst[0] = (datatype)t0; \
2437 dst1[0] = (datatype)t1; \
2440 dst[0] = (datatype)(src[n-1]*dct_wave->re); \
2445 #define ICV_DCT_INV( flavor, datatype ) \
2446 static CvStatus CV_STDCALL \
2447 icvDCT_inv_##flavor( const datatype* src, int src_step, datatype* dft_src,\
2448 datatype* dft_dst, datatype* dst, int dst_step, \
2449 int n, int nf, int* factors, const int* itab, \
2450 const CvComplex##flavor* dft_wave, \
2451 const CvComplex##flavor* dct_wave, \
2452 const void* spec, CvComplex##flavor* buf ) \
2454 int j, n2 = n >> 1; \
2456 src_step /= sizeof(src[0]); \
2457 dst_step /= sizeof(dst[0]); \
2458 const datatype* src1 = src + (n-1)*src_step; \
2466 dft_src[0] = (datatype)(src[0]*2*dct_wave->re*icv_sin_45); \
2468 for( j = 1, dct_wave++; j < n2; j++, dct_wave++, \
2469 src += src_step, src1 -= src_step ) \
2471 double t0 = dct_wave->re*src[0] - dct_wave->im*src1[0]; \
2472 double t1 = -dct_wave->im*src[0] - dct_wave->re*src1[0]; \
2473 dft_src[j*2-1] = (datatype)t0; \
2474 dft_src[j*2] = (datatype)t1; \
2477 dft_src[n-1] = (datatype)(src[0]*2*dct_wave->re); \
2478 icvCCSIDFT_##flavor( dft_src, dft_dst, n, nf, factors, itab, \
2479 dft_wave, n, spec, buf, CV_DXT_INVERSE, 1.0 ); \
2481 for( j = 0; j < n2; j++, dst += dst_step*2 ) \
2483 dst[0] = dft_dst[j]; \
2484 dst[dst_step] = dft_dst[n-j-1]; \
2490 ICV_DCT_FWD( 64f, double )
2491 ICV_DCT_INV( 64f, double )
2492 ICV_DCT_FWD( 32f, float )
2493 ICV_DCT_INV( 32f, float )
2496 icvDCTInit( int n, int elem_size, void* _wave, int inv )
2498 static const double icvDctScale[] =
2500 0.707106781186547570, 0.500000000000000000, 0.353553390593273790,
2501 0.250000000000000000, 0.176776695296636890, 0.125000000000000000,
2502 0.088388347648318447, 0.062500000000000000, 0.044194173824159223,
2503 0.031250000000000000, 0.022097086912079612, 0.015625000000000000,
2504 0.011048543456039806, 0.007812500000000000, 0.005524271728019903,
2505 0.003906250000000000, 0.002762135864009952, 0.001953125000000000,
2506 0.001381067932004976, 0.000976562500000000, 0.000690533966002488,
2507 0.000488281250000000, 0.000345266983001244, 0.000244140625000000,
2508 0.000172633491500622, 0.000122070312500000, 0.000086316745750311,
2509 0.000061035156250000, 0.000043158372875155, 0.000030517578125000
2519 assert( (n&1) == 0 );
2521 if( (n & (n - 1)) == 0 )
2524 scale = (!inv ? 2 : 1)*icvDctScale[m];
2525 w1.re = icvDxtTab[m+2][0];
2526 w1.im = -icvDxtTab[m+2][1];
2531 scale = (!inv ? 2 : 1)*sqrt(t);
2532 w1.im = sin(-CV_PI*t);
2533 w1.re = sqrt(1. - w1.im*w1.im);
2537 if( elem_size == sizeof(CvComplex64f) )
2539 CvComplex64f* wave = (CvComplex64f*)_wave;
2544 for( i = 0; i <= n; i++ )
2547 t = w.re*w1.re - w.im*w1.im;
2548 w.im = w.re*w1.im + w.im*w1.re;
2554 CvComplex32f* wave = (CvComplex32f*)_wave;
2555 assert( elem_size == sizeof(CvComplex32f) );
2557 w.re = (float)scale;
2560 for( i = 0; i <= n; i++ )
2562 wave[i].re = (float)w.re;
2563 wave[i].im = (float)w.im;
2564 t = w.re*w1.re - w.im*w1.im;
2565 w.im = w.re*w1.im + w.im*w1.re;
2572 typedef CvStatus (CV_STDCALL * CvDCTFunc)(
2573 const void* src, int src_step, void* dft_src,
2574 void* dft_dst, void* dst, int dst_step, int n,
2575 int nf, int* factors, const int* itab, const void* dft_wave,
2576 const void* dct_wave, const void* spec, void* buf );
2579 cvDCT( const CvArr* srcarr, CvArr* dstarr, int flags )
2581 static CvDCTFunc dct_tbl[4];
2582 static int inittab = 0;
2585 int local_alloc = 1;
2586 int inv = (flags & CV_DXT_INVERSE) != 0, depth = -1;
2587 void *spec_dft = 0, *spec = 0;
2589 CV_FUNCNAME( "cvDCT" );
2594 int prev_len = 0, buf_size = 0, nf = 0, stage, end_stage;
2595 CvMat *src = (CvMat*)srcarr, *dst = (CvMat*)dstarr;
2596 uchar *src_dft_buf = 0, *dst_dft_buf = 0;
2597 uchar *dft_wave = 0, *dct_wave = 0;
2600 CvMat srcstub, dststub;
2601 int complex_elem_size, elem_size;
2602 int factors[34], inplace_transform;
2608 dct_tbl[0] = (CvDCTFunc)icvDCT_fwd_32f;
2609 dct_tbl[1] = (CvDCTFunc)icvDCT_inv_32f;
2610 dct_tbl[2] = (CvDCTFunc)icvDCT_fwd_64f;
2611 dct_tbl[3] = (CvDCTFunc)icvDCT_inv_64f;
2615 if( !CV_IS_MAT( src ))
2618 CV_CALL( src = cvGetMat( src, &srcstub, &coi ));
2621 CV_ERROR( CV_BadCOI, "" );
2624 if( !CV_IS_MAT( dst ))
2627 CV_CALL( dst = cvGetMat( dst, &dststub, &coi ));
2630 CV_ERROR( CV_BadCOI, "" );
2633 depth = CV_MAT_DEPTH(src->type);
2634 elem_size = CV_ELEM_SIZE1(depth);
2635 complex_elem_size = elem_size*2;
2637 // check types and sizes
2638 if( !CV_ARE_TYPES_EQ(src, dst) )
2639 CV_ERROR( CV_StsUnmatchedFormats, "" );
2641 if( depth < CV_32F || CV_MAT_CN(src->type) != 1 )
2642 CV_ERROR( CV_StsUnsupportedFormat,
2643 "Only 32fC1 and 64fC1 formats are supported" );
2645 dct_func = dct_tbl[inv + (depth == CV_64F)*2];
2647 if( (flags & CV_DXT_ROWS) || src->rows == 1 ||
2648 src->cols == 1 && CV_IS_MAT_CONT(src->type & dst->type))
2650 stage = end_stage = 0;
2654 stage = src->cols == 1;
2658 for( ; stage <= end_stage; stage++ )
2660 uchar *sptr = src->data.ptr, *dptr = dst->data.ptr;
2661 int sstep0, sstep1, dstep0, dstep1;
2667 if( len == 1 && !(flags & CV_DXT_ROWS) )
2674 sstep1 = dstep1 = elem_size;
2682 sstep0 = dstep0 = elem_size;
2685 if( len != prev_len )
2689 if( len > 1 && (len & 1) )
2690 CV_ERROR( CV_StsNotImplemented, "Odd-size DCT\'s are not implemented" );
2693 sz += (len/2 + 1)*complex_elem_size;
2696 inplace_transform = 1;
2697 if( len*count >= 64 && icvDFTInitAlloc_R_32f_p )
2700 if( depth == CV_32F )
2703 IPPI_CALL( icvDFTFree_R_32f_p( spec_dft ));
2704 IPPI_CALL( icvDFTInitAlloc_R_32f_p( &spec_dft, len, 8, cvAlgHintNone ));
2705 IPPI_CALL( icvDFTGetBufSize_R_32f_p( spec_dft, &ipp_sz ));
2710 IPPI_CALL( icvDFTFree_R_64f_p( spec_dft ));
2711 IPPI_CALL( icvDFTInitAlloc_R_64f_p( &spec_dft, len, 8, cvAlgHintNone ));
2712 IPPI_CALL( icvDFTGetBufSize_R_64f_p( spec_dft, &ipp_sz ));
2719 sz += len*(complex_elem_size + sizeof(int)) + complex_elem_size;
2721 nf = icvDFTFactorize( len, factors );
2722 inplace_transform = factors[0] == factors[nf-1];
2724 i = nf > 1 && (factors[0] & 1) == 0;
2725 if( (factors[i] & 1) != 0 && factors[i] > 5 )
2726 sz += (factors[i]+1)*complex_elem_size;
2728 if( !inplace_transform )
2729 sz += len*elem_size;
2734 if( !local_alloc && buffer )
2736 if( sz <= CV_MAX_LOCAL_DFT_SIZE )
2738 buf_size = sz = CV_MAX_LOCAL_DFT_SIZE;
2739 buffer = cvStackAlloc(sz + 32);
2744 CV_CALL( buffer = cvAlloc(sz + 32) );
2750 ptr = (uchar*)buffer;
2754 ptr += len*complex_elem_size;
2756 ptr = (uchar*)cvAlignPtr( ptr + len*sizeof(int), 16 );
2757 icvDFTInit( len, nf, factors, itab, complex_elem_size, dft_wave, inv );
2761 ptr += (len/2 + 1)*complex_elem_size;
2762 src_dft_buf = dst_dft_buf = ptr;
2763 ptr += len*elem_size;
2764 if( !inplace_transform )
2767 ptr += len*elem_size;
2769 icvDCTInit( len, complex_elem_size, dct_wave, inv );
2774 // otherwise reuse the tables calculated on the previous stage
2776 for( i = 0; i < count; i++ )
2778 dct_func( sptr + i*sstep0, sstep1, src_dft_buf, dst_dft_buf,
2779 dptr + i*dstep0, dstep1, len, nf, factors,
2780 itab, dft_wave, dct_wave, spec, ptr );
2789 if( depth == CV_32F )
2790 icvDFTFree_R_32f_p( spec_dft );
2792 icvDFTFree_R_64f_p( spec_dft );
2795 if( buffer && !local_alloc )
2800 static const int icvOptimalDFTSize[] = {
2801 1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48,
2802 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, 128, 135, 144, 150, 160,
2803 162, 180, 192, 200, 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375,
2804 384, 400, 405, 432, 450, 480, 486, 500, 512, 540, 576, 600, 625, 640, 648, 675, 720,
2805 729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 1080, 1125, 1152, 1200,
2806 1215, 1250, 1280, 1296, 1350, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875,
2807 1920, 1944, 2000, 2025, 2048, 2160, 2187, 2250, 2304, 2400, 2430, 2500, 2560, 2592,
2808 2700, 2880, 2916, 3000, 3072, 3125, 3200, 3240, 3375, 3456, 3600, 3645, 3750, 3840,
2809 3888, 4000, 4050, 4096, 4320, 4374, 4500, 4608, 4800, 4860, 5000, 5120, 5184, 5400,
2810 5625, 5760, 5832, 6000, 6075, 6144, 6250, 6400, 6480, 6561, 6750, 6912, 7200, 7290,
2811 7500, 7680, 7776, 8000, 8100, 8192, 8640, 8748, 9000, 9216, 9375, 9600, 9720, 10000,
2812 10125, 10240, 10368, 10800, 10935, 11250, 11520, 11664, 12000, 12150, 12288, 12500,
2813 12800, 12960, 13122, 13500, 13824, 14400, 14580, 15000, 15360, 15552, 15625, 16000,
2814 16200, 16384, 16875, 17280, 17496, 18000, 18225, 18432, 18750, 19200, 19440, 19683,
2815 20000, 20250, 20480, 20736, 21600, 21870, 22500, 23040, 23328, 24000, 24300, 24576,
2816 25000, 25600, 25920, 26244, 27000, 27648, 28125, 28800, 29160, 30000, 30375, 30720,
2817 31104, 31250, 32000, 32400, 32768, 32805, 33750, 34560, 34992, 36000, 36450, 36864,
2818 37500, 38400, 38880, 39366, 40000, 40500, 40960, 41472, 43200, 43740, 45000, 46080,
2819 46656, 46875, 48000, 48600, 49152, 50000, 50625, 51200, 51840, 52488, 54000, 54675,
2820 55296, 56250, 57600, 58320, 59049, 60000, 60750, 61440, 62208, 62500, 64000, 64800,
2821 65536, 65610, 67500, 69120, 69984, 72000, 72900, 73728, 75000, 76800, 77760, 78125,
2822 78732, 80000, 81000, 81920, 82944, 84375, 86400, 87480, 90000, 91125, 92160, 93312,
2823 93750, 96000, 97200, 98304, 98415, 100000, 101250, 102400, 103680, 104976, 108000,
2824 109350, 110592, 112500, 115200, 116640, 118098, 120000, 121500, 122880, 124416, 125000,
2825 128000, 129600, 131072, 131220, 135000, 138240, 139968, 140625, 144000, 145800, 147456,
2826 150000, 151875, 153600, 155520, 156250, 157464, 160000, 162000, 163840, 164025, 165888,
2827 168750, 172800, 174960, 177147, 180000, 182250, 184320, 186624, 187500, 192000, 194400,
2828 196608, 196830, 200000, 202500, 204800, 207360, 209952, 216000, 218700, 221184, 225000,
2829 230400, 233280, 234375, 236196, 240000, 243000, 245760, 248832, 250000, 253125, 256000,
2830 259200, 262144, 262440, 270000, 273375, 276480, 279936, 281250, 288000, 291600, 294912,
2831 295245, 300000, 303750, 307200, 311040, 312500, 314928, 320000, 324000, 327680, 328050,
2832 331776, 337500, 345600, 349920, 354294, 360000, 364500, 368640, 373248, 375000, 384000,
2833 388800, 390625, 393216, 393660, 400000, 405000, 409600, 414720, 419904, 421875, 432000,
2834 437400, 442368, 450000, 455625, 460800, 466560, 468750, 472392, 480000, 486000, 491520,
2835 492075, 497664, 500000, 506250, 512000, 518400, 524288, 524880, 531441, 540000, 546750,
2836 552960, 559872, 562500, 576000, 583200, 589824, 590490, 600000, 607500, 614400, 622080,
2837 625000, 629856, 640000, 648000, 655360, 656100, 663552, 675000, 691200, 699840, 703125,
2838 708588, 720000, 729000, 737280, 746496, 750000, 759375, 768000, 777600, 781250, 786432,
2839 787320, 800000, 810000, 819200, 820125, 829440, 839808, 843750, 864000, 874800, 884736,
2840 885735, 900000, 911250, 921600, 933120, 937500, 944784, 960000, 972000, 983040, 984150,
2841 995328, 1000000, 1012500, 1024000, 1036800, 1048576, 1049760, 1062882, 1080000, 1093500,
2842 1105920, 1119744, 1125000, 1152000, 1166400, 1171875, 1179648, 1180980, 1200000,
2843 1215000, 1228800, 1244160, 1250000, 1259712, 1265625, 1280000, 1296000, 1310720,
2844 1312200, 1327104, 1350000, 1366875, 1382400, 1399680, 1406250, 1417176, 1440000,
2845 1458000, 1474560, 1476225, 1492992, 1500000, 1518750, 1536000, 1555200, 1562500,
2846 1572864, 1574640, 1594323, 1600000, 1620000, 1638400, 1640250, 1658880, 1679616,
2847 1687500, 1728000, 1749600, 1769472, 1771470, 1800000, 1822500, 1843200, 1866240,
2848 1875000, 1889568, 1920000, 1944000, 1953125, 1966080, 1968300, 1990656, 2000000,
2849 2025000, 2048000, 2073600, 2097152, 2099520, 2109375, 2125764, 2160000, 2187000,
2850 2211840, 2239488, 2250000, 2278125, 2304000, 2332800, 2343750, 2359296, 2361960,
2851 2400000, 2430000, 2457600, 2460375, 2488320, 2500000, 2519424, 2531250, 2560000,
2852 2592000, 2621440, 2624400, 2654208, 2657205, 2700000, 2733750, 2764800, 2799360,
2853 2812500, 2834352, 2880000, 2916000, 2949120, 2952450, 2985984, 3000000, 3037500,
2854 3072000, 3110400, 3125000, 3145728, 3149280, 3188646, 3200000, 3240000, 3276800,
2855 3280500, 3317760, 3359232, 3375000, 3456000, 3499200, 3515625, 3538944, 3542940,
2856 3600000, 3645000, 3686400, 3732480, 3750000, 3779136, 3796875, 3840000, 3888000,
2857 3906250, 3932160, 3936600, 3981312, 4000000, 4050000, 4096000, 4100625, 4147200,
2858 4194304, 4199040, 4218750, 4251528, 4320000, 4374000, 4423680, 4428675, 4478976,
2859 4500000, 4556250, 4608000, 4665600, 4687500, 4718592, 4723920, 4782969, 4800000,
2860 4860000, 4915200, 4920750, 4976640, 5000000, 5038848, 5062500, 5120000, 5184000,
2861 5242880, 5248800, 5308416, 5314410, 5400000, 5467500, 5529600, 5598720, 5625000,
2862 5668704, 5760000, 5832000, 5859375, 5898240, 5904900, 5971968, 6000000, 6075000,
2863 6144000, 6220800, 6250000, 6291456, 6298560, 6328125, 6377292, 6400000, 6480000,
2864 6553600, 6561000, 6635520, 6718464, 6750000, 6834375, 6912000, 6998400, 7031250,
2865 7077888, 7085880, 7200000, 7290000, 7372800, 7381125, 7464960, 7500000, 7558272,
2866 7593750, 7680000, 7776000, 7812500, 7864320, 7873200, 7962624, 7971615, 8000000,
2867 8100000, 8192000, 8201250, 8294400, 8388608, 8398080, 8437500, 8503056, 8640000,
2868 8748000, 8847360, 8857350, 8957952, 9000000, 9112500, 9216000, 9331200, 9375000,
2869 9437184, 9447840, 9565938, 9600000, 9720000, 9765625, 9830400, 9841500, 9953280,
2870 10000000, 10077696, 10125000, 10240000, 10368000, 10485760, 10497600, 10546875, 10616832,
2871 10628820, 10800000, 10935000, 11059200, 11197440, 11250000, 11337408, 11390625, 11520000,
2872 11664000, 11718750, 11796480, 11809800, 11943936, 12000000, 12150000, 12288000, 12301875,
2873 12441600, 12500000, 12582912, 12597120, 12656250, 12754584, 12800000, 12960000, 13107200,
2874 13122000, 13271040, 13286025, 13436928, 13500000, 13668750, 13824000, 13996800, 14062500,
2875 14155776, 14171760, 14400000, 14580000, 14745600, 14762250, 14929920, 15000000, 15116544,
2876 15187500, 15360000, 15552000, 15625000, 15728640, 15746400, 15925248, 15943230, 16000000,
2877 16200000, 16384000, 16402500, 16588800, 16777216, 16796160, 16875000, 17006112, 17280000,
2878 17496000, 17578125, 17694720, 17714700, 17915904, 18000000, 18225000, 18432000, 18662400,
2879 18750000, 18874368, 18895680, 18984375, 19131876, 19200000, 19440000, 19531250, 19660800,
2880 19683000, 19906560, 20000000, 20155392, 20250000, 20480000, 20503125, 20736000, 20971520,
2881 20995200, 21093750, 21233664, 21257640, 21600000, 21870000, 22118400, 22143375, 22394880,
2882 22500000, 22674816, 22781250, 23040000, 23328000, 23437500, 23592960, 23619600, 23887872,
2883 23914845, 24000000, 24300000, 24576000, 24603750, 24883200, 25000000, 25165824, 25194240,
2884 25312500, 25509168, 25600000, 25920000, 26214400, 26244000, 26542080, 26572050, 26873856,
2885 27000000, 27337500, 27648000, 27993600, 28125000, 28311552, 28343520, 28800000, 29160000,
2886 29296875, 29491200, 29524500, 29859840, 30000000, 30233088, 30375000, 30720000, 31104000,
2887 31250000, 31457280, 31492800, 31640625, 31850496, 31886460, 32000000, 32400000, 32768000,
2888 32805000, 33177600, 33554432, 33592320, 33750000, 34012224, 34171875, 34560000, 34992000,
2889 35156250, 35389440, 35429400, 35831808, 36000000, 36450000, 36864000, 36905625, 37324800,
2890 37500000, 37748736, 37791360, 37968750, 38263752, 38400000, 38880000, 39062500, 39321600,
2891 39366000, 39813120, 39858075, 40000000, 40310784, 40500000, 40960000, 41006250, 41472000,
2892 41943040, 41990400, 42187500, 42467328, 42515280, 43200000, 43740000, 44236800, 44286750,
2893 44789760, 45000000, 45349632, 45562500, 46080000, 46656000, 46875000, 47185920, 47239200,
2894 47775744, 47829690, 48000000, 48600000, 48828125, 49152000, 49207500, 49766400, 50000000,
2895 50331648, 50388480, 50625000, 51018336, 51200000, 51840000, 52428800, 52488000, 52734375,
2896 53084160, 53144100, 53747712, 54000000, 54675000, 55296000, 55987200, 56250000, 56623104,
2897 56687040, 56953125, 57600000, 58320000, 58593750, 58982400, 59049000, 59719680, 60000000,
2898 60466176, 60750000, 61440000, 61509375, 62208000, 62500000, 62914560, 62985600, 63281250,
2899 63700992, 63772920, 64000000, 64800000, 65536000, 65610000, 66355200, 66430125, 67108864,
2900 67184640, 67500000, 68024448, 68343750, 69120000, 69984000, 70312500, 70778880, 70858800,
2901 71663616, 72000000, 72900000, 73728000, 73811250, 74649600, 75000000, 75497472, 75582720,
2902 75937500, 76527504, 76800000, 77760000, 78125000, 78643200, 78732000, 79626240, 79716150,
2903 80000000, 80621568, 81000000, 81920000, 82012500, 82944000, 83886080, 83980800, 84375000,
2904 84934656, 85030560, 86400000, 87480000, 87890625, 88473600, 88573500, 89579520, 90000000,
2905 90699264, 91125000, 92160000, 93312000, 93750000, 94371840, 94478400, 94921875, 95551488,
2906 95659380, 96000000, 97200000, 97656250, 98304000, 98415000, 99532800, 100000000,
2907 100663296, 100776960, 101250000, 102036672, 102400000, 102515625, 103680000, 104857600,
2908 104976000, 105468750, 106168320, 106288200, 107495424, 108000000, 109350000, 110592000,
2909 110716875, 111974400, 112500000, 113246208, 113374080, 113906250, 115200000, 116640000,
2910 117187500, 117964800, 118098000, 119439360, 119574225, 120000000, 120932352, 121500000,
2911 122880000, 123018750, 124416000, 125000000, 125829120, 125971200, 126562500, 127401984,
2912 127545840, 128000000, 129600000, 131072000, 131220000, 132710400, 132860250, 134217728,
2913 134369280, 135000000, 136048896, 136687500, 138240000, 139968000, 140625000, 141557760,
2914 141717600, 143327232, 144000000, 145800000, 146484375, 147456000, 147622500, 149299200,
2915 150000000, 150994944, 151165440, 151875000, 153055008, 153600000, 155520000, 156250000,
2916 157286400, 157464000, 158203125, 159252480, 159432300, 160000000, 161243136, 162000000,
2917 163840000, 164025000, 165888000, 167772160, 167961600, 168750000, 169869312, 170061120,
2918 170859375, 172800000, 174960000, 175781250, 176947200, 177147000, 179159040, 180000000,
2919 181398528, 182250000, 184320000, 184528125, 186624000, 187500000, 188743680, 188956800,
2920 189843750, 191102976, 191318760, 192000000, 194400000, 195312500, 196608000, 196830000,
2921 199065600, 199290375, 200000000, 201326592, 201553920, 202500000, 204073344, 204800000,
2922 205031250, 207360000, 209715200, 209952000, 210937500, 212336640, 212576400, 214990848,
2923 216000000, 218700000, 221184000, 221433750, 223948800, 225000000, 226492416, 226748160,
2924 227812500, 230400000, 233280000, 234375000, 235929600, 236196000, 238878720, 239148450,
2925 240000000, 241864704, 243000000, 244140625, 245760000, 246037500, 248832000, 250000000,
2926 251658240, 251942400, 253125000, 254803968, 255091680, 256000000, 259200000, 262144000,
2927 262440000, 263671875, 265420800, 265720500, 268435456, 268738560, 270000000, 272097792,
2928 273375000, 276480000, 279936000, 281250000, 283115520, 283435200, 284765625, 286654464,
2929 288000000, 291600000, 292968750, 294912000, 295245000, 298598400, 300000000, 301989888,
2930 302330880, 303750000, 306110016, 307200000, 307546875, 311040000, 312500000, 314572800,
2931 314928000, 316406250, 318504960, 318864600, 320000000, 322486272, 324000000, 327680000,
2932 328050000, 331776000, 332150625, 335544320, 335923200, 337500000, 339738624, 340122240,
2933 341718750, 345600000, 349920000, 351562500, 353894400, 354294000, 358318080, 360000000,
2934 362797056, 364500000, 368640000, 369056250, 373248000, 375000000, 377487360, 377913600,
2935 379687500, 382205952, 382637520, 384000000, 388800000, 390625000, 393216000, 393660000,
2936 398131200, 398580750, 400000000, 402653184, 403107840, 405000000, 408146688, 409600000,
2937 410062500, 414720000, 419430400, 419904000, 421875000, 424673280, 425152800, 429981696,
2938 432000000, 437400000, 439453125, 442368000, 442867500, 447897600, 450000000, 452984832,
2939 453496320, 455625000, 460800000, 466560000, 468750000, 471859200, 472392000, 474609375,
2940 477757440, 478296900, 480000000, 483729408, 486000000, 488281250, 491520000, 492075000,
2941 497664000, 500000000, 503316480, 503884800, 506250000, 509607936, 510183360, 512000000,
2942 512578125, 518400000, 524288000, 524880000, 527343750, 530841600, 531441000, 536870912,
2943 537477120, 540000000, 544195584, 546750000, 552960000, 553584375, 559872000, 562500000,
2944 566231040, 566870400, 569531250, 573308928, 576000000, 583200000, 585937500, 589824000,
2945 590490000, 597196800, 597871125, 600000000, 603979776, 604661760, 607500000, 612220032,
2946 614400000, 615093750, 622080000, 625000000, 629145600, 629856000, 632812500, 637009920,
2947 637729200, 640000000, 644972544, 648000000, 655360000, 656100000, 663552000, 664301250,
2948 671088640, 671846400, 675000000, 679477248, 680244480, 683437500, 691200000, 699840000,
2949 703125000, 707788800, 708588000, 716636160, 720000000, 725594112, 729000000, 732421875,
2950 737280000, 738112500, 746496000, 750000000, 754974720, 755827200, 759375000, 764411904,
2951 765275040, 768000000, 777600000, 781250000, 786432000, 787320000, 791015625, 796262400,
2952 797161500, 800000000, 805306368, 806215680, 810000000, 816293376, 819200000, 820125000,
2953 829440000, 838860800, 839808000, 843750000, 849346560, 850305600, 854296875, 859963392,
2954 864000000, 874800000, 878906250, 884736000, 885735000, 895795200, 900000000, 905969664,
2955 906992640, 911250000, 921600000, 922640625, 933120000, 937500000, 943718400, 944784000,
2956 949218750, 955514880, 956593800, 960000000, 967458816, 972000000, 976562500, 983040000,
2957 984150000, 995328000, 996451875, 1000000000, 1006632960, 1007769600, 1012500000,
2958 1019215872, 1020366720, 1024000000, 1025156250, 1036800000, 1048576000, 1049760000,
2959 1054687500, 1061683200, 1062882000, 1073741824, 1074954240, 1080000000, 1088391168,
2960 1093500000, 1105920000, 1107168750, 1119744000, 1125000000, 1132462080, 1133740800,
2961 1139062500, 1146617856, 1152000000, 1166400000, 1171875000, 1179648000, 1180980000,
2962 1194393600, 1195742250, 1200000000, 1207959552, 1209323520, 1215000000, 1220703125,
2963 1224440064, 1228800000, 1230187500, 1244160000, 1250000000, 1258291200, 1259712000,
2964 1265625000, 1274019840, 1275458400, 1280000000, 1289945088, 1296000000, 1310720000,
2965 1312200000, 1318359375, 1327104000, 1328602500, 1342177280, 1343692800, 1350000000,
2966 1358954496, 1360488960, 1366875000, 1382400000, 1399680000, 1406250000, 1415577600,
2967 1417176000, 1423828125, 1433272320, 1440000000, 1451188224, 1458000000, 1464843750,
2968 1474560000, 1476225000, 1492992000, 1500000000, 1509949440, 1511654400, 1518750000,
2969 1528823808, 1530550080, 1536000000, 1537734375, 1555200000, 1562500000, 1572864000,
2970 1574640000, 1582031250, 1592524800, 1594323000, 1600000000, 1610612736, 1612431360,
2971 1620000000, 1632586752, 1638400000, 1640250000, 1658880000, 1660753125, 1677721600,
2972 1679616000, 1687500000, 1698693120, 1700611200, 1708593750, 1719926784, 1728000000,
2973 1749600000, 1757812500, 1769472000, 1771470000, 1791590400, 1800000000, 1811939328,
2974 1813985280, 1822500000, 1843200000, 1845281250, 1866240000, 1875000000, 1887436800,
2975 1889568000, 1898437500, 1911029760, 1913187600, 1920000000, 1934917632, 1944000000,
2976 1953125000, 1966080000, 1968300000, 1990656000, 1992903750, 2000000000, 2013265920,
2977 2015539200, 2025000000, 2038431744, 2040733440, 2048000000, 2050312500, 2073600000,
2978 2097152000, 2099520000, 2109375000, 2123366400, 2125764000
2983 cvGetOptimalDFTSize( int size0 )
2985 int a = 0, b = sizeof(icvOptimalDFTSize)/sizeof(icvOptimalDFTSize[0]) - 1;
2986 if( (unsigned)size0 >= (unsigned)icvOptimalDFTSize[b] )
2991 int c = (a + b) >> 1;
2992 if( size0 <= icvOptimalDFTSize[c] )
2998 return icvOptimalDFTSize[b];