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