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
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;
54
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;
60
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;
66
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;
72
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;
77
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;
82
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;
87
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;*/
92
93 /****************************************************************************************\
94 Discrete Fourier Transform
95 \****************************************************************************************/
96
97 #define CV_MAX_LOCAL_DFT_SIZE (1 << 15)
98
99 static const uchar log2tab[] = { 0, 0, 1, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0 };
icvlog2(int n)100 static 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
114 static 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
134 static 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
176 static int
icvDFTFactorize(int n,int * factors)177 icvDFTFactorize( 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
220 static void
icvDFTInit(int n0,int nf,int * factors,int * itab,int elem_size,void * _wave,int inv_itab)221 icvDFTInit( 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
421 static const double icv_sin_120 = 0.86602540378443864676372317075294;
422 static const double icv_sin_45 = 0.70710678118654752440084436210485;
423 static const double icv_fft5_2 = 0.559016994374947424102293417182819;
424 static const double icv_fft5_3 = -0.951056516295153572116439333379382;
425 static const double icv_fft5_4 = -1.538841768587626701285145288018455;
426 static 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
432 static CvStatus CV_STDCALL
icvDFT_64fc(const CvComplex64f * src,CvComplex64f * dst,int n,int nf,int * factors,const int * itab,const CvComplex64f * wave,int tab_size,const void * spec,CvComplex64f * buf,int flags,double scale)433 icvDFT_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
855 static CvStatus CV_STDCALL
icvDFT_32fc(const CvComplex32f * src,CvComplex32f * dst,int n,int nf,int * factors,const int * itab,const CvComplex32f * wave,int tab_size,const void * spec,CvComplex32f * buf,int flags,double scale)856 icvDFT_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 ) \
1288 static CvStatus CV_STDCALL \
1289 icvRealDFT_##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 \
1389 finalize: \
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 ) \
1407 static CvStatus CV_STDCALL \
1408 icvCCSIDFT_##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 \
1557 finalize: \
1558 if( complex_input ) \
1559 ((datatype*)src)[0] = (datatype)save_s1; \
1560 \
1561 return CV_OK; \
1562 }
1563
1564
1565 ICV_REAL_DFT( 64f, double )
1566 ICV_CCS_IDFT( 64f, double )
1567 ICV_REAL_DFT( 32f, float )
1568 ICV_CCS_IDFT( 32f, float )
1569
1570
1571 static void
icvCopyColumn(const uchar * _src,int src_step,uchar * _dst,int dst_step,int len,int elem_size)1572 icvCopyColumn( 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
1608 static void
icvCopyFrom2Columns(const uchar * _src,int src_step,uchar * _dst0,uchar * _dst1,int len,int elem_size)1609 icvCopyFrom2Columns( 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
1654 static void
icvCopyTo2Columns(const uchar * _src0,const uchar * _src1,uchar * _dst,int dst_step,int len,int elem_size)1655 icvCopyTo2Columns( 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
1700 static void
icvExpandCCS(uchar * _ptr,int len,int elem_size)1701 icvExpandCCS( 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
1737 typedef 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
1742 CV_IMPL void
cvDFT(const CvArr * srcarr,CvArr * dstarr,int flags,int nonzero_rows)1743 cvDFT( 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
2206 CV_IMPL void
cvMulSpectrums(const CvArr * srcAarr,const CvArr * srcBarr,CvArr * dstarr,int flags)2207 cvMulSpectrums( 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 ) \
2401 static CvStatus CV_STDCALL \
2402 icvDCT_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 ) \
2448 static CvStatus CV_STDCALL \
2449 icvDCT_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
2492 ICV_DCT_FWD( 64f, double )
2493 ICV_DCT_INV( 64f, double )
2494 ICV_DCT_FWD( 32f, float )
2495 ICV_DCT_INV( 32f, float )
2496
2497 static void
icvDCTInit(int n,int elem_size,void * _wave,int inv)2498 icvDCTInit( 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
2574 typedef 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
2580 CV_IMPL void
cvDCT(const CvArr * srcarr,CvArr * dstarr,int flags)2581 cvDCT( 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
2802 static const int icvOptimalDFTSize[] = {
2803 1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48,
2804 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, 128, 135, 144, 150, 160,
2805 162, 180, 192, 200, 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375,
2806 384, 400, 405, 432, 450, 480, 486, 500, 512, 540, 576, 600, 625, 640, 648, 675, 720,
2807 729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 1080, 1125, 1152, 1200,
2808 1215, 1250, 1280, 1296, 1350, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875,
2809 1920, 1944, 2000, 2025, 2048, 2160, 2187, 2250, 2304, 2400, 2430, 2500, 2560, 2592,
2810 2700, 2880, 2916, 3000, 3072, 3125, 3200, 3240, 3375, 3456, 3600, 3645, 3750, 3840,
2811 3888, 4000, 4050, 4096, 4320, 4374, 4500, 4608, 4800, 4860, 5000, 5120, 5184, 5400,
2812 5625, 5760, 5832, 6000, 6075, 6144, 6250, 6400, 6480, 6561, 6750, 6912, 7200, 7290,
2813 7500, 7680, 7776, 8000, 8100, 8192, 8640, 8748, 9000, 9216, 9375, 9600, 9720, 10000,
2814 10125, 10240, 10368, 10800, 10935, 11250, 11520, 11664, 12000, 12150, 12288, 12500,
2815 12800, 12960, 13122, 13500, 13824, 14400, 14580, 15000, 15360, 15552, 15625, 16000,
2816 16200, 16384, 16875, 17280, 17496, 18000, 18225, 18432, 18750, 19200, 19440, 19683,
2817 20000, 20250, 20480, 20736, 21600, 21870, 22500, 23040, 23328, 24000, 24300, 24576,
2818 25000, 25600, 25920, 26244, 27000, 27648, 28125, 28800, 29160, 30000, 30375, 30720,
2819 31104, 31250, 32000, 32400, 32768, 32805, 33750, 34560, 34992, 36000, 36450, 36864,
2820 37500, 38400, 38880, 39366, 40000, 40500, 40960, 41472, 43200, 43740, 45000, 46080,
2821 46656, 46875, 48000, 48600, 49152, 50000, 50625, 51200, 51840, 52488, 54000, 54675,
2822 55296, 56250, 57600, 58320, 59049, 60000, 60750, 61440, 62208, 62500, 64000, 64800,
2823 65536, 65610, 67500, 69120, 69984, 72000, 72900, 73728, 75000, 76800, 77760, 78125,
2824 78732, 80000, 81000, 81920, 82944, 84375, 86400, 87480, 90000, 91125, 92160, 93312,
2825 93750, 96000, 97200, 98304, 98415, 100000, 101250, 102400, 103680, 104976, 108000,
2826 109350, 110592, 112500, 115200, 116640, 118098, 120000, 121500, 122880, 124416, 125000,
2827 128000, 129600, 131072, 131220, 135000, 138240, 139968, 140625, 144000, 145800, 147456,
2828 150000, 151875, 153600, 155520, 156250, 157464, 160000, 162000, 163840, 164025, 165888,
2829 168750, 172800, 174960, 177147, 180000, 182250, 184320, 186624, 187500, 192000, 194400,
2830 196608, 196830, 200000, 202500, 204800, 207360, 209952, 216000, 218700, 221184, 225000,
2831 230400, 233280, 234375, 236196, 240000, 243000, 245760, 248832, 250000, 253125, 256000,
2832 259200, 262144, 262440, 270000, 273375, 276480, 279936, 281250, 288000, 291600, 294912,
2833 295245, 300000, 303750, 307200, 311040, 312500, 314928, 320000, 324000, 327680, 328050,
2834 331776, 337500, 345600, 349920, 354294, 360000, 364500, 368640, 373248, 375000, 384000,
2835 388800, 390625, 393216, 393660, 400000, 405000, 409600, 414720, 419904, 421875, 432000,
2836 437400, 442368, 450000, 455625, 460800, 466560, 468750, 472392, 480000, 486000, 491520,
2837 492075, 497664, 500000, 506250, 512000, 518400, 524288, 524880, 531441, 540000, 546750,
2838 552960, 559872, 562500, 576000, 583200, 589824, 590490, 600000, 607500, 614400, 622080,
2839 625000, 629856, 640000, 648000, 655360, 656100, 663552, 675000, 691200, 699840, 703125,
2840 708588, 720000, 729000, 737280, 746496, 750000, 759375, 768000, 777600, 781250, 786432,
2841 787320, 800000, 810000, 819200, 820125, 829440, 839808, 843750, 864000, 874800, 884736,
2842 885735, 900000, 911250, 921600, 933120, 937500, 944784, 960000, 972000, 983040, 984150,
2843 995328, 1000000, 1012500, 1024000, 1036800, 1048576, 1049760, 1062882, 1080000, 1093500,
2844 1105920, 1119744, 1125000, 1152000, 1166400, 1171875, 1179648, 1180980, 1200000,
2845 1215000, 1228800, 1244160, 1250000, 1259712, 1265625, 1280000, 1296000, 1310720,
2846 1312200, 1327104, 1350000, 1366875, 1382400, 1399680, 1406250, 1417176, 1440000,
2847 1458000, 1474560, 1476225, 1492992, 1500000, 1518750, 1536000, 1555200, 1562500,
2848 1572864, 1574640, 1594323, 1600000, 1620000, 1638400, 1640250, 1658880, 1679616,
2849 1687500, 1728000, 1749600, 1769472, 1771470, 1800000, 1822500, 1843200, 1866240,
2850 1875000, 1889568, 1920000, 1944000, 1953125, 1966080, 1968300, 1990656, 2000000,
2851 2025000, 2048000, 2073600, 2097152, 2099520, 2109375, 2125764, 2160000, 2187000,
2852 2211840, 2239488, 2250000, 2278125, 2304000, 2332800, 2343750, 2359296, 2361960,
2853 2400000, 2430000, 2457600, 2460375, 2488320, 2500000, 2519424, 2531250, 2560000,
2854 2592000, 2621440, 2624400, 2654208, 2657205, 2700000, 2733750, 2764800, 2799360,
2855 2812500, 2834352, 2880000, 2916000, 2949120, 2952450, 2985984, 3000000, 3037500,
2856 3072000, 3110400, 3125000, 3145728, 3149280, 3188646, 3200000, 3240000, 3276800,
2857 3280500, 3317760, 3359232, 3375000, 3456000, 3499200, 3515625, 3538944, 3542940,
2858 3600000, 3645000, 3686400, 3732480, 3750000, 3779136, 3796875, 3840000, 3888000,
2859 3906250, 3932160, 3936600, 3981312, 4000000, 4050000, 4096000, 4100625, 4147200,
2860 4194304, 4199040, 4218750, 4251528, 4320000, 4374000, 4423680, 4428675, 4478976,
2861 4500000, 4556250, 4608000, 4665600, 4687500, 4718592, 4723920, 4782969, 4800000,
2862 4860000, 4915200, 4920750, 4976640, 5000000, 5038848, 5062500, 5120000, 5184000,
2863 5242880, 5248800, 5308416, 5314410, 5400000, 5467500, 5529600, 5598720, 5625000,
2864 5668704, 5760000, 5832000, 5859375, 5898240, 5904900, 5971968, 6000000, 6075000,
2865 6144000, 6220800, 6250000, 6291456, 6298560, 6328125, 6377292, 6400000, 6480000,
2866 6553600, 6561000, 6635520, 6718464, 6750000, 6834375, 6912000, 6998400, 7031250,
2867 7077888, 7085880, 7200000, 7290000, 7372800, 7381125, 7464960, 7500000, 7558272,
2868 7593750, 7680000, 7776000, 7812500, 7864320, 7873200, 7962624, 7971615, 8000000,
2869 8100000, 8192000, 8201250, 8294400, 8388608, 8398080, 8437500, 8503056, 8640000,
2870 8748000, 8847360, 8857350, 8957952, 9000000, 9112500, 9216000, 9331200, 9375000,
2871 9437184, 9447840, 9565938, 9600000, 9720000, 9765625, 9830400, 9841500, 9953280,
2872 10000000, 10077696, 10125000, 10240000, 10368000, 10485760, 10497600, 10546875, 10616832,
2873 10628820, 10800000, 10935000, 11059200, 11197440, 11250000, 11337408, 11390625, 11520000,
2874 11664000, 11718750, 11796480, 11809800, 11943936, 12000000, 12150000, 12288000, 12301875,
2875 12441600, 12500000, 12582912, 12597120, 12656250, 12754584, 12800000, 12960000, 13107200,
2876 13122000, 13271040, 13286025, 13436928, 13500000, 13668750, 13824000, 13996800, 14062500,
2877 14155776, 14171760, 14400000, 14580000, 14745600, 14762250, 14929920, 15000000, 15116544,
2878 15187500, 15360000, 15552000, 15625000, 15728640, 15746400, 15925248, 15943230, 16000000,
2879 16200000, 16384000, 16402500, 16588800, 16777216, 16796160, 16875000, 17006112, 17280000,
2880 17496000, 17578125, 17694720, 17714700, 17915904, 18000000, 18225000, 18432000, 18662400,
2881 18750000, 18874368, 18895680, 18984375, 19131876, 19200000, 19440000, 19531250, 19660800,
2882 19683000, 19906560, 20000000, 20155392, 20250000, 20480000, 20503125, 20736000, 20971520,
2883 20995200, 21093750, 21233664, 21257640, 21600000, 21870000, 22118400, 22143375, 22394880,
2884 22500000, 22674816, 22781250, 23040000, 23328000, 23437500, 23592960, 23619600, 23887872,
2885 23914845, 24000000, 24300000, 24576000, 24603750, 24883200, 25000000, 25165824, 25194240,
2886 25312500, 25509168, 25600000, 25920000, 26214400, 26244000, 26542080, 26572050, 26873856,
2887 27000000, 27337500, 27648000, 27993600, 28125000, 28311552, 28343520, 28800000, 29160000,
2888 29296875, 29491200, 29524500, 29859840, 30000000, 30233088, 30375000, 30720000, 31104000,
2889 31250000, 31457280, 31492800, 31640625, 31850496, 31886460, 32000000, 32400000, 32768000,
2890 32805000, 33177600, 33554432, 33592320, 33750000, 34012224, 34171875, 34560000, 34992000,
2891 35156250, 35389440, 35429400, 35831808, 36000000, 36450000, 36864000, 36905625, 37324800,
2892 37500000, 37748736, 37791360, 37968750, 38263752, 38400000, 38880000, 39062500, 39321600,
2893 39366000, 39813120, 39858075, 40000000, 40310784, 40500000, 40960000, 41006250, 41472000,
2894 41943040, 41990400, 42187500, 42467328, 42515280, 43200000, 43740000, 44236800, 44286750,
2895 44789760, 45000000, 45349632, 45562500, 46080000, 46656000, 46875000, 47185920, 47239200,
2896 47775744, 47829690, 48000000, 48600000, 48828125, 49152000, 49207500, 49766400, 50000000,
2897 50331648, 50388480, 50625000, 51018336, 51200000, 51840000, 52428800, 52488000, 52734375,
2898 53084160, 53144100, 53747712, 54000000, 54675000, 55296000, 55987200, 56250000, 56623104,
2899 56687040, 56953125, 57600000, 58320000, 58593750, 58982400, 59049000, 59719680, 60000000,
2900 60466176, 60750000, 61440000, 61509375, 62208000, 62500000, 62914560, 62985600, 63281250,
2901 63700992, 63772920, 64000000, 64800000, 65536000, 65610000, 66355200, 66430125, 67108864,
2902 67184640, 67500000, 68024448, 68343750, 69120000, 69984000, 70312500, 70778880, 70858800,
2903 71663616, 72000000, 72900000, 73728000, 73811250, 74649600, 75000000, 75497472, 75582720,
2904 75937500, 76527504, 76800000, 77760000, 78125000, 78643200, 78732000, 79626240, 79716150,
2905 80000000, 80621568, 81000000, 81920000, 82012500, 82944000, 83886080, 83980800, 84375000,
2906 84934656, 85030560, 86400000, 87480000, 87890625, 88473600, 88573500, 89579520, 90000000,
2907 90699264, 91125000, 92160000, 93312000, 93750000, 94371840, 94478400, 94921875, 95551488,
2908 95659380, 96000000, 97200000, 97656250, 98304000, 98415000, 99532800, 100000000,
2909 100663296, 100776960, 101250000, 102036672, 102400000, 102515625, 103680000, 104857600,
2910 104976000, 105468750, 106168320, 106288200, 107495424, 108000000, 109350000, 110592000,
2911 110716875, 111974400, 112500000, 113246208, 113374080, 113906250, 115200000, 116640000,
2912 117187500, 117964800, 118098000, 119439360, 119574225, 120000000, 120932352, 121500000,
2913 122880000, 123018750, 124416000, 125000000, 125829120, 125971200, 126562500, 127401984,
2914 127545840, 128000000, 129600000, 131072000, 131220000, 132710400, 132860250, 134217728,
2915 134369280, 135000000, 136048896, 136687500, 138240000, 139968000, 140625000, 141557760,
2916 141717600, 143327232, 144000000, 145800000, 146484375, 147456000, 147622500, 149299200,
2917 150000000, 150994944, 151165440, 151875000, 153055008, 153600000, 155520000, 156250000,
2918 157286400, 157464000, 158203125, 159252480, 159432300, 160000000, 161243136, 162000000,
2919 163840000, 164025000, 165888000, 167772160, 167961600, 168750000, 169869312, 170061120,
2920 170859375, 172800000, 174960000, 175781250, 176947200, 177147000, 179159040, 180000000,
2921 181398528, 182250000, 184320000, 184528125, 186624000, 187500000, 188743680, 188956800,
2922 189843750, 191102976, 191318760, 192000000, 194400000, 195312500, 196608000, 196830000,
2923 199065600, 199290375, 200000000, 201326592, 201553920, 202500000, 204073344, 204800000,
2924 205031250, 207360000, 209715200, 209952000, 210937500, 212336640, 212576400, 214990848,
2925 216000000, 218700000, 221184000, 221433750, 223948800, 225000000, 226492416, 226748160,
2926 227812500, 230400000, 233280000, 234375000, 235929600, 236196000, 238878720, 239148450,
2927 240000000, 241864704, 243000000, 244140625, 245760000, 246037500, 248832000, 250000000,
2928 251658240, 251942400, 253125000, 254803968, 255091680, 256000000, 259200000, 262144000,
2929 262440000, 263671875, 265420800, 265720500, 268435456, 268738560, 270000000, 272097792,
2930 273375000, 276480000, 279936000, 281250000, 283115520, 283435200, 284765625, 286654464,
2931 288000000, 291600000, 292968750, 294912000, 295245000, 298598400, 300000000, 301989888,
2932 302330880, 303750000, 306110016, 307200000, 307546875, 311040000, 312500000, 314572800,
2933 314928000, 316406250, 318504960, 318864600, 320000000, 322486272, 324000000, 327680000,
2934 328050000, 331776000, 332150625, 335544320, 335923200, 337500000, 339738624, 340122240,
2935 341718750, 345600000, 349920000, 351562500, 353894400, 354294000, 358318080, 360000000,
2936 362797056, 364500000, 368640000, 369056250, 373248000, 375000000, 377487360, 377913600,
2937 379687500, 382205952, 382637520, 384000000, 388800000, 390625000, 393216000, 393660000,
2938 398131200, 398580750, 400000000, 402653184, 403107840, 405000000, 408146688, 409600000,
2939 410062500, 414720000, 419430400, 419904000, 421875000, 424673280, 425152800, 429981696,
2940 432000000, 437400000, 439453125, 442368000, 442867500, 447897600, 450000000, 452984832,
2941 453496320, 455625000, 460800000, 466560000, 468750000, 471859200, 472392000, 474609375,
2942 477757440, 478296900, 480000000, 483729408, 486000000, 488281250, 491520000, 492075000,
2943 497664000, 500000000, 503316480, 503884800, 506250000, 509607936, 510183360, 512000000,
2944 512578125, 518400000, 524288000, 524880000, 527343750, 530841600, 531441000, 536870912,
2945 537477120, 540000000, 544195584, 546750000, 552960000, 553584375, 559872000, 562500000,
2946 566231040, 566870400, 569531250, 573308928, 576000000, 583200000, 585937500, 589824000,
2947 590490000, 597196800, 597871125, 600000000, 603979776, 604661760, 607500000, 612220032,
2948 614400000, 615093750, 622080000, 625000000, 629145600, 629856000, 632812500, 637009920,
2949 637729200, 640000000, 644972544, 648000000, 655360000, 656100000, 663552000, 664301250,
2950 671088640, 671846400, 675000000, 679477248, 680244480, 683437500, 691200000, 699840000,
2951 703125000, 707788800, 708588000, 716636160, 720000000, 725594112, 729000000, 732421875,
2952 737280000, 738112500, 746496000, 750000000, 754974720, 755827200, 759375000, 764411904,
2953 765275040, 768000000, 777600000, 781250000, 786432000, 787320000, 791015625, 796262400,
2954 797161500, 800000000, 805306368, 806215680, 810000000, 816293376, 819200000, 820125000,
2955 829440000, 838860800, 839808000, 843750000, 849346560, 850305600, 854296875, 859963392,
2956 864000000, 874800000, 878906250, 884736000, 885735000, 895795200, 900000000, 905969664,
2957 906992640, 911250000, 921600000, 922640625, 933120000, 937500000, 943718400, 944784000,
2958 949218750, 955514880, 956593800, 960000000, 967458816, 972000000, 976562500, 983040000,
2959 984150000, 995328000, 996451875, 1000000000, 1006632960, 1007769600, 1012500000,
2960 1019215872, 1020366720, 1024000000, 1025156250, 1036800000, 1048576000, 1049760000,
2961 1054687500, 1061683200, 1062882000, 1073741824, 1074954240, 1080000000, 1088391168,
2962 1093500000, 1105920000, 1107168750, 1119744000, 1125000000, 1132462080, 1133740800,
2963 1139062500, 1146617856, 1152000000, 1166400000, 1171875000, 1179648000, 1180980000,
2964 1194393600, 1195742250, 1200000000, 1207959552, 1209323520, 1215000000, 1220703125,
2965 1224440064, 1228800000, 1230187500, 1244160000, 1250000000, 1258291200, 1259712000,
2966 1265625000, 1274019840, 1275458400, 1280000000, 1289945088, 1296000000, 1310720000,
2967 1312200000, 1318359375, 1327104000, 1328602500, 1342177280, 1343692800, 1350000000,
2968 1358954496, 1360488960, 1366875000, 1382400000, 1399680000, 1406250000, 1415577600,
2969 1417176000, 1423828125, 1433272320, 1440000000, 1451188224, 1458000000, 1464843750,
2970 1474560000, 1476225000, 1492992000, 1500000000, 1509949440, 1511654400, 1518750000,
2971 1528823808, 1530550080, 1536000000, 1537734375, 1555200000, 1562500000, 1572864000,
2972 1574640000, 1582031250, 1592524800, 1594323000, 1600000000, 1610612736, 1612431360,
2973 1620000000, 1632586752, 1638400000, 1640250000, 1658880000, 1660753125, 1677721600,
2974 1679616000, 1687500000, 1698693120, 1700611200, 1708593750, 1719926784, 1728000000,
2975 1749600000, 1757812500, 1769472000, 1771470000, 1791590400, 1800000000, 1811939328,
2976 1813985280, 1822500000, 1843200000, 1845281250, 1866240000, 1875000000, 1887436800,
2977 1889568000, 1898437500, 1911029760, 1913187600, 1920000000, 1934917632, 1944000000,
2978 1953125000, 1966080000, 1968300000, 1990656000, 1992903750, 2000000000, 2013265920,
2979 2015539200, 2025000000, 2038431744, 2040733440, 2048000000, 2050312500, 2073600000,
2980 2097152000, 2099520000, 2109375000, 2123366400, 2125764000
2981 };
2982
2983
2984 CV_IMPL int
cvGetOptimalDFTSize(int size0)2985 cvGetOptimalDFTSize( 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