cvdxt.cpp.svn-base

来自「非结构化路识别」· SVN-BASE 代码 · 共 1,339 行 · 第 1/5 页

SVN-BASE
1,339
字号
/*M///////////////////////////////////////////////////////////////////////////////////////
//
//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
//
//  By downloading, copying, installing or using the software you agree to this license.
//  If you do not agree to this license, do not download, install,
//  copy or use the software.
//
//
//                        Intel License Agreement
//                For Open Source Computer Vision Library
//
// Copyright (C) 2000, Intel Corporation, all rights reserved.
// Third party copyrights are property of their respective owners.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//
//   * The name of Intel Corporation may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
//
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
// the use of this software, even if advised of the possibility of such damage.
//
//M*/

#include "_cv.h"

typedef struct CvPoint2D64f
{
    double x, y;
}
CvPoint2D64f;


/****************************************************************************************\
                               Discrete Fourier Transform
\****************************************************************************************/


/****************************************************************************************\
    Power-2 DFT in this file is a translation to C of Sorensen's split radix FFT.
    See comments below.
\****************************************************************************************/

static int icvlog2( int n )
{
    int m = 0;
    while(n>1)
        m++, n >>= 1;
    return m;
}


const double icvDxtTab[][2] =
{
{ 1.00000000000000000, 0.00000000000000000 },
{-1.00000000000000000, 0.00000000000000000 },
{ 0.00000000000000000, 1.00000000000000000 },
{ 0.70710678118654757, 0.70710678118654746 },
{ 0.92387953251128674, 0.38268343236508978 },
{ 0.98078528040323043, 0.19509032201612825 },
{ 0.99518472667219693, 0.09801714032956060 },
{ 0.99879545620517241, 0.04906767432741802 },
{ 0.99969881869620425, 0.02454122852291229 },
{ 0.99992470183914450, 0.01227153828571993 },
{ 0.99998117528260111, 0.00613588464915448 },
{ 0.99999529380957619, 0.00306795676296598 },
{ 0.99999882345170188, 0.00153398018628477 },
{ 0.99999970586288223, 0.00076699031874270 },
{ 0.99999992646571789, 0.00038349518757140 },
{ 0.99999998161642933, 0.00019174759731070 },
{ 0.99999999540410733, 0.00009587379909598 },
{ 0.99999999885102686, 0.00004793689960307 },
{ 0.99999999971275666, 0.00002396844980842 },
{ 0.99999999992818922, 0.00001198422490507 },
{ 0.99999999998204725, 0.00000599211245264 },
{ 0.99999999999551181, 0.00000299605622633 },
{ 0.99999999999887801, 0.00000149802811317 },
{ 0.99999999999971945, 0.00000074901405658 },
{ 0.99999999999992983, 0.00000037450702829 },
{ 0.99999999999998246, 0.00000018725351415 },
{ 0.99999999999999567, 0.00000009362675707 },
{ 0.99999999999999889, 0.00000004681337854 },
{ 0.99999999999999978, 0.00000002340668927 },
{ 0.99999999999999989, 0.00000001170334463 },
{ 1.00000000000000000, 0.00000000585167232 },
{ 1.00000000000000000, 0.00000000292583616 }
};

static unsigned char icvRevTable[] =
{
  0x00,0x80,0x40,0xc0,0x20,0xa0,0x60,0xe0,0x10,0x90,0x50,0xd0,0x30,0xb0,0x70,0xf0,
  0x08,0x88,0x48,0xc8,0x28,0xa8,0x68,0xe8,0x18,0x98,0x58,0xd8,0x38,0xb8,0x78,0xf8,
  0x04,0x84,0x44,0xc4,0x24,0xa4,0x64,0xe4,0x14,0x94,0x54,0xd4,0x34,0xb4,0x74,0xf4,
  0x0c,0x8c,0x4c,0xcc,0x2c,0xac,0x6c,0xec,0x1c,0x9c,0x5c,0xdc,0x3c,0xbc,0x7c,0xfc,
  0x02,0x82,0x42,0xc2,0x22,0xa2,0x62,0xe2,0x12,0x92,0x52,0xd2,0x32,0xb2,0x72,0xf2,
  0x0a,0x8a,0x4a,0xca,0x2a,0xaa,0x6a,0xea,0x1a,0x9a,0x5a,0xda,0x3a,0xba,0x7a,0xfa,
  0x06,0x86,0x46,0xc6,0x26,0xa6,0x66,0xe6,0x16,0x96,0x56,0xd6,0x36,0xb6,0x76,0xf6,
  0x0e,0x8e,0x4e,0xce,0x2e,0xae,0x6e,0xee,0x1e,0x9e,0x5e,0xde,0x3e,0xbe,0x7e,0xfe,
  0x01,0x81,0x41,0xc1,0x21,0xa1,0x61,0xe1,0x11,0x91,0x51,0xd1,0x31,0xb1,0x71,0xf1,
  0x09,0x89,0x49,0xc9,0x29,0xa9,0x69,0xe9,0x19,0x99,0x59,0xd9,0x39,0xb9,0x79,0xf9,
  0x05,0x85,0x45,0xc5,0x25,0xa5,0x65,0xe5,0x15,0x95,0x55,0xd5,0x35,0xb5,0x75,0xf5,
  0x0d,0x8d,0x4d,0xcd,0x2d,0xad,0x6d,0xed,0x1d,0x9d,0x5d,0xdd,0x3d,0xbd,0x7d,0xfd,
  0x03,0x83,0x43,0xc3,0x23,0xa3,0x63,0xe3,0x13,0x93,0x53,0xd3,0x33,0xb3,0x73,0xf3,
  0x0b,0x8b,0x4b,0xcb,0x2b,0xab,0x6b,0xeb,0x1b,0x9b,0x5b,0xdb,0x3b,0xbb,0x7b,0xfb,
  0x07,0x87,0x47,0xc7,0x27,0xa7,0x67,0xe7,0x17,0x97,0x57,0xd7,0x37,0xb7,0x77,0xf7,
  0x0f,0x8f,0x4f,0xcf,0x2f,0xaf,0x6f,0xef,0x1f,0x9f,0x5f,0xdf,0x3f,0xbf,0x7f,0xff
};

#define icvBitRev(i,shift) \
   ((int)((((unsigned)icvRevTable[(i)&255] << 24)+ \
           ((unsigned)icvRevTable[((i)>> 8)&255] << 16)+ \
           ((unsigned)icvRevTable[((i)>>16)&255] <<  8)+ \
           ((unsigned)icvRevTable[((i)>>24)])) >> (shift)))


static CvStatus icvInitBitRevTab( int* itab, int m )
{
    int i, n = (1 << m) >> 2, s = 34 - m;
    for( i = 0; i < n; i++ )
        itab[i] = icvBitRev(i,s)*4;
    return CV_OK;
}

#define ICV_BITREV_FUNC( flavor, datatype )                         \
static CvStatus icvBitRev_##flavor( datatype* v, int m, int* itab ) \
{                                                                   \
    datatype t;                                                     \
    int  n = 1 << m;                                                \
    int  j, k, l = n>>2;                                            \
    int  s = 34 - m;                                                \
    datatype *tmp0 = v, *tmp1 = v+l, *tmp2 = v+l*2, *tmp3 = v+l*3;  \
                                                                    \
    if( n >= 16 )                                                   \
    {                                                               \
        for( j = 0; j < l; j += 4 )                                 \
        {                                                           \
            k = itab ? itab[j] : icvBitRev(j,s)*4;                  \
                                                                    \
            CV_SWAP(tmp0[j+2],tmp1[k],t);                           \
            CV_SWAP(tmp0[j+1],tmp2[k],t);                           \
            CV_SWAP(tmp0[j+3],tmp3[k],t);                           \
            CV_SWAP(tmp1[j+1],tmp2[k+2],t);                         \
            CV_SWAP(tmp1[j+3],tmp3[k+2],t);                         \
            CV_SWAP(tmp2[j+3],tmp3[k+1],t);                         \
            if( k > j )                                             \
            {                                                       \
                CV_SWAP(tmp0[j  ],tmp0[k  ],t);                     \
                CV_SWAP(tmp1[j+2],tmp1[k+2],t);                     \
                CV_SWAP(tmp2[j+1],tmp2[k+1],t);                     \
                CV_SWAP(tmp3[j+3],tmp3[k+3],t);                     \
            }                                                       \
        }                                                           \
    }                                                               \
    else if( n > 2 )                                                \
    {                                                               \
        CV_SWAP(tmp0[1],tmp0[n>>1],t);                              \
        if( n == 8 )                                                \
            CV_SWAP(tmp0[3],tmp0[6],t);                             \
    }                                                               \
    return CV_OK;                                                   \
}


ICV_BITREV_FUNC( 32fc, CvPoint2D32f )
ICV_BITREV_FUNC( 64fc, CvPoint2D64f )

/*
CC=================================================================CC
CC                                                                 CC
CC  Subroutine CFFTSR(X,Y,M):                                      CC
CC      An in-place, split-radix complex FFT program               CC
CC      Decimation-in-frequency, cos/sin in second loop            CC
CC      and is computed recursively                                CC
CC      The program is based on Tran ASSP Feb 1986 pp152-156       CC
CC                                                                 CC
CC  Input/output                                                   CC
CC      X    Array of real part of input/output (length >= n)      CC
CC      Y    Array of imaginary part of input/output (length >= n) CC
CC      M    Transform length is n=2**M                            CC
CC                                                                 CC
CC  Calls:                                                         CC
CC      CSTAGE,CBITREV                                             CC
CC                                                                 CC
CC  Author:                                                        CC
CC      H.V. Sorensen,   University of Pennsylvania,  Dec. 1984    CC
CC                       Arpa address: hvs@ee.upenn.edu            CC
CC  Modified:                                                      CC
CC      H.V. Sorensen,   University of Pennsylvania,  Jul. 1987    CC
CC                                                                 CC
CC  Reference:                                                     CC
CC      Sorensen, Heideman, Burrus :"On computing the split-radix  CC
CC      FFT", IEEE Tran. ASSP, Vol. ASSP-34, No. 1, pp. 152-156    CC
CC      Feb. 1986                                                  CC
CC      Mitra&Kaiser: "Digital Signal Processing Handbook, Chap.   CC
CC      8, page 491-610, John Wiley&Sons, 1993                     CC
CC                                                                 CC
CC      This program may be used and distributed freely as         CC
CC      as long as this header is included                         CC
CC                                                                 CC
CC=================================================================CC
*/
#define X1(i)  v[(i)].x
#define Y1(i)  v[(i)].y
#define X2(i)  v[(i)+n4].x
#define Y2(i)  v[(i)+n4].y
#define X3(i)  v3[(i)].x
#define Y3(i)  v3[(i)].y
#define X4(i)  v3[(i)+n4].x
#define Y4(i)  v3[(i)+n4].y
#define CV_SIN_45 0.70710678118654752440084436210485

static CvStatus icvFFT_fwd_32fc( CvPoint2D32f* v, int m, int* itab )
{
    int k, n = 1 << m, n2 = n;
    int is, id, i1, i2;

    /*C-----L shaped butterflies------------------------------------------C*/
    for( k = m - 1; k > 0; k--, n2 >>= 1 )
    {
        int n4 = n2 >> 2;
        int n8 = n4 >> 1;
        CvPoint2D32f* v3 = v + 2*n4;

        /* CSTAGE */
        /*C-------Zero butterfly----------------------------------------------C*/

        for( is = 0, id = 2*n2; is < n; is = 2*id - n2, id *= 4 )
        {
            for( i1 = is; i1 < n; i1 += id )
            {
                double T1, T2;
            
                T1     = X1(i1) - X3(i1);
                X1(i1) = X1(i1) + X3(i1);
                T2     = Y2(i1) - Y4(i1);
                Y2(i1) = Y2(i1) + Y4(i1);
                X3(i1) = (float)(T1 + T2);
                T2     = T1 - T2;
                T1     = X2(i1) - X4(i1);
                X2(i1) = X2(i1) + X4(i1);
                X4(i1) = (float)T2;
                T2     = Y1(i1) - Y3(i1);
                Y1(i1) = Y1(i1) + Y3(i1);
                Y3(i1) = (float)(T2 - T1);
                Y4(i1) = (float)(T2 + T1);
            }
        }

        if( n4 <= 1 )
            continue;

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?