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
1462 void cv::dft( InputArray _src0, OutputArray _dst, int flags, int nonzero_rows )
1463 {
1464     static DFTFunc dft_tbl[6] =
1465     {
1466         (DFTFunc)DFT_32f,
1467         (DFTFunc)RealDFT_32f,
1468         (DFTFunc)CCSIDFT_32f,
1469         (DFTFunc)DFT_64f,
1470         (DFTFunc)RealDFT_64f,
1471         (DFTFunc)CCSIDFT_64f
1472     };
1473
1474     AutoBuffer<uchar> buf;
1475     void *spec = 0;
1476
1477     Mat src0 = _src0.getMat(), src = src0;
1478     int prev_len = 0, stage = 0;
1479     bool inv = (flags & DFT_INVERSE) != 0;
1480     int nf = 0, real_transform = src.channels() == 1 || (inv && (flags & DFT_REAL_OUTPUT)!=0);
1481     int type = src.type(), depth = src.depth();
1482     int elem_size = (int)src.elemSize1(), complex_elem_size = elem_size*2;
1483     int factors[34];
1484     bool inplace_transform = false;
1485 #ifdef HAVE_IPP
1486     void *spec_r = 0, *spec_c = 0;
1487     int ipp_norm_flag = !(flags & DFT_SCALE) ? 8 : inv ? 2 : 1;
1488 #endif
1489
1490     CV_Assert( type == CV_32FC1 || type == CV_32FC2 || type == CV_64FC1 || type == CV_64FC2 );
1491
1492     if( !inv && src.channels() == 1 && (flags & DFT_COMPLEX_OUTPUT) )
1493         _dst.create( src.size(), CV_MAKETYPE(depth, 2) );
1494     else if( inv && src.channels() == 2 && (flags & DFT_REAL_OUTPUT) )
1495         _dst.create( src.size(), depth );
1496     else
1497         _dst.create( src.size(), type );
1498
1499     Mat dst = _dst.getMat();
1500
1501     if( !real_transform )
1502         elem_size = complex_elem_size;
1503
1504     if( src.cols == 1 && nonzero_rows > 0 )
1505         CV_Error( CV_StsNotImplemented,
1506         "This mode (using nonzero_rows with a single-column matrix) breaks the function's logic, so it is prohibited.\n"
1507         "For fast convolution/correlation use 2-column matrix or single-row matrix instead" );
1508
1509     // determine, which transform to do first - row-wise
1510     // (stage 0) or column-wise (stage 1) transform
1511     if( !(flags & DFT_ROWS) && src.rows > 1 &&
1512         ((src.cols == 1 && (!src.isContinuous() || !dst.isContinuous())) ||
1513          (src.cols > 1 && inv && real_transform)) )
1514         stage = 1;
1515
1516     for(;;)
1517     {
1518         double scale = 1;
1519         uchar* wave = 0;
1520         int* itab = 0;
1521         uchar* ptr;
1522         int i, len, count, sz = 0;
1523         int use_buf = 0, odd_real = 0;
1524         DFTFunc dft_func;
1525
1526         if( stage == 0 ) // row-wise transform
1527         {
1528             len = !inv ? src.cols : dst.cols;
1529             count = src.rows;
1530             if( len == 1 && !(flags & DFT_ROWS) )
1531             {
1532                 len = !inv ? src.rows : dst.rows;
1533                 count = 1;
1534             }
1535             odd_real = real_transform && (len & 1);
1536         }
1537         else
1538         {
1539             len = dst.rows;
1540             count = !inv ? src0.cols : dst.cols;
1541             sz = 2*len*complex_elem_size;
1542         }
1543
1544         spec = 0;
1545 #ifdef HAVE_IPP
1546         if( len*count >= 64 ) // use IPP DFT if available
1547         {
1548             int ipp_sz = 0;
1549
1550             if( real_transform && stage == 0 )
1551             {
1552                 if( depth == CV_32F )
1553                 {
1554                     if( spec_r )
1555                         IPPI_CALL( ippsDFTFree_R_32f( (IppsDFTSpec_R_32f*)spec_r ));
1556                     IPPI_CALL( ippsDFTInitAlloc_R_32f(
1557                         (IppsDFTSpec_R_32f**)&spec_r, len, ipp_norm_flag, ippAlgHintNone ));
1558                     IPPI_CALL( ippsDFTGetBufSize_R_32f( (IppsDFTSpec_R_32f*)spec_r, &ipp_sz ));
1559                 }
1560                 else
1561                 {
1562                     if( spec_r )
1563                         IPPI_CALL( ippsDFTFree_R_64f( (IppsDFTSpec_R_64f*)spec_r ));
1564                     IPPI_CALL( ippsDFTInitAlloc_R_64f(
1565                         (IppsDFTSpec_R_64f**)&spec_r, len, ipp_norm_flag, ippAlgHintNone ));
1566                     IPPI_CALL( ippsDFTGetBufSize_R_64f( (IppsDFTSpec_R_64f*)spec_r, &ipp_sz ));
1567                 }
1568                 spec = spec_r;
1569             }
1570             else
1571             {
1572                 if( depth == CV_32F )
1573                 {
1574                     if( spec_c )
1575                         IPPI_CALL( ippsDFTFree_C_32fc( (IppsDFTSpec_C_32fc*)spec_c ));
1576                     IPPI_CALL( ippsDFTInitAlloc_C_32fc(
1577                         (IppsDFTSpec_C_32fc**)&spec_c, len, ipp_norm_flag, ippAlgHintNone ));
1578                     IPPI_CALL( ippsDFTGetBufSize_C_32fc( (IppsDFTSpec_C_32fc*)spec_c, &ipp_sz ));
1579                 }
1580                 else
1581                 {
1582                     if( spec_c )
1583                         IPPI_CALL( ippsDFTFree_C_64fc( (IppsDFTSpec_C_64fc*)spec_c ));
1584                     IPPI_CALL( ippsDFTInitAlloc_C_64fc(
1585                         (IppsDFTSpec_C_64fc**)&spec_c, len, ipp_norm_flag, ippAlgHintNone ));
1586                     IPPI_CALL( ippsDFTGetBufSize_C_64fc( (IppsDFTSpec_C_64fc*)spec_c, &ipp_sz ));
1587                 }
1588                 spec = spec_c;
1589             }
1590
1591             sz += ipp_sz;
1592         }
1593         else
1594 #endif
1595         {
1596             if( len != prev_len )
1597                 nf = DFTFactorize( len, factors );
1598
1599             inplace_transform = factors[0] == factors[nf-1];
1600             sz += len*(complex_elem_size + sizeof(int));
1601             i = nf > 1 && (factors[0] & 1) == 0;
1602             if( (factors[i] & 1) != 0 && factors[i] > 5 )
1603                 sz += (factors[i]+1)*complex_elem_size;
1604
1605             if( (stage == 0 && ((src.data == dst.data && !inplace_transform) || odd_real)) ||
1606                 (stage == 1 && !inplace_transform) )
1607             {
1608                 use_buf = 1;
1609                 sz += len*complex_elem_size;
1610             }
1611         }
1612
1613         ptr = (uchar*)buf;
1614         buf.allocate( sz + 32 );
1615         if( ptr != (uchar*)buf )
1616             prev_len = 0; // because we release the buffer,
1617                           // force recalculation of
1618                           // twiddle factors and permutation table
1619         ptr = (uchar*)buf;
1620         if( !spec )
1621         {
1622             wave = ptr;
1623             ptr += len*complex_elem_size;
1624             itab = (int*)ptr;
1625             ptr = (uchar*)cvAlignPtr( ptr + len*sizeof(int), 16 );
1626
1627             if( len != prev_len || (!inplace_transform && inv && real_transform))
1628                 DFTInit( len, nf, factors, itab, complex_elem_size,
1629                             wave, stage == 0 && inv && real_transform );
1630             // otherwise reuse the tables calculated on the previous stage
1631         }
1632
1633         if( stage == 0 )
1634         {
1635             uchar* tmp_buf = 0;
1636             int dptr_offset = 0;
1637             int dst_full_len = len*elem_size;
1638             int _flags = (int)inv + (src.channels() != dst.channels() ?
1639                          DFT_COMPLEX_INPUT_OR_OUTPUT : 0);
1640             if( use_buf )
1641             {
1642                 tmp_buf = ptr;
1643                 ptr += len*complex_elem_size;
1644                 if( odd_real && !inv && len > 1 &&
1645                     !(_flags & DFT_COMPLEX_INPUT_OR_OUTPUT))
1646                     dptr_offset = elem_size;
1647             }
1648
1649             if( !inv && (_flags & DFT_COMPLEX_INPUT_OR_OUTPUT) )
1650                 dst_full_len += (len & 1) ? elem_size : complex_elem_size;
1651
1652             dft_func = dft_tbl[(!real_transform ? 0 : !inv ? 1 : 2) + (depth == CV_64F)*3];
1653
1654             if( count > 1 && !(flags & DFT_ROWS) && (!inv || !real_transform) )
1655                 stage = 1;
1656             else if( flags & CV_DXT_SCALE )
1657                 scale = 1./(len * (flags & DFT_ROWS ? 1 : count));
1658
1659             if( nonzero_rows <= 0 || nonzero_rows > count )
1660                 nonzero_rows = count;
1661
1662             for( i = 0; i < nonzero_rows; i++ )
1663             {
1664                 uchar* sptr = src.data + i*src.step;
1665                 uchar* dptr0 = dst.data + i*dst.step;
1666                 uchar* dptr = dptr0;
1667
1668                 if( tmp_buf )
1669                     dptr = tmp_buf;
1670
1671                 dft_func( sptr, dptr, len, nf, factors, itab, wave, len, spec, ptr, _flags, scale );
1672                 if( dptr != dptr0 )
1673                     memcpy( dptr0, dptr + dptr_offset, dst_full_len );
1674             }
1675
1676             for( ; i < count; i++ )
1677             {
1678                 uchar* dptr0 = dst.data + i*dst.step;
1679                 memset( dptr0, 0, dst_full_len );
1680             }
1681
1682             if( stage != 1 )
1683                 break;
1684             src = dst;
1685         }
1686         else
1687         {
1688             int a = 0, b = count;
1689             uchar *buf0, *buf1, *dbuf0, *dbuf1;
1690             uchar* sptr0 = src.data;
1691             uchar* dptr0 = dst.data;
1692             buf0 = ptr;
1693             ptr += len*complex_elem_size;
1694             buf1 = ptr;
1695             ptr += len*complex_elem_size;
1696             dbuf0 = buf0, dbuf1 = buf1;
1697
1698             if( use_buf )
1699             {
1700                 dbuf1 = ptr;
1701                 dbuf0 = buf1;
1702                 ptr += len*complex_elem_size;
1703             }
1704
1705             dft_func = dft_tbl[(depth == CV_64F)*3];
1706
1707             if( real_transform && inv && src.cols > 1 )
1708                 stage = 0;
1709             else if( flags & CV_DXT_SCALE )
1710                 scale = 1./(len * count);
1711
1712             if( real_transform )
1713             {
1714                 int even;
1715                 a = 1;
1716                 even = (count & 1) == 0;
1717                 b = (count+1)/2;
1718                 if( !inv )
1719                 {
1720                     memset( buf0, 0, len*complex_elem_size );
1721                     CopyColumn( sptr0, src.step, buf0, complex_elem_size, len, elem_size );
1722                     sptr0 += dst.channels()*elem_size;
1723                     if( even )
1724                     {
1725                         memset( buf1, 0, len*complex_elem_size );
1726                         CopyColumn( sptr0 + (count-2)*elem_size, src.step,
1727                                     buf1, complex_elem_size, len, elem_size );
1728                     }
1729                 }
1730                 else if( src.channels() == 1 )
1731                 {
1732                     CopyColumn( sptr0, src.step, buf0, elem_size, len, elem_size );
1733                     ExpandCCS( buf0, len, elem_size );
1734                     if( even )
1735                     {
1736                         CopyColumn( sptr0 + (count-1)*elem_size, src.step,
1737                                     buf1, elem_size, len, elem_size );
1738                         ExpandCCS( buf1, len, elem_size );
1739                     }
1740                     sptr0 += elem_size;
1741                 }
1742                 else
1743                 {
1744                     CopyColumn( sptr0, src.step, buf0, complex_elem_size, len, complex_elem_size );
1745                     if( even )
1746                     {
1747                         CopyColumn( sptr0 + b*complex_elem_size, src.step,
1748                                        buf1, complex_elem_size, len, complex_elem_size );
1749                     }
1750                     sptr0 += complex_elem_size;
1751                 }
1752
1753                 if( even )
1754                     dft_func( buf1, dbuf1, len, nf, factors, itab,
1755                               wave, len, spec, ptr, inv, scale );
1756                 dft_func( buf0, dbuf0, len, nf, factors, itab,
1757                           wave, len, spec, ptr, inv, scale );
1758
1759                 if( dst.channels() == 1 )
1760                 {
1761                     if( !inv )
1762                     {
1763                         // copy the half of output vector to the first/last column.
1764                         // before doing that, defgragment the vector
1765                         memcpy( dbuf0 + elem_size, dbuf0, elem_size );
1766                         CopyColumn( dbuf0 + elem_size, elem_size, dptr0,
1767                                        dst.step, len, elem_size );
1768                         if( even )
1769                         {
1770                             memcpy( dbuf1 + elem_size, dbuf1, elem_size );
1771                             CopyColumn( dbuf1 + elem_size, elem_size,
1772                                            dptr0 + (count-1)*elem_size,
1773                                            dst.step, len, elem_size );
1774                         }
1775                         dptr0 += elem_size;
1776                     }
1777                     else
1778                     {
1779                         // copy the real part of the complex vector to the first/last column
1780                         CopyColumn( dbuf0, complex_elem_size, dptr0, dst.step, len, elem_size );
1781                         if( even )
1782                             CopyColumn( dbuf1, complex_elem_size, dptr0 + (count-1)*elem_size,
1783                                            dst.step, len, elem_size );
1784                         dptr0 += elem_size;
1785                     }
1786                 }
1787                 else
1788                 {
1789                     assert( !inv );
1790                     CopyColumn( dbuf0, complex_elem_size, dptr0,
1791                                    dst.step, len, complex_elem_size );
1792                     if( even )
1793                         CopyColumn( dbuf1, complex_elem_size,
1794                                        dptr0 + b*complex_elem_size,
1795                                        dst.step, len, complex_elem_size );
1796                     dptr0 += complex_elem_size;
1797                 }
1798             }
1799
1800             for( i = a; i < b; i += 2 )
1801             {
1802                 if( i+1 < b )
1803                 {
1804                     CopyFrom2Columns( sptr0, src.step, buf0, buf1, len, complex_elem_size );
1805                     dft_func( buf1, dbuf1, len, nf, factors, itab,
1806                               wave, len, spec, ptr, inv, scale );
1807                 }
1808                 else
1809                     CopyColumn( sptr0, src.step, buf0, complex_elem_size, len, complex_elem_size );
1810
1811                 dft_func( buf0, dbuf0, len, nf, factors, itab,
1812                           wave, len, spec, ptr, inv, scale );
1813
1814                 if( i+1 < b )
1815                     CopyTo2Columns( dbuf0, dbuf1, dptr0, dst.step, len, complex_elem_size );
1816                 else
1817                     CopyColumn( dbuf0, complex_elem_size, dptr0, dst.step, len, complex_elem_size );
1818                 sptr0 += 2*complex_elem_size;
1819                 dptr0 += 2*complex_elem_size;
1820             }
1821
1822             if( stage != 0 )
1823             {
1824                 if( !inv && real_transform && dst.channels() == 2 && len > 1 )
1825                 {
1826                     int n = dst.cols;
1827                     if( elem_size == (int)sizeof(float) )
1828                     {
1829                         float* p0 = (float*)dst.data;
1830                         size_t dstep = dst.step/sizeof(p0[0]);
1831                         for( i = 0; i < len; i++ )
1832                         {
1833                             float* p = p0 + dstep*i;
1834                             float* q = i == 0 || i*2 == len ? p : p0 + dstep*(len-i);
1835
1836                             for( int j = 1; j < (n+1)/2; j++ )
1837                             {
1838                                 p[(n-j)*2] = q[j*2];
1839                                 p[(n-j)*2+1] = -q[j*2+1];
1840                             }
1841                         }
1842                     }
1843                     else
1844                     {
1845                         double* p0 = (double*)dst.data;
1846                         size_t dstep = dst.step/sizeof(p0[0]);
1847                         for( i = 0; i < len; i++ )
1848                         {
1849                             double* p = p0 + dstep*i;
1850                             double* q = i == 0 || i*2 == len ? p : p0 + dstep*(len-i);
1851
1852                             for( int j = 1; j < (n+1)/2; j++ )
1853                             {
1854                                 p[(n-j)*2] = q[j*2];
1855                                 p[(n-j)*2+1] = -q[j*2+1];
1856                             }
1857                         }
1858                     }
1859                 }
1860                 break;
1861             }
1862             src = dst;
1863         }
1864     }
1865
1866 #ifdef HAVE_IPP
1867     if( spec_c )
1868     {
1869         if( depth == CV_32F )
1870             ippsDFTFree_C_32fc( (IppsDFTSpec_C_32fc*)spec_c );
1871         else
1872             ippsDFTFree_C_64fc( (IppsDFTSpec_C_64fc*)spec_c );
1873     }
1874
1875     if( spec_r )
1876     {
1877         if( depth == CV_32F )
1878             ippsDFTFree_R_32f( (IppsDFTSpec_R_32f*)spec_r );
1879         else
1880             ippsDFTFree_R_64f( (IppsDFTSpec_R_64f*)spec_r );
1881     }
1882 #endif
1883 }
1884
1885
1886 void cv::idft( InputArray src, OutputArray dst, int flags, int nonzero_rows )
1887 {
1888     dft( src, dst, flags | DFT_INVERSE, nonzero_rows );
1889 }
1890
1891 void cv::mulSpectrums( InputArray _srcA, InputArray _srcB,
1892                        OutputArray _dst, int flags, bool conjB )
1893 {
1894     Mat srcA = _srcA.getMat(), srcB = _srcB.getMat();
1895     int depth = srcA.depth(), cn = srcA.channels(), type = srcA.type();
1896     int rows = srcA.rows, cols = srcA.cols;
1897     int j, k;
1898
1899     CV_Assert( type == srcB.type() && srcA.size() == srcB.size() );
1900     CV_Assert( type == CV_32FC1 || type == CV_32FC2 || type == CV_64FC1 || type == CV_64FC2 );
1901
1902     _dst.create( srcA.rows, srcA.cols, type );
1903     Mat dst = _dst.getMat();
1904
1905     bool is_1d = (flags & DFT_ROWS) || (rows == 1 || (cols == 1 &&
1906              srcA.isContinuous() && srcB.isContinuous() && dst.isContinuous()));
1907
1908     if( is_1d && !(flags & DFT_ROWS) )
1909         cols = cols + rows - 1, rows = 1;
1910
1911     int ncols = cols*cn;
1912     int j0 = cn == 1;
1913     int j1 = ncols - (cols % 2 == 0 && cn == 1);
1914
1915     if( depth == CV_32F )
1916     {
1917         const float* dataA = (const float*)srcA.data;
1918         const float* dataB = (const float*)srcB.data;
1919         float* dataC = (float*)dst.data;
1920
1921         size_t stepA = srcA.step/sizeof(dataA[0]);
1922         size_t stepB = srcB.step/sizeof(dataB[0]);
1923         size_t stepC = dst.step/sizeof(dataC[0]);
1924
1925         if( !is_1d && cn == 1 )
1926         {
1927             for( k = 0; k < (cols % 2 ? 1 : 2); k++ )
1928             {
1929                 if( k == 1 )
1930                     dataA += cols - 1, dataB += cols - 1, dataC += cols - 1;
1931                 dataC[0] = dataA[0]*dataB[0];
1932                 if( rows % 2 == 0 )
1933                     dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA]*dataB[(rows-1)*stepB];
1934                 if( !conjB )
1935                     for( j = 1; j <= rows - 2; j += 2 )
1936                     {
1937                         double re = (double)dataA[j*stepA]*dataB[j*stepB] -
1938                                     (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
1939                         double im = (double)dataA[j*stepA]*dataB[(j+1)*stepB] +
1940                                     (double)dataA[(j+1)*stepA]*dataB[j*stepB];
1941                         dataC[j*stepC] = (float)re; dataC[(j+1)*stepC] = (float)im;
1942                     }
1943                 else
1944                     for( j = 1; j <= rows - 2; j += 2 )
1945                     {
1946                         double re = (double)dataA[j*stepA]*dataB[j*stepB] +
1947                                     (double)dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
1948                         double im = (double)dataA[(j+1)*stepA]*dataB[j*stepB] -
1949                                     (double)dataA[j*stepA]*dataB[(j+1)*stepB];
1950                         dataC[j*stepC] = (float)re; dataC[(j+1)*stepC] = (float)im;
1951                     }
1952                 if( k == 1 )
1953                     dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1;
1954             }
1955         }
1956
1957         for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC )
1958         {
1959             if( is_1d && cn == 1 )
1960             {
1961                 dataC[0] = dataA[0]*dataB[0];
1962                 if( cols % 2 == 0 )
1963                     dataC[j1] = dataA[j1]*dataB[j1];
1964             }
1965
1966             if( !conjB )
1967                 for( j = j0; j < j1; j += 2 )
1968                 {
1969                     double re = (double)dataA[j]*dataB[j] - (double)dataA[j+1]*dataB[j+1];
1970                     double im = (double)dataA[j+1]*dataB[j] + (double)dataA[j]*dataB[j+1];
1971                     dataC[j] = (float)re; dataC[j+1] = (float)im;
1972                 }
1973             else
1974                 for( j = j0; j < j1; j += 2 )
1975                 {
1976                     double re = (double)dataA[j]*dataB[j] + (double)dataA[j+1]*dataB[j+1];
1977                     double im = (double)dataA[j+1]*dataB[j] - (double)dataA[j]*dataB[j+1];
1978                     dataC[j] = (float)re; dataC[j+1] = (float)im;
1979                 }
1980         }
1981     }
1982     else
1983     {
1984         const double* dataA = (const double*)srcA.data;
1985         const double* dataB = (const double*)srcB.data;
1986         double* dataC = (double*)dst.data;
1987
1988         size_t stepA = srcA.step/sizeof(dataA[0]);
1989         size_t stepB = srcB.step/sizeof(dataB[0]);
1990         size_t stepC = dst.step/sizeof(dataC[0]);
1991
1992         if( !is_1d && cn == 1 )
1993         {
1994             for( k = 0; k < (cols % 2 ? 1 : 2); k++ )
1995             {
1996                 if( k == 1 )
1997                     dataA += cols - 1, dataB += cols - 1, dataC += cols - 1;
1998                 dataC[0] = dataA[0]*dataB[0];
1999                 if( rows % 2 == 0 )
2000                     dataC[(rows-1)*stepC] = dataA[(rows-1)*stepA]*dataB[(rows-1)*stepB];
2001                 if( !conjB )
2002                     for( j = 1; j <= rows - 2; j += 2 )
2003                     {
2004                         double re = dataA[j*stepA]*dataB[j*stepB] -
2005                                     dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
2006                         double im = dataA[j*stepA]*dataB[(j+1)*stepB] +
2007                                     dataA[(j+1)*stepA]*dataB[j*stepB];
2008                         dataC[j*stepC] = re; dataC[(j+1)*stepC] = im;
2009                     }
2010                 else
2011                     for( j = 1; j <= rows - 2; j += 2 )
2012                     {
2013                         double re = dataA[j*stepA]*dataB[j*stepB] +
2014                                     dataA[(j+1)*stepA]*dataB[(j+1)*stepB];
2015                         double im = dataA[(j+1)*stepA]*dataB[j*stepB] -
2016                                     dataA[j*stepA]*dataB[(j+1)*stepB];
2017                         dataC[j*stepC] = re; dataC[(j+1)*stepC] = im;
2018                     }
2019                 if( k == 1 )
2020                     dataA -= cols - 1, dataB -= cols - 1, dataC -= cols - 1;
2021             }
2022         }
2023
2024         for( ; rows--; dataA += stepA, dataB += stepB, dataC += stepC )
2025         {
2026             if( is_1d && cn == 1 )
2027             {
2028                 dataC[0] = dataA[0]*dataB[0];
2029                 if( cols % 2 == 0 )
2030                     dataC[j1] = dataA[j1]*dataB[j1];
2031             }
2032
2033             if( !conjB )
2034                 for( j = j0; j < j1; j += 2 )
2035                 {
2036                     double re = dataA[j]*dataB[j] - dataA[j+1]*dataB[j+1];
2037                     double im = dataA[j+1]*dataB[j] + dataA[j]*dataB[j+1];
2038                     dataC[j] = re; dataC[j+1] = im;
2039                 }
2040             else
2041                 for( j = j0; j < j1; j += 2 )
2042                 {
2043                     double re = dataA[j]*dataB[j] + dataA[j+1]*dataB[j+1];
2044                     double im = dataA[j+1]*dataB[j] - dataA[j]*dataB[j+1];
2045                     dataC[j] = re; dataC[j+1] = im;
2046                 }
2047         }
2048     }
2049 }
2050
2051
2052 /****************************************************************************************\
2053                                Discrete Cosine Transform
2054 \****************************************************************************************/
2055
2056 namespace cv
2057 {
2058
2059 /* DCT is calculated using DFT, as described here:
2060    http://www.ece.utexas.edu/~bevans/courses/ee381k/lectures/09_DCT/lecture9/:
2061 */
2062 template<typename T> static void
2063 DCT( const T* src, int src_step, T* dft_src, T* dft_dst, T* dst, int dst_step,
2064      int n, int nf, int* factors, const int* itab, const Complex<T>* dft_wave,
2065      const Complex<T>* dct_wave, const void* spec, Complex<T>* buf )
2066 {
2067     static const T sin_45 = (T)0.70710678118654752440084436210485;
2068     int j, n2 = n >> 1;
2069
2070     src_step /= sizeof(src[0]);
2071     dst_step /= sizeof(dst[0]);
2072     T* dst1 = dst + (n-1)*dst_step;
2073
2074     if( n == 1 )
2075     {
2076         dst[0] = src[0];
2077         return;
2078     }
2079
2080     for( j = 0; j < n2; j++, src += src_step*2 )
2081     {
2082         dft_src[j] = src[0];
2083         dft_src[n-j-1] = src[src_step];
2084     }
2085
2086     RealDFT( dft_src, dft_dst, n, nf, factors,
2087              itab, dft_wave, n, spec, buf, 0, 1.0 );
2088     src = dft_dst;
2089
2090     dst[0] = (T)(src[0]*dct_wave->re*sin_45);
2091     dst += dst_step;
2092     for( j = 1, dct_wave++; j < n2; j++, dct_wave++,
2093                                     dst += dst_step, dst1 -= dst_step )
2094     {
2095         T t0 = dct_wave->re*src[j*2-1] - dct_wave->im*src[j*2];
2096         T t1 = -dct_wave->im*src[j*2-1] - dct_wave->re*src[j*2];
2097         dst[0] = t0;
2098         dst1[0] = t1;
2099     }
2100
2101     dst[0] = src[n-1]*dct_wave->re;
2102 }
2103
2104
2105 template<typename T> static void
2106 IDCT( const T* src, int src_step, T* dft_src, T* dft_dst, T* dst, int dst_step,
2107       int n, int nf, int* factors, const int* itab, const Complex<T>* dft_wave,
2108       const Complex<T>* dct_wave, const void* spec, Complex<T>* buf )
2109 {
2110     static const T sin_45 = (T)0.70710678118654752440084436210485;
2111     int j, n2 = n >> 1;
2112
2113     src_step /= sizeof(src[0]);
2114     dst_step /= sizeof(dst[0]);
2115     const T* src1 = src + (n-1)*src_step;
2116
2117     if( n == 1 )
2118     {
2119         dst[0] = src[0];
2120         return;
2121     }
2122
2123     dft_src[0] = (T)(src[0]*2*dct_wave->re*sin_45);
2124     src += src_step;
2125     for( j = 1, dct_wave++; j < n2; j++, dct_wave++,
2126                                     src += src_step, src1 -= src_step )
2127     {
2128         T t0 = dct_wave->re*src[0] - dct_wave->im*src1[0];
2129         T t1 = -dct_wave->im*src[0] - dct_wave->re*src1[0];
2130         dft_src[j*2-1] = t0;
2131         dft_src[j*2] = t1;
2132     }
2133
2134     dft_src[n-1] = (T)(src[0]*2*dct_wave->re);
2135     CCSIDFT( dft_src, dft_dst, n, nf, factors, itab,
2136              dft_wave, n, spec, buf, 0, 1.0 );
2137
2138     for( j = 0; j < n2; j++, dst += dst_step*2 )
2139     {
2140         dst[0] = dft_dst[j];
2141         dst[dst_step] = dft_dst[n-j-1];
2142     }
2143 }
2144
2145
2146 static void
2147 DCTInit( int n, int elem_size, void* _wave, int inv )
2148 {
2149     static const double DctScale[] =
2150     {
2151     0.707106781186547570, 0.500000000000000000, 0.353553390593273790,
2152     0.250000000000000000, 0.176776695296636890, 0.125000000000000000,
2153     0.088388347648318447, 0.062500000000000000, 0.044194173824159223,
2154     0.031250000000000000, 0.022097086912079612, 0.015625000000000000,
2155     0.011048543456039806, 0.007812500000000000, 0.005524271728019903,
2156     0.003906250000000000, 0.002762135864009952, 0.001953125000000000,
2157     0.001381067932004976, 0.000976562500000000, 0.000690533966002488,
2158     0.000488281250000000, 0.000345266983001244, 0.000244140625000000,
2159     0.000172633491500622, 0.000122070312500000, 0.000086316745750311,
2160     0.000061035156250000, 0.000043158372875155, 0.000030517578125000
2161     };
2162
2163     int i;
2164     Complex<double> w, w1;
2165     double t, scale;
2166
2167     if( n == 1 )
2168         return;
2169
2170     assert( (n&1) == 0 );
2171
2172     if( (n & (n - 1)) == 0 )
2173     {
2174         int m;
2175         for( m = 0; (unsigned)(1 << m) < (unsigned)n; m++ )
2176             ;
2177         scale = (!inv ? 2 : 1)*DctScale[m];
2178         w1.re = DFTTab[m+2][0];
2179         w1.im = -DFTTab[m+2][1];
2180     }
2181     else
2182     {
2183         t = 1./(2*n);
2184         scale = (!inv ? 2 : 1)*std::sqrt(t);
2185         w1.im = sin(-CV_PI*t);
2186         w1.re = std::sqrt(1. - w1.im*w1.im);
2187     }
2188     n >>= 1;
2189
2190     if( elem_size == sizeof(Complex<double>) )
2191     {
2192         Complex<double>* wave = (Complex<double>*)_wave;
2193
2194         w.re = scale;
2195         w.im = 0.;
2196
2197         for( i = 0; i <= n; i++ )
2198         {
2199             wave[i] = w;
2200             t = w.re*w1.re - w.im*w1.im;
2201             w.im = w.re*w1.im + w.im*w1.re;
2202             w.re = t;
2203         }
2204     }
2205     else
2206     {
2207         Complex<float>* wave = (Complex<float>*)_wave;
2208         assert( elem_size == sizeof(Complex<float>) );
2209
2210         w.re = (float)scale;
2211         w.im = 0.f;
2212
2213         for( i = 0; i <= n; i++ )
2214         {
2215             wave[i].re = (float)w.re;
2216             wave[i].im = (float)w.im;
2217             t = w.re*w1.re - w.im*w1.im;
2218             w.im = w.re*w1.im + w.im*w1.re;
2219             w.re = t;
2220         }
2221     }
2222 }
2223
2224
2225 typedef void (*DCTFunc)(const void* src, int src_step, void* dft_src,
2226                         void* dft_dst, void* dst, int dst_step, int n,
2227                         int nf, int* factors, const int* itab, const void* dft_wave,
2228                         const void* dct_wave, const void* spec, void* buf );
2229
2230 static void DCT_32f(const float* src, int src_step, float* dft_src, float* dft_dst,
2231                     float* dst, int dst_step, int n, int nf, int* factors, const int* itab,
2232                     const Complexf* dft_wave, const Complexf* dct_wave, const void* spec, Complexf* buf )
2233 {
2234     DCT(src, src_step, dft_src, dft_dst, dst, dst_step,
2235         n, nf, factors, itab, dft_wave, dct_wave, spec, buf);
2236 }
2237
2238 static void IDCT_32f(const float* src, int src_step, float* dft_src, float* dft_dst,
2239                     float* dst, int dst_step, int n, int nf, int* factors, const int* itab,
2240                     const Complexf* dft_wave, const Complexf* dct_wave, const void* spec, Complexf* buf )
2241 {
2242     IDCT(src, src_step, dft_src, dft_dst, dst, dst_step,
2243          n, nf, factors, itab, dft_wave, dct_wave, spec, buf);
2244 }
2245
2246 static void DCT_64f(const double* src, int src_step, double* dft_src, double* dft_dst,
2247                     double* dst, int dst_step, int n, int nf, int* factors, const int* itab,
2248                     const Complexd* dft_wave, const Complexd* dct_wave, const void* spec, Complexd* buf )
2249 {
2250     DCT(src, src_step, dft_src, dft_dst, dst, dst_step,
2251         n, nf, factors, itab, dft_wave, dct_wave, spec, buf);
2252 }
2253
2254 static void IDCT_64f(const double* src, int src_step, double* dft_src, double* dft_dst,
2255                      double* dst, int dst_step, int n, int nf, int* factors, const int* itab,
2256                      const Complexd* dft_wave, const Complexd* dct_wave, const void* spec, Complexd* buf )
2257 {
2258     IDCT(src, src_step, dft_src, dft_dst, dst, dst_step,
2259          n, nf, factors, itab, dft_wave, dct_wave, spec, buf);
2260 }
2261
2262 }
2263
2264 void cv::dct( InputArray _src0, OutputArray _dst, int flags )
2265 {
2266     static DCTFunc dct_tbl[4] =
2267     {
2268         (DCTFunc)DCT_32f,
2269         (DCTFunc)IDCT_32f,
2270         (DCTFunc)DCT_64f,
2271         (DCTFunc)IDCT_64f
2272     };
2273
2274     bool inv = (flags & DCT_INVERSE) != 0;
2275     Mat src0 = _src0.getMat(), src = src0;
2276     int type = src.type(), depth = src.depth();
2277     void /* *spec_dft = 0, */ *spec = 0;
2278
2279     double scale = 1.;
2280     int prev_len = 0, nf = 0, stage, end_stage;
2281     uchar *src_dft_buf = 0, *dst_dft_buf = 0;
2282     uchar *dft_wave = 0, *dct_wave = 0;
2283     int* itab = 0;
2284     uchar* ptr = 0;
2285     int elem_size = (int)src.elemSize(), complex_elem_size = elem_size*2;
2286     int factors[34], inplace_transform;
2287     int i, len, count;
2288     AutoBuffer<uchar> buf;
2289
2290     CV_Assert( type == CV_32FC1 || type == CV_64FC1 );
2291     _dst.create( src.rows, src.cols, type );
2292     Mat dst = _dst.getMat();
2293
2294     DCTFunc dct_func = dct_tbl[(int)inv + (depth == CV_64F)*2];
2295
2296     if( (flags & DFT_ROWS) || src.rows == 1 ||
2297         (src.cols == 1 && (src.isContinuous() && dst.isContinuous())))
2298     {
2299         stage = end_stage = 0;
2300     }
2301     else
2302     {
2303         stage = src.cols == 1;
2304         end_stage = 1;
2305     }
2306
2307     for( ; stage <= end_stage; stage++ )
2308     {
2309         uchar *sptr = src.data, *dptr = dst.data;
2310         size_t sstep0, sstep1, dstep0, dstep1;
2311
2312         if( stage == 0 )
2313         {
2314             len = src.cols;
2315             count = src.rows;
2316             if( len == 1 && !(flags & DFT_ROWS) )
2317             {
2318                 len = src.rows;
2319                 count = 1;
2320             }
2321             sstep0 = src.step;
2322             dstep0 = dst.step;
2323             sstep1 = dstep1 = elem_size;
2324         }
2325         else
2326         {
2327             len = dst.rows;
2328             count = dst.cols;
2329             sstep1 = src.step;
2330             dstep1 = dst.step;
2331             sstep0 = dstep0 = elem_size;
2332         }
2333
2334         if( len != prev_len )
2335         {
2336             int sz;
2337
2338             if( len > 1 && (len & 1) )
2339                 CV_Error( CV_StsNotImplemented, "Odd-size DCT\'s are not implemented" );
2340
2341             sz = len*elem_size;
2342             sz += (len/2 + 1)*complex_elem_size;
2343
2344             spec = 0;
2345             inplace_transform = 1;
2346             /*if( len*count >= 64 && DFTInitAlloc_R_32f_p )
2347             {
2348                 int ipp_sz = 0;
2349                 if( depth == CV_32F )
2350                 {
2351                     if( spec_dft )
2352                         IPPI_CALL( DFTFree_R_32f_p( spec_dft ));
2353                     IPPI_CALL( DFTInitAlloc_R_32f_p( &spec_dft, len, 8, cvAlgHintNone ));
2354                     IPPI_CALL( DFTGetBufSize_R_32f_p( spec_dft, &ipp_sz ));
2355                 }
2356                 else
2357                 {
2358                     if( spec_dft )
2359                         IPPI_CALL( DFTFree_R_64f_p( spec_dft ));
2360                     IPPI_CALL( DFTInitAlloc_R_64f_p( &spec_dft, len, 8, cvAlgHintNone ));
2361                     IPPI_CALL( DFTGetBufSize_R_64f_p( spec_dft, &ipp_sz ));
2362                 }
2363                 spec = spec_dft;
2364                 sz += ipp_sz;
2365             }
2366             else*/
2367             {
2368                 sz += len*(complex_elem_size + sizeof(int)) + complex_elem_size;
2369
2370                 nf = DFTFactorize( len, factors );
2371                 inplace_transform = factors[0] == factors[nf-1];
2372
2373                 i = nf > 1 && (factors[0] & 1) == 0;
2374                 if( (factors[i] & 1) != 0 && factors[i] > 5 )
2375                     sz += (factors[i]+1)*complex_elem_size;
2376
2377                 if( !inplace_transform )
2378                     sz += len*elem_size;
2379             }
2380
2381             buf.allocate( sz + 32 );
2382             ptr = (uchar*)buf;
2383
2384             if( !spec )
2385             {
2386                 dft_wave = ptr;
2387                 ptr += len*complex_elem_size;
2388                 itab = (int*)ptr;
2389                 ptr = (uchar*)cvAlignPtr( ptr + len*sizeof(int), 16 );
2390                 DFTInit( len, nf, factors, itab, complex_elem_size, dft_wave, inv );
2391             }
2392
2393             dct_wave = ptr;
2394             ptr += (len/2 + 1)*complex_elem_size;
2395             src_dft_buf = dst_dft_buf = ptr;
2396             ptr += len*elem_size;
2397             if( !inplace_transform )
2398             {
2399                 dst_dft_buf = ptr;
2400                 ptr += len*elem_size;
2401             }
2402             DCTInit( len, complex_elem_size, dct_wave, inv );
2403             if( !inv )
2404                 scale += scale;
2405             prev_len = len;
2406         }
2407         // otherwise reuse the tables calculated on the previous stage
2408         for( i = 0; i < count; i++ )
2409         {
2410             dct_func( sptr + i*sstep0, (int)sstep1, src_dft_buf, dst_dft_buf,
2411                       dptr + i*dstep0, (int)dstep1, len, nf, factors,
2412                       itab, dft_wave, dct_wave, spec, ptr );
2413         }
2414         src = dst;
2415     }
2416 }
2417
2418
2419 void cv::idct( InputArray src, OutputArray dst, int flags )
2420 {
2421     dct( src, dst, flags | DCT_INVERSE );
2422 }
2423
2424 namespace cv
2425 {
2426
2427 static const int optimalDFTSizeTab[] = {
2428 1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48,
2429 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 125, 128, 135, 144, 150, 160,
2430 162, 180, 192, 200, 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375,
2431 384, 400, 405, 432, 450, 480, 486, 500, 512, 540, 576, 600, 625, 640, 648, 675, 720,
2432 729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 1080, 1125, 1152, 1200,
2433 1215, 1250, 1280, 1296, 1350, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875,
2434 1920, 1944, 2000, 2025, 2048, 2160, 2187, 2250, 2304, 2400, 2430, 2500, 2560, 2592,
2435 2700, 2880, 2916, 3000, 3072, 3125, 3200, 3240, 3375, 3456, 3600, 3645, 3750, 3840,
2436 3888, 4000, 4050, 4096, 4320, 4374, 4500, 4608, 4800, 4860, 5000, 5120, 5184, 5400,
2437 5625, 5760, 5832, 6000, 6075, 6144, 6250, 6400, 6480, 6561, 6750, 6912, 7200, 7290,
2438 7500, 7680, 7776, 8000, 8100, 8192, 8640, 8748, 9000, 9216, 9375, 9600, 9720, 10000,
2439 10125, 10240, 10368, 10800, 10935, 11250, 11520, 11664, 12000, 12150, 12288, 12500,
2440 12800, 12960, 13122, 13500, 13824, 14400, 14580, 15000, 15360, 15552, 15625, 16000,
2441 16200, 16384, 16875, 17280, 17496, 18000, 18225, 18432, 18750, 19200, 19440, 19683,
2442 20000, 20250, 20480, 20736, 21600, 21870, 22500, 23040, 23328, 24000, 24300, 24576,
2443 25000, 25600, 25920, 26244, 27000, 27648, 28125, 28800, 29160, 30000, 30375, 30720,
2444 31104, 31250, 32000, 32400, 32768, 32805, 33750, 34560, 34992, 36000, 36450, 36864,
2445 37500, 38400, 38880, 39366, 40000, 40500, 40960, 41472, 43200, 43740, 45000, 46080,
2446 46656, 46875, 48000, 48600, 49152, 50000, 50625, 51200, 51840, 52488, 54000, 54675,
2447 55296, 56250, 57600, 58320, 59049, 60000, 60750, 61440, 62208, 62500, 64000, 64800,
2448 65536, 65610, 67500, 69120, 69984, 72000, 72900, 73728, 75000, 76800, 77760, 78125,
2449 78732, 80000, 81000, 81920, 82944, 84375, 86400, 87480, 90000, 91125, 92160, 93312,
2450 93750, 96000, 97200, 98304, 98415, 100000, 101250, 102400, 103680, 104976, 108000,
2451 109350, 110592, 112500, 115200, 116640, 118098, 120000, 121500, 122880, 124416, 125000,
2452 128000, 129600, 131072, 131220, 135000, 138240, 139968, 140625, 144000, 145800, 147456,
2453 150000, 151875, 153600, 155520, 156250, 157464, 160000, 162000, 163840, 164025, 165888,
2454 168750, 172800, 174960, 177147, 180000, 182250, 184320, 186624, 187500, 192000, 194400,
2455 196608, 196830, 200000, 202500, 204800, 207360, 209952, 216000, 218700, 221184, 225000,
2456 230400, 233280, 234375, 236196, 240000, 243000, 245760, 248832, 250000, 253125, 256000,
2457 259200, 262144, 262440, 270000, 273375, 276480, 279936, 281250, 288000, 291600, 294912,
2458 295245, 300000, 303750, 307200, 311040, 312500, 314928, 320000, 324000, 327680, 328050,
2459 331776, 337500, 345600, 349920, 354294, 360000, 364500, 368640, 373248, 375000, 384000,
2460 388800, 390625, 393216, 393660, 400000, 405000, 409600, 414720, 419904, 421875, 432000,
2461 437400, 442368, 450000, 455625, 460800, 466560, 468750, 472392, 480000, 486000, 491520,
2462 492075, 497664, 500000, 506250, 512000, 518400, 524288, 524880, 531441, 540000, 546750,
2463 552960, 559872, 562500, 576000, 583200, 589824, 590490, 600000, 607500, 614400, 622080,
2464 625000, 629856, 640000, 648000, 655360, 656100, 663552, 675000, 691200, 699840, 703125,
2465 708588, 720000, 729000, 737280, 746496, 750000, 759375, 768000, 777600, 781250, 786432,
2466 787320, 800000, 810000, 819200, 820125, 829440, 839808, 843750, 864000, 874800, 884736,
2467 885735, 900000, 911250, 921600, 933120, 937500, 944784, 960000, 972000, 983040, 984150,
2468 995328, 1000000, 1012500, 1024000, 1036800, 1048576, 1049760, 1062882, 1080000, 1093500,
2469 1105920, 1119744, 1125000, 1152000, 1166400, 1171875, 1179648, 1180980, 1200000,
2470 1215000, 1228800, 1244160, 1250000, 1259712, 1265625, 1280000, 1296000, 1310720,
2471 1312200, 1327104, 1350000, 1366875, 1382400, 1399680, 1406250, 1417176, 1440000,
2472 1458000, 1474560, 1476225, 1492992, 1500000, 1518750, 1536000, 1555200, 1562500,
2473 1572864, 1574640, 1594323, 1600000, 1620000, 1638400, 1640250, 1658880, 1679616,
2474 1687500, 1728000, 1749600, 1769472, 1771470, 1800000, 1822500, 1843200, 1866240,
2475 1875000, 1889568, 1920000, 1944000, 1953125, 1966080, 1968300, 1990656, 2000000,
2476 2025000, 2048000, 2073600, 2097152, 2099520, 2109375, 2125764, 2160000, 2187000,
2477 2211840, 2239488, 2250000, 2278125, 2304000, 2332800, 2343750, 2359296, 2361960,
2478 2400000, 2430000, 2457600, 2460375, 2488320, 2500000, 2519424, 2531250, 2560000,
2479 2592000, 2621440, 2624400, 2654208, 2657205, 2700000, 2733750, 2764800, 2799360,
2480 2812500, 2834352, 2880000, 2916000, 2949120, 2952450, 2985984, 3000000, 3037500,
2481 3072000, 3110400, 3125000, 3145728, 3149280, 3188646, 3200000, 3240000, 3276800,
2482 3280500, 3317760, 3359232, 3375000, 3456000, 3499200, 3515625, 3538944, 3542940,
2483 3600000, 3645000, 3686400, 3732480, 3750000, 3779136, 3796875, 3840000, 3888000,
2484 3906250, 3932160, 3936600, 3981312, 4000000, 4050000, 4096000, 4100625, 4147200,
2485 4194304, 4199040, 4218750, 4251528, 4320000, 4374000, 4423680, 4428675, 4478976,
2486 4500000, 4556250, 4608000, 4665600, 4687500, 4718592, 4723920, 4782969, 4800000,
2487 4860000, 4915200, 4920750, 4976640, 5000000, 5038848, 5062500, 5120000, 5184000,
2488 5242880, 5248800, 5308416, 5314410, 5400000, 5467500, 5529600, 5598720, 5625000,
2489 5668704, 5760000, 5832000, 5859375, 5898240, 5904900, 5971968, 6000000, 6075000,
2490 6144000, 6220800, 6250000, 6291456, 6298560, 6328125, 6377292, 6400000, 6480000,
2491 6553600, 6561000, 6635520, 6718464, 6750000, 6834375, 6912000, 6998400, 7031250,
2492 7077888, 7085880, 7200000, 7290000, 7372800, 7381125, 7464960, 7500000, 7558272,
2493 7593750, 7680000, 7776000, 7812500, 7864320, 7873200, 7962624, 7971615, 8000000,
2494 8100000, 8192000, 8201250, 8294400, 8388608, 8398080, 8437500, 8503056, 8640000,
2495 8748000, 8847360, 8857350, 8957952, 9000000, 9112500, 9216000, 9331200, 9375000,
2496 9437184, 9447840, 9565938, 9600000, 9720000, 9765625, 9830400, 9841500, 9953280,
2497 10000000, 10077696, 10125000, 10240000, 10368000, 10485760, 10497600, 10546875, 10616832,
2498 10628820, 10800000, 10935000, 11059200, 11197440, 11250000, 11337408, 11390625, 11520000,
2499 11664000, 11718750, 11796480, 11809800, 11943936, 12000000, 12150000, 12288000, 12301875,
2500 12441600, 12500000, 12582912, 12597120, 12656250, 12754584, 12800000, 12960000, 13107200,
2501 13122000, 13271040, 13286025, 13436928, 13500000, 13668750, 13824000, 13996800, 14062500,
2502 14155776, 14171760, 14400000, 14580000, 14745600, 14762250, 14929920, 15000000, 15116544,
2503 15187500, 15360000, 15552000, 15625000, 15728640, 15746400, 15925248, 15943230, 16000000,
2504 16200000, 16384000, 16402500, 16588800, 16777216, 16796160, 16875000, 17006112, 17280000,
2505 17496000, 17578125, 17694720, 17714700, 17915904, 18000000, 18225000, 18432000, 18662400,
2506 18750000, 18874368, 18895680, 18984375, 19131876, 19200000, 19440000, 19531250, 19660800,
2507 19683000, 19906560, 20000000, 20155392, 20250000, 20480000, 20503125, 20736000, 20971520,
2508 20995200, 21093750, 21233664, 21257640, 21600000, 21870000, 22118400, 22143375, 22394880,
2509 22500000, 22674816, 22781250, 23040000, 23328000, 23437500, 23592960, 23619600, 23887872,
2510 23914845, 24000000, 24300000, 24576000, 24603750, 24883200, 25000000, 25165824, 25194240,
2511 25312500, 25509168, 25600000, 25920000, 26214400, 26244000, 26542080, 26572050, 26873856,
2512 27000000, 27337500, 27648000, 27993600, 28125000, 28311552, 28343520, 28800000, 29160000,
2513 29296875, 29491200, 29524500, 29859840, 30000000, 30233088, 30375000, 30720000, 31104000,
2514 31250000, 31457280, 31492800, 31640625, 31850496, 31886460, 32000000, 32400000, 32768000,
2515 32805000, 33177600, 33554432, 33592320, 33750000, 34012224, 34171875, 34560000, 34992000,
2516 35156250, 35389440, 35429400, 35831808, 36000000, 36450000, 36864000, 36905625, 37324800,
2517 37500000, 37748736, 37791360, 37968750, 38263752, 38400000, 38880000, 39062500, 39321600,
2518 39366000, 39813120, 39858075, 40000000, 40310784, 40500000, 40960000, 41006250, 41472000,
2519 41943040, 41990400, 42187500, 42467328, 42515280, 43200000, 43740000, 44236800, 44286750,
2520 44789760, 45000000, 45349632, 45562500, 46080000, 46656000, 46875000, 47185920, 47239200,
2521 47775744, 47829690, 48000000, 48600000, 48828125, 49152000, 49207500, 49766400, 50000000,
2522 50331648, 50388480, 50625000, 51018336, 51200000, 51840000, 52428800, 52488000, 52734375,
2523 53084160, 53144100, 53747712, 54000000, 54675000, 55296000, 55987200, 56250000, 56623104,
2524 56687040, 56953125, 57600000, 58320000, 58593750, 58982400, 59049000, 59719680, 60000000,
2525 60466176, 60750000, 61440000, 61509375, 62208000, 62500000, 62914560, 62985600, 63281250,
2526 63700992, 63772920, 64000000, 64800000, 65536000, 65610000, 66355200, 66430125, 67108864,
2527 67184640, 67500000, 68024448, 68343750, 69120000, 69984000, 70312500, 70778880, 70858800,
2528 71663616, 72000000, 72900000, 73728000, 73811250, 74649600, 75000000, 75497472, 75582720,
2529 75937500, 76527504, 76800000, 77760000, 78125000, 78643200, 78732000, 79626240, 79716150,
2530 80000000, 80621568, 81000000, 81920000, 82012500, 82944000, 83886080, 83980800, 84375000,
2531 84934656, 85030560, 86400000, 87480000, 87890625, 88473600, 88573500, 89579520, 90000000,
2532 90699264, 91125000, 92160000, 93312000, 93750000, 94371840, 94478400, 94921875, 95551488,
2533 95659380, 96000000, 97200000, 97656250, 98304000, 98415000, 99532800, 100000000,
2534 100663296, 100776960, 101250000, 102036672, 102400000, 102515625, 103680000, 104857600,
2535 104976000, 105468750, 106168320, 106288200, 107495424, 108000000, 109350000, 110592000,
2536 110716875, 111974400, 112500000, 113246208, 113374080, 113906250, 115200000, 116640000,
2537 117187500, 117964800, 118098000, 119439360, 119574225, 120000000, 120932352, 121500000,
2538 122880000, 123018750, 124416000, 125000000, 125829120, 125971200, 126562500, 127401984,
2539 127545840, 128000000, 129600000, 131072000, 131220000, 132710400, 132860250, 134217728,
2540 134369280, 135000000, 136048896, 136687500, 138240000, 139968000, 140625000, 141557760,
2541 141717600, 143327232, 144000000, 145800000, 146484375, 147456000, 147622500, 149299200,
2542 150000000, 150994944, 151165440, 151875000, 153055008, 153600000, 155520000, 156250000,
2543 157286400, 157464000, 158203125, 159252480, 159432300, 160000000, 161243136, 162000000,
2544 163840000, 164025000, 165888000, 167772160, 167961600, 168750000, 169869312, 170061120,
2545 170859375, 172800000, 174960000, 175781250, 176947200, 177147000, 179159040, 180000000,
2546 181398528, 182250000, 184320000, 184528125, 186624000, 187500000, 188743680, 188956800,
2547 189843750, 191102976, 191318760, 192000000, 194400000, 195312500, 196608000, 196830000,
2548 199065600, 199290375, 200000000, 201326592, 201553920, 202500000, 204073344, 204800000,
2549 205031250, 207360000, 209715200, 209952000, 210937500, 212336640, 212576400, 214990848,
2550 216000000, 218700000, 221184000, 221433750, 223948800, 225000000, 226492416, 226748160,
2551 227812500, 230400000, 233280000, 234375000, 235929600, 236196000, 238878720, 239148450,
2552 240000000, 241864704, 243000000, 244140625, 245760000, 246037500, 248832000, 250000000,
2553 251658240, 251942400, 253125000, 254803968, 255091680, 256000000, 259200000, 262144000,
2554 262440000, 263671875, 265420800, 265720500, 268435456, 268738560, 270000000, 272097792,
2555 273375000, 276480000, 279936000, 281250000, 283115520, 283435200, 284765625, 286654464,
2556 288000000, 291600000, 292968750, 294912000, 295245000, 298598400, 300000000, 301989888,
2557 302330880, 303750000, 306110016, 307200000, 307546875, 311040000, 312500000, 314572800,
2558 314928000, 316406250, 318504960, 318864600, 320000000, 322486272, 324000000, 327680000,
2559 328050000, 331776000, 332150625, 335544320, 335923200, 337500000, 339738624, 340122240,
2560 341718750, 345600000, 349920000, 351562500, 353894400, 354294000, 358318080, 360000000,
2561 362797056, 364500000, 368640000, 369056250, 373248000, 375000000, 377487360, 377913600,
2562 379687500, 382205952, 382637520, 384000000, 388800000, 390625000, 393216000, 393660000,
2563 398131200, 398580750, 400000000, 402653184, 403107840, 405000000, 408146688, 409600000,
2564 410062500, 414720000, 419430400, 419904000, 421875000, 424673280, 425152800, 429981696,
2565 432000000, 437400000, 439453125, 442368000, 442867500, 447897600, 450000000, 452984832,
2566 453496320, 455625000, 460800000, 466560000, 468750000, 471859200, 472392000, 474609375,
2567 477757440, 478296900, 480000000, 483729408, 486000000, 488281250, 491520000, 492075000,
2568 497664000, 500000000, 503316480, 503884800, 506250000, 509607936, 510183360, 512000000,
2569 512578125, 518400000, 524288000, 524880000, 527343750, 530841600, 531441000, 536870912,
2570 537477120, 540000000, 544195584, 546750000, 552960000, 553584375, 559872000, 562500000,
2571 566231040, 566870400, 569531250, 573308928, 576000000, 583200000, 585937500, 589824000,
2572 590490000, 597196800, 597871125, 600000000, 603979776, 604661760, 607500000, 612220032,
2573 614400000, 615093750, 622080000, 625000000, 629145600, 629856000, 632812500, 637009920,
2574 637729200, 640000000, 644972544, 648000000, 655360000, 656100000, 663552000, 664301250,
2575 671088640, 671846400, 675000000, 679477248, 680244480, 683437500, 691200000, 699840000,
2576 703125000, 707788800, 708588000, 716636160, 720000000, 725594112, 729000000, 732421875,
2577 737280000, 738112500, 746496000, 750000000, 754974720, 755827200, 759375000, 764411904,
2578 765275040, 768000000, 777600000, 781250000, 786432000, 787320000, 791015625, 796262400,
2579 797161500, 800000000, 805306368, 806215680, 810000000, 816293376, 819200000, 820125000,
2580 829440000, 838860800, 839808000, 843750000, 849346560, 850305600, 854296875, 859963392,
2581 864000000, 874800000, 878906250, 884736000, 885735000, 895795200, 900000000, 905969664,
2582 906992640, 911250000, 921600000, 922640625, 933120000, 937500000, 943718400, 944784000,
2583 949218750, 955514880, 956593800, 960000000, 967458816, 972000000, 976562500, 983040000,
2584 984150000, 995328000, 996451875, 1000000000, 1006632960, 1007769600, 1012500000,
2585 1019215872, 1020366720, 1024000000, 1025156250, 1036800000, 1048576000, 1049760000,
2586 1054687500, 1061683200, 1062882000, 1073741824, 1074954240, 1080000000, 1088391168,
2587 1093500000, 1105920000, 1107168750, 1119744000, 1125000000, 1132462080, 1133740800,
2588 1139062500, 1146617856, 1152000000, 1166400000, 1171875000, 1179648000, 1180980000,
2589 1194393600, 1195742250, 1200000000, 1207959552, 1209323520, 1215000000, 1220703125,
2590 1224440064, 1228800000, 1230187500, 1244160000, 1250000000, 1258291200, 1259712000,
2591 1265625000, 1274019840, 1275458400, 1280000000, 1289945088, 1296000000, 1310720000,
2592 1312200000, 1318359375, 1327104000, 1328602500, 1342177280, 1343692800, 1350000000,
2593 1358954496, 1360488960, 1366875000, 1382400000, 1399680000, 1406250000, 1415577600,
2594 1417176000, 1423828125, 1433272320, 1440000000, 1451188224, 1458000000, 1464843750,
2595 1474560000, 1476225000, 1492992000, 1500000000, 1509949440, 1511654400, 1518750000,
2596 1528823808, 1530550080, 1536000000, 1537734375, 1555200000, 1562500000, 1572864000,
2597 1574640000, 1582031250, 1592524800, 1594323000, 1600000000, 1610612736, 1612431360,
2598 1620000000, 1632586752, 1638400000, 1640250000, 1658880000, 1660753125, 1677721600,
2599 1679616000, 1687500000, 1698693120, 1700611200, 1708593750, 1719926784, 1728000000,
2600 1749600000, 1757812500, 1769472000, 1771470000, 1791590400, 1800000000, 1811939328,
2601 1813985280, 1822500000, 1843200000, 1845281250, 1866240000, 1875000000, 1887436800,
2602 1889568000, 1898437500, 1911029760, 1913187600, 1920000000, 1934917632, 1944000000,
2603 1953125000, 1966080000, 1968300000, 1990656000, 1992903750, 2000000000, 2013265920,
2604 2015539200, 2025000000, 2038431744, 2040733440, 2048000000, 2050312500, 2073600000,
2605 2097152000, 2099520000, 2109375000, 2123366400, 2125764000
2606 };
2607
2608 }
2609
2610 int cv::getOptimalDFTSize( int size0 )
2611 {
2612     int a = 0, b = sizeof(optimalDFTSizeTab)/sizeof(optimalDFTSizeTab[0]) - 1;
2613     if( (unsigned)size0 >= (unsigned)optimalDFTSizeTab[b] )
2614         return -1;
2615
2616     while( a < b )
2617     {
2618         int c = (a + b) >> 1;
2619         if( size0 <= optimalDFTSizeTab[c] )
2620             b = c;
2621         else
2622             a = c+1;
2623     }
2624
2625     return optimalDFTSizeTab[b];
2626 }
2627
2628 CV_IMPL void
2629 cvDFT( const CvArr* srcarr, CvArr* dstarr, int flags, int nonzero_rows )
2630 {
2631     cv::Mat src = cv::cvarrToMat(srcarr), dst0 = cv::cvarrToMat(dstarr), dst = dst0;
2632     int _flags = ((flags & CV_DXT_INVERSE) ? cv::DFT_INVERSE : 0) |
2633         ((flags & CV_DXT_SCALE) ? cv::DFT_SCALE : 0) |
2634         ((flags & CV_DXT_ROWS) ? cv::DFT_ROWS : 0);
2635
2636     CV_Assert( src.size == dst.size );
2637
2638     if( src.type() != dst.type() )
2639     {
2640         if( dst.channels() == 2 )
2641             _flags |= cv::DFT_COMPLEX_OUTPUT;
2642         else
2643             _flags |= cv::DFT_REAL_OUTPUT;
2644     }
2645
2646     cv::dft( src, dst, _flags, nonzero_rows );
2647     CV_Assert( dst.data == dst0.data ); // otherwise it means that the destination size or type was incorrect
2648 }
2649
2650
2651 CV_IMPL void
2652 cvMulSpectrums( const CvArr* srcAarr, const CvArr* srcBarr,
2653                 CvArr* dstarr, int flags )
2654 {
2655     cv::Mat srcA = cv::cvarrToMat(srcAarr),
2656         srcB = cv::cvarrToMat(srcBarr),
2657         dst = cv::cvarrToMat(dstarr);
2658     CV_Assert( srcA.size == dst.size && srcA.type() == dst.type() );
2659
2660     cv::mulSpectrums(srcA, srcB, dst,
2661         (flags & CV_DXT_ROWS) ? cv::DFT_ROWS : 0,
2662         (flags & CV_DXT_MUL_CONJ) != 0 );
2663 }
2664
2665
2666 CV_IMPL void
2667 cvDCT( const CvArr* srcarr, CvArr* dstarr, int flags )
2668 {
2669     cv::Mat src = cv::cvarrToMat(srcarr), dst = cv::cvarrToMat(dstarr);
2670     CV_Assert( src.size == dst.size && src.type() == dst.type() );
2671     int _flags = ((flags & CV_DXT_INVERSE) ? cv::DCT_INVERSE : 0) |
2672             ((flags & CV_DXT_ROWS) ? cv::DCT_ROWS : 0);
2673     cv::dct( src, dst, _flags );
2674 }
2675
2676
2677 CV_IMPL int
2678 cvGetOptimalDFTSize( int size0 )
2679 {
2680     return cv::getOptimalDFTSize(size0);
2681 }
2682
2683 /* End of file. */