www.pudn.com > hongmoyuan.rar > GlobalApi.cpp
#include "stdafx.h"
#include "GlobalApi.h"
#include "Cdib.H"
#include
#include
#include
/*************************************************************************
*
* \函数名称：
* MakeGauss()
*
* \输入参数:
* double sigma - 高斯函数的标准差
* double **pdKernel - 指向高斯数据数组的指针
* int *pnWindowSize - 数据的长度
*
* \返回值:
* 无
*
* \说明:
* 这个函数可以生成一个一维的高斯函数的数字数据，理论上高斯数据的长度应
* 该是无限长的，但是为了计算的简单和速度，实际的高斯数据只能是有限长的
* pnWindowSize就是数据长度
*
**************************************************************************/
void MakeGauss(double sigma, double **pdKernel, int *pnWindowSize)
{
// 循环控制变量
int i ;
// 数组的中心点
int nCenter;
// 数组的某一点到中心点的距离
double dDis ;
double PI = 3.14159;
// 中间变量
double dValue;
double dSum ;
dSum = 0 ;
// 数组长度，根据概率论的知识，选取[-3*sigma, 3*sigma]以内的数据。
// 这些数据会覆盖绝大部分的滤波系数
*pnWindowSize = 1 + 2 * ceil(3 * sigma);
// 中心
nCenter = (*pnWindowSize) / 2;
// 分配内存
*pdKernel = new double[*pnWindowSize] ;
for(i=0; i< (*pnWindowSize); i++)
{
dDis = (double)(i - nCenter);
dValue = exp(-(1/2)*dDis*dDis/(sigma*sigma)) / (sqrt(2 * PI) * sigma );
(*pdKernel)[i] = dValue ;
dSum += dValue;
}
// 归一化
for(i=0; i<(*pnWindowSize) ; i++)
{
(*pdKernel)[i] /= dSum;
}
}
/*************************************************************************
*
* \函数名称：
* GaussianSmooth()
*
* \输入参数:
* unsigned char * pUnchImg - 指向图象数据的指针
* int nWidth - 图象数据宽度
* int nHeight - 图象数据高度
* double dSigma - 高斯函数的标准差
* unsigned char * pUnchSmthdImg - 指向经过平滑之后的图象数据
*
* \返回值:
* 无
*
* \说明:
* 为了抑止噪声，采用高斯滤波对图象进行滤波，滤波先对x方向进行，然后对
* y方向进行。
*
*************************************************************************
*/
void GaussianSmooth(unsigned char *pUnchImg, int nWidth, int nHeight,
double sigma, unsigned char * pUnchSmthdImg)
{
// 循环控制变量
int y;
int x;
int i;
// 高斯滤波器的数组长度
int nWindowSize;
// 窗口长度的1/2
int nHalfLen;
// 一维高斯数据滤波器
double *pdKernel ;
// 高斯系数与图象数据的点乘
double dDotMul ;
// 高斯滤波系数的总和
double dWeightSum ;
// 中间变量
double * pdTmp ;
// 分配内存
pdTmp = new double[nWidth*nHeight];
// 产生一维高斯数据滤波器
// MakeGauss(sigma, &dKernel, &nWindowSize);
MakeGauss(sigma, &pdKernel, &nWindowSize) ;
// MakeGauss返回窗口的长度，利用此变量计算窗口的半长
nHalfLen = nWindowSize / 2;
// x方向进行滤波
for(y=0; y= 0 && (i+x) < nWidth )
{
dDotMul += (double)pUnchImg[y*nWidth + (i+x)] * pdKernel[nHalfLen+i];
dWeightSum += pdKernel[nHalfLen+i];
}
}
pdTmp[y*nWidth + x] = dDotMul/dWeightSum ;
}
}
// y方向进行滤波
for(x=0; x= 0 && (i+y) < nHeight )
{
dDotMul += (double)pdTmp[(y+i)*nWidth + x] * pdKernel[nHalfLen+i];
dWeightSum += pdKernel[nHalfLen+i];
}
}
pUnchSmthdImg[y*nWidth + x] = (unsigned char)(int)dDotMul/dWeightSum ;
}
}
// 释放内存
delete []pdKernel;
pdKernel = NULL ;
delete []pdTmp;
pdTmp = NULL;
}
/*************************************************************************
* 函数名称：
* DirGrad()
* 输入参数:
* unsigned char *pUnchSmthdImg - 经过高斯滤波后的图象
* int nWidht - 图象宽度
* int nHeight - 图象高度
* int *pnGradX - x方向的方向导数
* int *pnGradY - y方向的方向导数
* 说明:
* 这个函数计算方向倒数，采用的微分算子是(-1 0 1) 和 (-1 0 1)'(转置)
* 计算的时候对边界象素采用了特殊处理
**************************************************************************/
void DirGrad(unsigned char *pUnchSmthdImg, int nWidth, int nHeight,
int *pnGradX , int *pnGradY)
{
// 循环控制变量
int y ;
int x ;
// 计算x方向的方向导数，在边界出进行了处理，防止要访问的象素出界
for(y=0; y abs(gx))
{
// 计算插值的比例
weight = fabs(gx)/fabs(gy);
g2 = pnMag[nPos-nWidth] ;
g4 = pnMag[nPos+nWidth] ;
// 如果x，y两个方向的方向导数的符号相同
// C是当前象素，与g1-g4的位置关系为：
// g1 g2
// C
// g4 g3
if (gx*gy > 0)
{
g1 = pnMag[nPos-nWidth-1] ;
g3 = pnMag[nPos+nWidth+1] ;
}
// 如果x，y两个方向的方向导数的符号相反
// C是当前象素，与g1-g4的位置关系为：
// g2 g1
// C
// g3 g4
else
{
g1 = pnMag[nPos-nWidth+1] ;
g3 = pnMag[nPos+nWidth-1] ;
}
}
// 如果方向导数x分量比y分量大，说明导数的方向更加“趋向”于x分量
// 这个判断语句包含了x分量和y分量相等的情况
else
{
// 计算插值的比例
weight = fabs(gy)/fabs(gx);
g2 = pnMag[nPos+1] ;
g4 = pnMag[nPos-1] ;
// 如果x，y两个方向的方向导数的符号相同
// C是当前象素，与g1-g4的位置关系为：
// g3
// g4 C g2
// g1
if (gx*gy > 0)
{
g1 = pnMag[nPos+nWidth+1] ;
g3 = pnMag[nPos-nWidth-1] ;
}
// 如果x，y两个方向的方向导数的符号相反
// C是当前象素，与g1-g4的位置关系为：
// g1
// g4 C g2
// g3
else
{
g1 = pnMag[nPos-nWidth+1] ;
g3 = pnMag[nPos+nWidth-1] ;
}
}
// 下面利用g1-g4对梯度进行插值
{
dTmp1 = weight*g1 + (1-weight)*g2 ;
dTmp2 = weight*g3 + (1-weight)*g4 ;
// 当前象素的梯度是局部的最大值
// 该点可能是个边界点
if(dTmp>=dTmp1 && dTmp>=dTmp2)
{
pUnchRst[nPos] = 128 ;
}
else
{
// 不可能是边界点
pUnchRst[nPos] = 0 ;
}
}
} //else
} // for
}
}
/*************************************************************************
* 函数名称：
* TraceEdge()
* 输入参数:
* int x - 跟踪起点的x坐标
* int y - 跟踪起点的y坐标
* int nLowThd - 判断一个点是否为边界点的低阈值
* unsigned char *pUnchEdge - 记录边界点的缓冲区
* int *pnMag - 梯度幅度图
* int nWidth - 图象数据宽度
* 说明:
* 递归调用
* 从(x,y)坐标出发，进行边界点的跟踪，跟踪只考虑pUnchEdge中没有处理并且
* 可能是边界点的那些象素(=128),象素值为0表明该点不可能是边界点，象素值
* 为255表明该点已经被设置为边界点，不必再考虑
**************************************************************************/
void TraceEdge (int y, int x, int nLowThd, unsigned char *pUnchEdge, int *pnMag, int nWidth)
{
// 对8邻域象素进行查询
int xNb[8] = {1, 1, 0,-1,-1,-1, 0, 1} ;
int yNb[8] = {0, 1, 1, 1,0 ,-1,-1,-1} ;
int yy ;
int xx ;
int k ;
for(k=0; k<8; k++)
{
yy = y + yNb[k] ;
xx = x + xNb[k] ;
// 如果该象素为可能的边界点，又没有处理过
// 并且梯度大于阈值
if(pUnchEdge[yy*nWidth+xx] == 128 && pnMag[yy*nWidth+xx]>=nLowThd)
{
// 把该点设置成为边界点
pUnchEdge[yy*nWidth+xx] = 255 ;
// 以该点为中心进行跟踪
TraceEdge(yy, xx, nLowThd, pUnchEdge, pnMag, nWidth);
}
}
}
/*************************************************************************
* 函数名称：
* EstimateThreshold()
* 输入参数:
* int *pnMag - 梯度幅度图
* int nWidth - 图象数据宽度
* int nHeight - 图象数据高度
* int *pnThdHigh - 高阈值
* int *pnThdLow - 低阈值
* double dRatioLow - 低阈值和高阈值之间的比例
* double dRatioHigh - 高阈值占图象象素总数的比例
* unsigned char *pUnchEdge - 经过non-maximum处理后的数据
* 说明:
* 经过non-maximum处理后的数据pUnchEdge，统计pnMag的直方图，确定阈值。
* 本函数中只是统计pUnchEdge中可能为边界点的那些象素。然后利用直方图，
* 根据dRatioHigh设置高阈值，存储到pnThdHigh。利用dRationLow和高阈值，
* 设置低阈值，存储到*pnThdLow。dRatioHigh是一种比例：表明梯度小于
* *pnThdHigh的象素数目占象素总数目的比例。dRationLow表明*pnThdHigh
* 和*pnThdLow的比例，这个比例在canny算法的原文里，作者给出了一个区间。
*
**************************************************************************/
void EstimateThreshold(int *pnMag, int nWidth, int nHeight, int *pnThdHigh,int *pnThdLow,
unsigned char * pUnchEdge, double dRatioHigh, double dRationLow)
{
// 循环控制变量
int y;
int x;
int k;
// 该数组的大小和梯度值的范围有关，如果采用本程序的算法，那么梯度的范围不会超过pow(2,10)
int nHist[1024] ;
// 可能的边界数目
int nEdgeNb ;
// 最大梯度值
int nMaxMag ;
int nHighCount ;
nMaxMag = 0 ;
// 初始化
for(k=0; k<1024; k++)
{
nHist[k] = 0;
}
// 统计直方图，然后利用直方图计算阈值
for(y=0; y= nThdHigh))
{
// 设置该点为边界点
pUnchEdge[nPos] = 255;
TraceEdge(y, x, nThdLow, pUnchEdge, pnMag, nWidth);
}
}
}
// 那些还没有被设置为边界点的象素已经不可能成为边界点
for(y=0; y0 && a0 && b MaxCount)
{
MaxCount = Count[i*nWidth+j][k];
Gao_temp1 = i; // 圆心纵坐标
Gao_temp2 = j; // 圆心横坐标
Gao_temp3 = k; // 半径
}
}
}
}
// 释放内存
delete []Count;
Count = NULL;
}
/*************************************************************************
* 函数名称：
* Bresenham()
* 输入参数:
* int Gao_temp1 - 记录圆心纵坐标
* int Gao_temp2 - 记录圆心横坐标
* int Gao_temp3 - 记录半径
* unsigned char *pUnchEdge - 记录原图像的缓冲区
**************************************************************************/
void Bresenham(int Gao_temp1, int Gao_temp2,int Gao_temp3,int nWidth,int nHeight,unsigned char *HoughImage)
{
int Bresenham_x = 0;
int Bresenham_y = Gao_temp3;
int Delta_D;
int Delta_1;
int Delta_2;
int direction;
Delta_D = (Bresenham_x +1)*(Bresenham_x +1) + (Bresenham_y -1)*(Bresenham_y -1) - Gao_temp3*Gao_temp3;
while(Bresenham_y >= 0)
{
HoughImage[(Gao_temp1+Bresenham_y)*nWidth+(Gao_temp2+Bresenham_x)] = 255;
HoughImage[(Gao_temp1-Bresenham_y)*nWidth+(Gao_temp2+Bresenham_x)] = 255;
HoughImage[(Gao_temp1-Bresenham_y)*nWidth+(Gao_temp2-Bresenham_x)] = 255;
HoughImage[(Gao_temp1+Bresenham_y)*nWidth+(Gao_temp2-Bresenham_x)] = 255;
if(Delta_D<0)
{
Delta_1 = 2*(Delta_D+Bresenham_y)-1;
if(Delta_1<=0)
direction = 1;
else
direction = 2;
}
else if(Delta_D>0)
{
Delta_2 = 2*(Delta_D-Bresenham_x)-1;
if(Delta_2<=0)
direction = 2;
else
direction = 3;
}
else
{
direction = 2;
}
switch(direction)
{
case 1:
Bresenham_x ++;
Delta_D += 2*Bresenham_x + 1;
break;
case 2:
Bresenham_x ++;
Bresenham_y --;
Delta_D += 2*(Bresenham_x - Bresenham_y + 1);
break;
case 3:
Bresenham_y --;
Delta_D += ((-2)*Bresenham_y + 1);
break;
}
}
}
/*************************************************************************
* 函数名称：
* Gradient()
* 输入参数:
* int Gao_temp1 - 记录内圆圆心纵坐标
* int Gao_temp2 - 记录内圆圆心横坐标
* int Gao_temp3 - 记录内圆半径
* int Gao_temp4 - 记录外圆圆心纵坐标
* int Gao_temp5 - 记录外圆圆心横坐标
* int Gao_temp6 - 记录外圆半径
* unsigned char *pUnchEdge - 记录原图像的缓冲区
**************************************************************************/
void Gradient(int Gao_temp1,int Gao_temp2,int Gao_temp3,int nWidth,int nHeight,unsigned char *HoughImage,int &Gao_temp4,int &Gao_temp5,int &Gao_temp6)
{
//循环变量
int i = 0;
int j = 0;
int k = 0;
int Big_R = Gao_temp3 + 55;
double MaxCount = 0.0;
double (*grey)[120] = new double[nHeight*nWidth][120];
double (*grad)[119] = new double[nHeight*nWidth][119];
for(i=0; i= 0)
{
sum = sum + double(HoughImage[(i+Bresenham_y)*nWidth+(j+Bresenham_x)] + HoughImage[(i-Bresenham_y)*nWidth+(j+Bresenham_x)]
+ HoughImage[(i-Bresenham_y)*nWidth+(j-Bresenham_x)] + HoughImage[(i+Bresenham_y)*nWidth+(j-Bresenham_x)]);
count = count + 4;
if(Delta_D<0)
{
Delta_1 = 2*(Delta_D+Bresenham_y)-1;
if(Delta_1<=0)
direction = 1;
else
direction = 2;
}
else if(Delta_D>0)
{
Delta_2 = 2*(Delta_D-Bresenham_x)-1;
if(Delta_2<=0)
direction = 2;
else
direction = 3;
}
else
{
direction = 2;
}
switch(direction)
{
case 1:
Bresenham_x ++;
Delta_D += 2*Bresenham_x + 1;
break;
case 2:
Bresenham_x ++;
Bresenham_y --;
Delta_D += 2*(Bresenham_x - Bresenham_y + 1);
break;
case 3:
Bresenham_y --;
Delta_D += ((-2)*Bresenham_y + 1);
break;
}
}
grey[i*nWidth+j][k] = sum/count;
}
}
}
}
for(i=0; i= MaxCount)
{
MaxCount = fabs(grad[i*nWidth+j][k]);
Gao_temp4 = i; // 虹膜圆心纵坐标
Gao_temp5 = j; // 虹膜圆心横坐标
Gao_temp6 = k; // 虹膜半径
}
}
}
}
//释放内存
delete []grey;
grey = NULL;
delete []grad;
grad = NULL;
}
/*************************************************************************
* 函数名称：
* Normalization()
* 输入参数:
* int Gao_temp1 - 记录瞳孔圆心纵坐标
* int Gao_temp2 - 记录瞳孔圆心横坐标
* int Gao_temp3 - 记录瞳孔半径
* int Gao_temp4 - 记录虹膜圆心纵坐标
* int Gao_temp5 - 记录虹膜圆心横坐标
* int Gao_temp6 - 记录虹膜半径
* double x - 插值元素的x坐标
* double y - 插值元素的y坐标
************************************************************************/
void Normalization(unsigned char *NormalizeImage,unsigned char *HoughImage,
int Gao_temp1,int Gao_temp2,int Gao_temp3,int Gao_temp4,
int Gao_temp5,int Gao_temp6,int nWidth,int nHeight)
{
double alfa = 0.0;
double Theta ;
double angle_OIA = 0.0;
double angle_OAI = 0.0;
double angle_IOA = 0.0;
double OI = 0.0;
double IA = 0.0;
double r = 0.0;
double xp = 0.0;
double yp = 0.0;
double xi = 0.0;
double yi = 0.0;
double x = 0.0;
double y = 0.0;
for(Theta=0.0001; Theta<2*pi; Theta+=2*pi/1024)
{
alfa = atan(fabs((double)(Gao_temp1-Gao_temp4)/(double)(Gao_temp2-Gao_temp5)));
if(Theta2*pi/3)
{
if(Theta==alfa||Theta==alfa+pi)
{
xp = Gao_temp2+Gao_temp3*cos(Theta);
yp = Gao_temp1+Gao_temp3*sin(Theta);
xi = Gao_temp5+Gao_temp6*cos(Theta);
yi = Gao_temp4+Gao_temp6*sin(Theta);
}
else
{
angle_OIA = pi- Theta + alfa;
OI = sqrt((Gao_temp1-Gao_temp4)*(Gao_temp1-Gao_temp4)+(Gao_temp2-Gao_temp5)*(Gao_temp2-Gao_temp5));
angle_OAI = asin(OI*sin(angle_OIA)/Gao_temp6);
angle_IOA = pi - angle_OIA - angle_OAI;
IA = sqrt(OI*OI+Gao_temp6*Gao_temp6-2*OI*Gao_temp6*cos(angle_IOA));
xp = Gao_temp2 + Gao_temp3 * cos(Theta);
yp = Gao_temp1 + Gao_temp3 * sin(Theta);
xi = Gao_temp2 + IA * cos(Theta);
yi = Gao_temp1 + IA * sin(Theta);
}
for(r=0.0;r<1.0;r+=1.0/128)
{
x = (1-r)*xp + r*xi;
y = (1-r)*yp + r*yi;
NormalizeImage[int(128*r)*1024+int(Theta*1024/(2*pi))] = Interpolation(HoughImage, nWidth, nHeight, x, y);
}
}
}
}
/*************************************************************************
* 函数名称：
* Interpolation()
* 输入参数:
* unsigned char *InterpolaImage - 记录插值图像的缓冲区
* int nWidth - 图象数据宽度
* int nHeight - 图象数据高度
* float x - 插值元素的x坐标
* float y - 插值元素的y坐标
************************************************************************/
unsigned char Interpolation(unsigned char *InterpolaImage, int nWidth, int nHeight, double x, double y)
{
// 四个最临近象素的坐标(i1, j1), (i2, j1), (i1, j2), (i2, j2)
int i1, i2;
int j1, j2;
// 四个最临近象素值
char f1,f2,f3,f4;
// 二个插值中间值
char f12, f34;
// 计算四个最临近象素的坐标
i1 = (int) x;
i2 = i1 + 1;
j1 = (int) y;
j2 = j1 + 1;
// 根据不同情况分别处理
if( (x < 0) || (x > nWidth - 1) || (y < 0) || (y > nHeight - 1))
{
return 255;
}
else
{
if (x == nWidth - 1)
{
// 要计算的点在图像右边缘上
if (y == 0)
{
// 要计算的点正好是图像最右下角那一个象素，直接返回该点象素值
f1 = InterpolaImage[j1*nWidth+i1];
return f1;
}
else
{
// 在图像右边缘上且不是最后一点，直接一次插值即可
f1 = InterpolaImage[j1*nWidth+i1];
f3 = InterpolaImage[j1*nWidth+i2];
// 返回插值结果
return ((unsigned char)(f1 + (y -j1) * (f3 - f1)));
}
}
else if (y == 0)
{
// 要计算的点在图像下边缘上且不是最后一点，直接一次插值即可
f1 = InterpolaImage[j1*nWidth+i1];
f2 = InterpolaImage[j2*nWidth+i1];
// 返回插值结果
return ((unsigned char)(f1 + (x -i1) * (f2 - f1)));
}
else
{
// 计算四个最临近象素值
f1 = InterpolaImage[j1*nWidth+i1];
f2 = InterpolaImage[j2*nWidth+i1];
f3 = InterpolaImage[j1*nWidth+i2];
f4 = InterpolaImage[j2*nWidth+i2];
// 插值1
f12 = (unsigned char) (f1 + (x - i1) * (f2 - f1));
// 插值2
f34 = (unsigned char) (f3 + (x - i1) * (f4 - f3));
// 插值3
return ((unsigned char) (f12 + (y -j1) * (f34 - f12)));
}
}
}
/*************************************************************************
* 函数名称：
* InteEqualize()
* 参数:
* InteEqualizeImg - 指向图象数据的指针
* int nWidth - 图象数据宽度
* int nHeight - 图象数据高度
************************************************************************/
void InteEqualize(unsigned char *InteEqualizeImg, int nWidth,int nHeight)
{
// 循环变量
int i = 0;
int j = 0;
// 临时变量
int Temp = 0;
// 灰度映射表
BYTE bMap[256];
// 灰度映射表
int Count[256];
// 清零
for (i = 0; i < 256;i++)
{
Count[i] = 0;
}
// 计算各个灰度值的计数
for (i = 0; i < nHeight; i ++)
{
for (j = 0; j < nWidth; j ++)
{
Count[InteEqualizeImg[i*nWidth+j]]++;
}
}
// 计算灰度映射表
for (i = 0; i < 256; i++)
{
// 初始为0
Temp = 0;
for (j = 0; j <= i ;j++)
{
Temp += Count[j];
}
// 计算对应的新灰度值
bMap[i] = (BYTE) (Temp * 255 / nHeight / nWidth);
}
for(i = 0; i < nHeight; i++)
{
for(j = 0; j < nWidth; j++)
{
// 计算新的灰度值
InteEqualizeImg[i*nWidth+j] = bMap[InteEqualizeImg[i*nWidth+j]];
}
}
}
/*************************************************************************
* 函数名称：
* FFT()
* 输入参数:
* Cdib *FourierImg - 图像对象
*************************************************************************/
void FFT( unsigned char *FourierImg)
{
// 循环控制变量
int y;
int x;
int nWidth = 1024;
int nHeight= 128;
// 临时变量
double dTmpOne;
double dTmpTwo;
// 傅立叶变换竖直方向点数
int nTransHeight;
// 傅立叶变换水平方向点数
int nTransWidth;
// 计算进行傅立叶变换的点数 （2的整数次幂）
dTmpOne = log(nWidth)/log(2);
dTmpTwo = ceil(dTmpOne);
dTmpTwo = pow(2,dTmpTwo);
nTransWidth = (int) dTmpTwo;
// 计算进行傅立叶变换的点数 （2的整数次幂）
dTmpOne = log(nHeight)/log(2);
dTmpTwo = ceil(dTmpOne);
dTmpTwo = pow(2,dTmpTwo);
nTransHeight = (int) dTmpTwo;
// 计算图象数据存储每行需要的字节数
// BMP文件的每行数据存储是DWORD对齐的
int nSaveWidth;
nSaveWidth = ( (nWidth << 3) + 31)/32 * 4;
// 图象象素值
unsigned char unchValue;
// 指向时域数据的指针
complex * pCTData;
// 指向频域数据的指针
complex * pCFData;
// 分配内存
pCTData=new complex[nTransWidth * nTransHeight];
pCFData=new complex[nTransWidth * nTransHeight];
// 初始化
// 图象数据的宽和高不一定是2的整数次幂，所以pCTData
// 有一部分数据需要补0
for(y=0; y(0,0);
}
}
// 把图象数据传给pCTData
for(y=0; y(unchValue,0);
}
}
// 傅立叶正变换
DIBFFT_2D(pCTData, nWidth, nHeight, pCFData) ;
// 临时变量
double dTmp;
for(y=0; y * pCTData - 图像数据
* int nWidth - 数据宽度
* int nHeight - 数据高度
* complex * pCFData - 傅立叶变换后的结果
*************************************************************************/
VOID WINAPI DIBFFT_2D(complex * pCTData, int nWidth, int nHeight, complex * pCFData)
{
// 循环控制变量
int x;
int y;
// 临时变量
double dTmpOne;
double dTmpTwo;
// 进行傅立叶变换的宽度和高度，（2的整数次幂）
// 图像的宽度和高度不一定为2的整数次幂
int nTransWidth;
int nTransHeight;
// 计算进行傅立叶变换的宽度 （2的整数次幂）
dTmpOne = log(nWidth)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransWidth = (int) dTmpTwo ;
// 计算进行傅立叶变换的高度 （2的整数次幂）
dTmpOne = log(nHeight)/log(2);
dTmpTwo = ceil(dTmpOne) ;
dTmpTwo = pow(2,dTmpTwo) ;
nTransHeight = (int) dTmpTwo ;
// x，y（行列）方向上的迭代次数
int nXLev;
int nYLev;
// 计算x，y（行列）方向上的迭代次数
nXLev = (int) ( log(nTransWidth)/log(2) + 0.5 );
nYLev = (int) ( log(nTransHeight)/log(2) + 0.5 );
for(y = 0; y < nTransHeight; y++)
{
// x方向进行快速傅立叶变换
FFT_1D(&pCTData[nTransWidth * y], &pCFData[nTransWidth * y], nXLev);
}
// pCFData中目前存储了pCTData经过行变换的结果
// 为了直接利用FFT_1D，需要把pCFData的二维数据转置，再一次利用FFT_1D进行
// 傅立叶行变换（实际上相当于对列进行傅立叶变换）
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTData[nTransHeight * x + y] = pCFData[nTransWidth * y + x];
}
}
for(x = 0; x < nTransWidth; x++)
{
// 对x方向进行快速傅立叶变换，实际上相当于对原来的图象数据进行列方向的
// 傅立叶变换
FFT_1D(&pCTData[x * nTransHeight], &pCFData[x * nTransHeight], nYLev);
}
// pCFData中目前存储了pCTData经过二维傅立叶变换的结果，但是为了方便列方向
// 的傅立叶变换，对其进行了转置，现在把结果转置回来
for(y = 0; y < nTransHeight; y++)
{
for(x = 0; x < nTransWidth; x++)
{
pCTData[nTransWidth * y + x] = pCFData[nTransHeight * x + y];
}
}
memcpy(pCTData, pCFData, sizeof(complex) * nTransHeight * nTransWidth );
}
/*************************************************************************
* 函数名称：
* FFT_1D()
*
* 输入参数:
* complex * pCTData - 指向时域数据的指针，输入的需要变换的数据
* complex * pCFData - 指向频域数据的指针，输出的经过变换的数据
* nLevel －傅立叶变换蝶形算法的级数，2的幂数，
**************************************************************************/
VOID WINAPI FFT_1D(complex * pCTData, complex * pCFData, int nLevel)
{
// 循环控制变量
int i;
int j;
int k;
double PI = 3.1415926;
// 傅立叶变换点数
int nCount =0 ;
// 计算傅立叶变换点数
nCount =(int)pow(2,nLevel) ;
// 某一级的长度
int nBtFlyLen;
nBtFlyLen = 0 ;
// 变换系数的角度 ＝2 * PI * i / nCount
double dAngle;
complex *pCW ;
// 分配内存，存储傅立叶变化需要的系数表
pCW = new complex[nCount / 2];
// 计算傅立叶变换的系数
for(i = 0; i < nCount / 2; i++)
{
dAngle = -2 * PI * i / nCount;
pCW[i] = complex ( cos(dAngle), sin(dAngle) );
}
// 变换需要的工作空间
complex *pCWork1,*pCWork2;
// 分配工作空间
pCWork1 = new complex[nCount];
pCWork2 = new complex[nCount];
// 临时变量
complex *pCTmp;
// 初始化，写入数据
memcpy(pCWork1, pCTData, sizeof(complex) * nCount);
// 临时变量
int nInter;
nInter = 0;
// 蝶形算法进行快速傅立叶变换
for(k = 0; k < nLevel; k++)
{
for(j = 0; j < (int)pow(2,k); j++)
{
//计算长度
nBtFlyLen = (int)pow( 2,(nLevel-k) );
//倒序重排，加权计算
for(i = 0; i < nBtFlyLen/2; i++)
{
nInter = j * nBtFlyLen;
pCWork2[i + nInter] =
pCWork1[i + nInter] + pCWork1[i + nInter + nBtFlyLen / 2];
pCWork2[i + nInter + nBtFlyLen / 2] =
(pCWork1[i + nInter] - pCWork1[i + nInter + nBtFlyLen / 2])
* pCW[(int)(i * pow(2,k))];
}
}
// 交换 pCWork1和pCWork2的数据
pCTmp = pCWork1 ;
pCWork1 = pCWork2 ;
pCWork2 = pCTmp ;
}
// 重新排序
for(j = 0; j < nCount; j++)
{
nInter = 0;
for(i = 0; i < nLevel; i++)
{
if ( j&(1<