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