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

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