vcs.maemo.org Git - opencv/blob - src/cxcore/cxdxt.cpp

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