图像的正交变换----傅立叶变换

/*************************************************************************

*

* 函数名称：

*   FFT()

*

* 参数:

*   complex<double> * TD	- 指向时域数组的指针

*   complex<double> * FD	- 指向频域数组的指针

*   r						－2的幂数，即迭代次数

*

* 返回值:

*   无。

*

* 说明:

*   该函数用来实现快速付立叶变换。

*

************************************************************************/

void FFT(std::complex<double>*TD,std::complex<double>*FD,int r);

/*************************************************************************

*

* 函数名称：

*   IFFT()

*

* 参数:

*   complex<double> * FD	- 指向频域值的指针

*   complex<double> * TD	- 指向时域值的指针

*   r						－2的幂数

*

* 返回值:

*   无。

*

* 说明:

*   该函数用来实现快速付立叶逆变换。

*

************************************************************************/

void IFFT(std::complex<double>* FD, std::complex<double>* TD,int r);

/*************************************************************************

*

* 函数名称：

*   Fourier()

*

* 参数:

*   complex* TD		- 输入的时域序列

*	 LONG lWidth		- 图象宽度

*	 LONG lHeight		- 图象高度

*	 complex* FD		- 输出的频域序列

*

* 返回值:

*   BOOL               - 成功返回TRUE，否则返回FALSE。

*

* 说明:

*   该函数进行二维快速付立叶变换。

*

************************************************************************/

BOOL Fourier(std::complex<double>* TD, LONG width, LONG height, std::complex<double>* FD);

/*************************************************************************

*

* 函数名称：

*   IFourier()

*

* 参数:

*   LPBYTE TD			- 返回的空域数据

*   LONG lWidth		- 空域图象的宽度

*	 LONG lHeight		- 空域图象的高度

*	 complex* FD		- 输入的频域数据

*

* 返回值:

*   BOOL               - 成功返回TRUE，否则返回FALSE。

*

* 说明:

*   该函数进行二维快速付立叶逆变换。

*

************************************************************************/

BOOL IFourier(LPBYTE TD, LONG width, LONG height, std::complex<double>* FD);

/*************************************************************************

*

* 函数名称：

*   BmpFourier()

*

* 参数:

*   BYTE *bmp --------待处理的图像

LONG width,LONG height-------图像的宽度和高度

*

* 返回值:

*   BOOL               - 成功返回TRUE，否则返回FALSE。

*

* 说明:

*   该函数用来对图像进行付立叶变换。

*

************************************************************************/

BOOL BmpFourier(BYTE* bmp,LONG width,LONG height);

voidMyProcess::FFT(std::complex<double>*TD,std::complex<double>*FD,int r)

{

// 循环变量

LONG	i;

LONG	j;

LONG	k;

// 中间变量

int		p;

// 角度

double	angle;

complex<double>*W,*X1,*X2,*X;

// 计算付立叶变换点数

LONG N =1<< r;

// 分配运算所需存储器

W  =new complex<double>[N /2];

X1 =new complex<double>[N];

X2 =new complex<double>[N];

// 计算加权系数

for(i =0; i < N /2; i++)

{

angle =-i * PI *2/ N;

W[i]= complex<double>(cos(angle), sin(angle));

}

// 将时域点写入X1

memcpy(X1, TD,sizeof(complex<double>)* N);

// 采用蝶形算法进行快速付立叶变换

for(k =0; k < r; k++)

{

for(j =0; j <1<< k; j++)

{

for(i =0; i <1<<(r - k -1); i++)

{

p = j *(1<<(r - k));

X2[i + p]= X1[i + p]+ X1[i + p +(int)(1<<(r - k -1))];

X2[i + p +(int)(1<<(r - k -1))]=(X1[i + p]- X1[i + p +(int)(1<<(r - k -1))])* W[i *(1<<k)];

}

}

X  = X1;

X1 = X2;

X2 = X;

}

// 重新排序

for(j =0; j < N; j++)

{

p =0;

for(i =0; i < r; i++)

{

if(j&(1<<i))

{

p+=1<<(r - i -1);

}

}

FD[j]= X1[p];

}

// 释放内存

delete W;

delete X1;

delete X2;

}

voidMyProcess::IFFT(complex<double>* FD, complex<double>* TD,int r)

{

// 循环变量

int		i;

complex<double>*X;

// 计算付立叶变换点数

LONG N =1<<r;

// 分配运算所需存储器

X =new complex<double>[N];

// 将频域点写入X

memcpy(X, FD,sizeof(complex<double>)* N);

// 求共轭

for(i =0; i < N; i++)

{

X[i]= complex<double>(X[i].real(),-X[i].imag());

}

// 调用快速付立叶变换

FFT(X, TD, r);

// 求时域点的共轭

for(i =0; i < N; i++)

{

TD[i]= complex<double>(TD[i].real()/ N,-TD[i].imag()/ N);

}

// 释放内存

delete X;

}

BOOL MyProcess::Fourier(std::complex<double>* TD, LONG width, LONG height, std::complex<double>* FD)

{

// 循环变量

LONG	i;

LONG	j;

LONG    k;

// 进行付立叶变换的宽度和高度（2的整数次方）

LONG w =1;

LONG h =1;

int wp =0;

int hp =0;

// 计算进行付立叶变换的宽度和高度（2的整数次方）

while(w < width/3)

{

w *=2;

wp++;

}

while(h < height)

{

h *=2;

hp++;

}

// 分配内存

complex<double>*TempT,*TempF;

TempT=new complex<double>[h];

TempF=new complex<double>[h];

// 对y方向进行快速付立叶变换

for(i =0; i < w *3; i++)

{

// 抽取数据

for(j =0; j < h; j++)

TempT[j]= TD[j * w *3+ i];

// 一维快速傅立叶变换

FFT(TempT,TempF, hp);

// 保存变换结果

for(j =0; j < h; j++)

TD[j * w *3+ i]=TempF[j];

}

// 释放内存

deleteTempT;

deleteTempF;

// 分配内存

TempT=new complex<double>[w];

TempF=new complex<double>[w];

// 对x方向进行快速付立叶变换

for(i =0; i < h; i++)

{

for(k =0; k <3; k++)

{

// 抽取数据

for(j =0; j < w; j++)

TempT[j]= TD[i * w *3+ j *3+ k];

// 一维快速傅立叶变换

FFT(TempT,TempF, wp);

// 保存变换结果

for(j =0; j < w; j++)

FD[i * w *3+ j *3+ k]=TempF[j];

}

}

// 释放内存

deleteTempT;

deleteTempF;

return TRUE;

}

BOOL MyProcess::IFourier(BYTE *TD, LONG width, LONG height, std::complex<double>* FD)

{

// 循环变量

LONG	i;

LONG	j;

LONG    k;

// 进行付立叶变换的宽度和高度（2的整数次方）

LONG w =1;

LONG h =1;

int wp =0;

int hp =0;

// 计算进行付立叶变换的宽度和高度（2的整数次方）

while(w < width/3)

{

w *=2;

wp++;

}

while(h < height)

{

h *=2;

hp++;

}

// 计算图像每行的字节数

// 分配内存

complex<double>*TempT,*TempF;

TempT=new complex<double>[w];

TempF=new complex<double>[w];

// 对x方向进行快速付立叶变换

for(i =0; i < h; i++)

{

for(k =0; k <3; k++)

{

// 抽取数据

for(j =0; j < w; j++)

TempF[j]= FD[i * w *3+ j *3+ k];

// 一维快速傅立叶变换

IFFT(TempF,TempT, wp);

// 保存变换结果

for(j =0; j < w; j++)

FD[i * w *3+ j *3+ k]=TempT[j];

}

}

// 释放内存

deleteTempT;

deleteTempF;

TempT=new complex<double>[h];

TempF=new complex<double>[h];

// 对y方向进行快速付立叶变换

for(i =0; i < w *3; i++)

{

// 抽取数据

for(j =0; j < h; j++)

TempF[j]= FD[j * w *3+ i];

// 一维快速傅立叶变换

IFFT(TempF,TempT, hp);

// 保存变换结果

for(j =0; j < h; j++)

FD[j * w *3+ i]=TempT[j];

}

// 释放内存

deleteTempT;

deleteTempF;

for(i =0; i < h; i++)

{

for(j =0; j < w *3; j++)

{

if(i < height && j < width)

*(TD + i * width + j)= FD[i * w *3+ j].real();

}

}

return TRUE;

}

BOOL MyProcess::BmpFourier(BYTE* bmp,LONG width,LONG height)

{

LONG i,j;//循环变量

// 中间变量

double	dTemp;

LONG TI,TJ;

// 进行付立叶变换的宽度和高度（2的整数次方）

LONG w =1;

LONG h =1;

int wp =0;

int hp =0;

// 计算进行付立叶变换的宽度和高度（2的整数次方）

while(w < width/3)

{

w *=2;

wp++;

}

while(h < height)

{

h *=2;

hp++;

}// 分配内存

complex<double>*FD,*TD,*TempD;

FD =new complex<double>[w * h *3];

TD =new complex<double>[w * h *3];

TempD=new complex<double>[w * h *3];

// 行

for(i =0; i < h; i++)

{

// 列

for(j =0; j <3* w; j++)

{

if(i < height && j < width)

{

// 获取时域数值

unsignedcharValue= bmp[i*width+j];

// 时域赋值

TD[w * i *3+ j]= complex<double>(Value,0.0f);

}

else

{

// 否则补零

TD[w * i *3+ j]= complex<double>(0.0f,0.0f);

}

}

}

// 进行频谱分析

if(Fourier(TD, width, height, FD)== FALSE)

return FALSE;

// 释放内存

delete[]TD;

// 将原点放置于图像中心位置

for(i =0; i < h; i++)

{

for(j =0; j <3* w; j++)

{

if(i < h /2)

TI = i + h /2;

else

TI = i - h /2;

if(j < w *3/2)

TJ = j +3* w /2;

else

TJ = j -3* w /2;

// 保存转换后的频谱图像

TempD[i * w *3+ j]= FD[TI * w *3+ TJ];

}

}

// 行

for(i =(int)(h - height)/2; i <(int)(h + height)/2; i++)

{

// 列

for(j =(int)(w *3- width)/2; j <(int)(w *3+ width)/2; j +=3)

{

// 计算频谱

dTemp = sqrt(TempD[w *3* i + j].real()*TempD[w *3* i + j].real()+

TempD[w *3* i + j].imag()*TempD[w *3* i + j].imag())/100;

// 判断是否超过255

if(dTemp >255)

{

// 对于超过的，直接设置为255

dTemp =255;

}

// 限制为原图大小范围

TI = i -(h - height)/2;

TJ = j /3-(w - width/3)/2;

// 对应象素指针

LONG p = width * TI + TJ *3;

// 更新源图像

bmp[p]=(BYTE)(dTemp);

bmp[p+1]=(BYTE)(dTemp);

bmp[p+2]=(BYTE)(dTemp);

}

}

// 删除临时变量

delete[]FD;

delete[]TempD;

return TRUE;

}