傅里叶变换与数字图像处理

标签：
it
傅里叶变换是将时域信号分解为不同频率的正弦和/余弦和的形式。

傅里叶变换是
数字图像处理技术的基础，其通过在时域和频域来回切换图像，对图像的信息特征进行提取和分析。

一维傅里叶变换及其反变换
单变量连续函数，f(x)的傅里叶变换F(u)定义为等式：
u=0, 1, 2,…，M
1
相瓯给定F®),通过溥里叶氏变换可以訣得K幼
J—tZ
这两个等式姐咸了傭里叶变换对9它fl指出在前一节中提到的董要尊实.即一t■函鞍可以从它的反变换中重新获得。

: 尸仏V)二[[/(砂)尹诚旳如刼
J 一8 J —OC
类愎地，反变換为：
/(X, v) = f I F(阳)/咖+诃
J—00 00
单变童离欝函b N・r M-l)的傅里叶变换由以下等式绪也
I M-1
.v=0
同样，给出F(u)，能用反DFT来获得原函数:
川门二工尸(心」皿卫
i/=0
x=Of 1? 2i■> MT
从瓯拉公式中得至花
""-cos^+ /sin^
律出:
A/-1
1 'I -1
=一V f ix\Qos(—2mtx}/ M + / sin(—2^?/A i
f ( A- /cos 2.77a / M -j sin 2 mtx / M)
其中，u=0, 1, 2,…，M 一1。

因此，我们看到傅里叶变换的每项［即对于每个u 值，F(u)的值由f(x)函数所有值的和组成。

f(x)的值则与各种频率的正弦值和余弦值相乘。

F(u)值的范围覆盖的域(u的值)称为频率域，因为u决定了变换的频率成分(x 也作用于频率，但它们相加，对每个u值有相同的贡献)。

F(u)的M项中的每一个被称为变换的频率分量。

使用术语频率域”和频率成分”与时间域”和时间成分”没有差别，如果x是一个时间变量，可以用它来表示f(x)的域和值。

二维DFT及其反变换
一维离散傅里叶变换及其反变换向二维扩展是简单明了的。

一个图像尺寸为MX N
的函数f(x，y)的离散傅里叶变换由以下等式给出：
像在一维中的情形一样，此表达式必须对u值(u=0, 1, 2,…，M-1)和v值(v=0 , 1, 2,…，N-1)计算。

同样，给出F(u,v)，可以通过反傅里叶变换获得，f(x , y),
由表达式给出：
M -1 N-1
»= 0 v = 0
其中，x=0, 1, 2,…，M-1 , y=0, 1, 2,…，N-1。

变量u和v是变换或频率变量，x和y 是空间或图像变量。

正如在一维中的情形那样，常量1/MN的位置并不重要，有时它在反变换之前。

其他时候，它被分为两个相等的常数1/根号MN，分别乘在
变换和反变换的式子前。

定义傅里叶谱、相角和频率谱：
—arc tan
并且其功率谱为：
P(?/,v)=|F(z/,v)|
2=X W>V)2十A",讦
其中，R(u,v)和I(u,v)分别是F(u,v)的实部和虚部。

通常在进行傅里叶变换之前用(-1)x+y乘以输入的图像函数。

由于指数的性质，很容易看出:
其中Z ［表示引文中的傅里叶变换。

这个等式说明f(x,y)(-1)x+y傅里叶变换的原点［即
F(0,0)］被设置在u二M /2和v二N/2上。

换句话说，用(-1)x+y乘以f(x,y)将F(u,v)原点变
换到频率坐标下的(M/2; N/2)，它是二维DFT设置的MX N区域的中心。

我们将此频率域的范围指定为频率矩形，它从u=0到u=M-1从v=0到v=N-1(u和v是整数)。

为了确保移动后的坐标为整数，要求M和N为偶数。

当在计算机中使用傅里叶变换时，总和的范围为u从1到M , v从1到N。

实际的变换中心将为u=(M /2)+1 和v=(N / 2)十1。

当(u,v)=(0,0)时的变换值为：
1 M-1N-1
MN铝仁
即f(x,y)的平均值。

换句话说，如果f(x,y)是一幅图像，在原点的傅里叶变换即等于图像的平均灰度级。

因为在原点处常常为零，F(0,0)有时称做频率谱的直流成分。

了解了傅里叶变化，下面看看为什么要在频率域研究图像增强。

1•可以利用频率成分和图像外表之间的对应关系。

一些在空间域表述困难的增强任务，在频率域中变得非常普通。

2•滤波在频率域更为直观，它可以解释空间域滤波的某些性质。

3•可以在频率域指定滤波器，做反变换，然后在空间域使用结果滤波器作为空间域滤波器的指导。

4•一旦通过频率域试验选择了空间滤波，通常实施都在空间域进行。

在冈萨雷斯版＜数字图像处理＞里面的解释就非常的形象：一个恰当的比喻是将傅里叶变换比作一个玻璃棱镜。

棱镜是可以将光分解为不同颜色的物理仪器，每个成分的颜色由波长(或频率)来决定。

傅里叶变换可以看作是数学上的棱镜，将函数基于频率分解为不同的成分。

当我们考虑光时，讨论它的光谱或频率谱。

同样，傅立叶变换使我们能通过频率成分来分析一个函数。

经过变换之后的图像在原点平移之前四角是低频， 最亮，平移之后中间部分是低频, 最亮，亮度大说明低频的能量大。

傅立叶变换在图像处理中有非常非常的作用。

因为不仅傅立叶分析涉及图像处理的
很多方面，傅立叶的改进算法， 比如离散余弦变换，gabor 与小波在图像处理中也有重要的分量 印象中，傅立叶变换在图像处理以下几个话题都有重要作用：
1•图像增强与图像去噪
绝大部分噪音都是图像的高频分量，通过低通滤波器来滤除高频
是图像的高频分量，可以通过添加高频分量来增强原始图像的边缘
;噪声；边缘也
2•图像分割之边缘检测
提取图像高频分量
3•图像特征提取：
形状特征：傅里叶描述子
纹理特征：直接通过傅里叶系数来计算纹理特征
其他特征：将提取的特征值进行傅里叶变换来使特征具有平移、伸缩、旋转不变性
4•图像压缩
可以直接通过傅里叶系数来压缩数据；常用的离散余弦变换是傅立叶变换的实变换；
傅立叶变换
傅里叶变换是将时域信号分解为不同频率的正弦信号或余弦函数叠加之和。

连续情况下要求原始信号在一个周期内满足绝对可积条件。

离散情况下，傅里叶变换一定存在。

冈萨雷斯版＜图像处理＞里面的解释非常形象：一个恰当的比喻是将傅里叶变换比作一个玻璃棱镜。

棱镜是可以将光分解为不同颜色的物理仪器，每个成分的颜色由波长（或频率）来决定。

傅里叶变换可以看作是数学上的棱镜，将函数基于频率分解为不同的成分。

当我们考虑光时，讨论它的光谱或频率谱。

同样,傅立叶变换使我们能通过频率成分来分析一个函数。

傅立叶变换有很多优良的性质。

比如线性，对称性（可以用在计算信号的傅里叶变
换里面）；
时移性：函数在时域中的时移，对应于其在频率域中附加产生的相移，而幅度频谱则保持不变；
频移性：函数在时域中乘以e A jwt，可以使整个频谱搬移w。

这个也叫调制定理，通讯里面信号的频分复用需要用到这个特性（将不同的信号调制到不同的频段上同时传输）；
卷积定理：时域卷积等于频域乘积；时域乘积等于频域卷积（附加一个系数）。

（图像处理里面这个是个重点）
信号在频率域的表现
在频域中，频率越大说明原始信号变化速度越快；频率越小说明原始信号越平缓。

当频率为0时，表示直流信号，没有变化。

因此，频率的大小反应了信号的变化快
慢。

高频分量解释信号的突变部分，而低频分量决定信号的整体形象。

在图像处理中，频域反应了图像在空域灰度变化剧烈程度，也就是图像灰度的变化速度，也就
是图像的梯度大小。

对图像而言，图像的边缘部分是突变部分，变化较快，因此反应在频域上是高频分量；图像的噪声大部分情况下是高频部分；图像平缓变化部分则为低频分量。

也就是说，傅立叶变换提供另外一个角度来观察图像，可以将图像从灰度分布转化到频率分布上来观察图像的特征。

书面一点说就是，傅里叶变换提供了一条从空域到频率自由转换的途径。

对图像处理而言，以下概念非常的重要：
图像高频分量：图像突变部分；在某些情况下指图像边缘信息，某些情况下指噪声, 更多是两者的混合；
低频分量：图像变化平缓的部分，也就是图像轮廓信息
高通滤波器：让图像使低频分量抑制，高频分量通过
低通滤波器：与高通相反，让图像使高频分量抑制，低频分量通过
带通滤波器：使图像在某一部分的频率信息通过，其他过低或过高都抑制
还有个带阻滤波器，是带通的反。

模板运算与卷积定理
在时域内做模板运算，实际上就是对图像进行卷积。

模板运算是图像处理一个很重要的处理过程，很多图像处理过程，比如增强/去噪（这两个分不清楚），边缘检测中普遍用到。

根据卷积定理，时域卷积等价与频域乘积。

因此，在时域内对图像做模板运算就等效于在频域内对图像做滤波处理。

比如说一个均值模板，其频域响应为一个低通滤波器；在时域内对图像作均值滤波就等效于在频域内对图像用均值模板的频域响应对图像的频域响应作一个低通滤波。

图像去噪
图像去噪就是压制图像的噪音部分。

因此，如果噪音是高频额，从频域的角度来看, 就是需要用一个低通滤波器对图像进行处理。

通过低通滤波器可以抑制图像的高频分量。

但是这种情况下常常会造成边缘信息的抑制。

常见的去噪模板有均值模板，高斯模板等。

这两种滤波器都是在局部区域抑制图像的高频分量，模糊图像边缘的
同时也抑制了噪声。

还有一种非线性滤波-中值滤波器。

中值滤波器对脉冲型噪声有很好的去掉。

因为脉冲点都是突变的点，排序以后输出中值，那么那些最大点和最小点就可以去掉了。

中值滤波对高斯噪音效果较差。

椒盐噪声：对于椒盐采用中值滤波可以很好的去除。

用均值也可以取得一定的效果，但是会引起边缘的模糊。

高斯白噪声：白噪音在整个频域的都有分布，好像比较困难。

冈萨雷斯版图像处理P185 :算术均值滤波器和几何均值滤波器(尤其是后者)更适
合于处理高斯或者均匀的随机噪声。

谐波均值滤波器更适合于处理脉冲噪声。

图像增强
有时候感觉图像增强与图像去噪是一对矛盾的过程，图像增强经常是需要增强图像的边缘，以获得更好的显示效果，这就需要增加图像的高频分量。

而图像去噪是为了消除图像的噪音，也就是需要抑制高频分量。

有时候这两个又是指类似的事情。

比如说，消除噪音的同时图像的显示效果显著的提升了，那么，这时候就是同样的意思了。

常见的图像增强方法有对比度拉伸，直方图均衡化，图像锐化等。

前面两个是在空域进行基于像素点的变换，后面一个是在频域处理。

我理解的锐化就是直接在图像上加上图像高通滤波后的分量，也就是图像的边缘效果。

对比度拉伸和直方图均衡化都是为了提高图像的对比度，也就是使图像看起来差异更明显一些，我想，经过这样的处理以后，图像也应该增强了图像的高频分量，使得图像的细节上差异更大。

同时也引入了一些噪音。

* /函数名称：
* LowPassFilter()
* /输入参数:
* LPSTR IpDIBBits - 指向需要滤波的图像像素指针
* int nWidth - 数据宽度
* int nHeight - 数据高度
* int nRadius -理想低通滤波的滤波半径
* /返回值:
* BOOL //成功返回TRUE，否则返回FALSE。

* /说明: * IpDIBBits 是指向需要滤波的图像像素指针。

*经过低通滤波的数据仍然存储在IpDIBBits 当中。

**************************************************************************/
BOOL LowPassFiIter(LPSTR IpDIBBits, LONG nWidth, LONG nHeight 」nt n Radius) {
un sig ned char* lpSrc; //
指向源图像的指针 int y ;
//循环控制变量 int x ; //循环控制变量
double dTmp One ; double dTmpTwo ;
int nTran sWidth ; int nTran sHeight;
double un chValue; complex<double> * pCTData ; // complex<double> * pCFData ; //
//计算进行傅立叶变换的点数—横向
dTmpO ne = log( nWidth)/log(2);
dTmpTwo = ceil(dTmp One);
dTmpTwo = pow(2,dTmpTwo); nTran
sWidth = (int) dTmpTwo;
//计算进行傅立叶变换的点数—纵向 dTmpO ne = log( nH eight)/log(2); dTmpTwo = ceil(dTmp One); dTmpTwo = pow(2,dTmpTwo); nTran sHeight = (int) dTmpTwo;
double dReal;
// 傅立叶变换的实部 double dImag; // 傅立叶变换的虚部
//低通滤波的半径不能超过频域的最大半径
if （n Radius> （nTran sWidth/2） || n Radius> （nTran sHeight/2））
return (false); // 返回 FALSE
}
pCTData 二new complex<double>[ nTran sWidth * nTran sHeight]; // 分配内存
//存放临时变量
//存放临时变量 //傅立叶变换的宽度（2的整数次幕）
//傅立叶变换的高度（2的整数次幕） // 存贮图像各像素灰度的临时变量
指向时域数据的指针
指向频域数据的指针 （2的整数次幕） （2的整数次幕）
pCFData 二new complex<double>[ nTran sWidth * nTran sHeight]; // 分配内存//图像数据的宽和高不一定是2的整数次幕，所以pCTData有一部分数据需要补
for(y=0; yvnTran sHeight; y++)
{
for(x=0; xvnTran sWidth; x++)
{
pCTData[y* nTran sWidth + x]=complex<double>(0,0); // 补零
}
}
〃把图像数据传给pCTData
for(y=0; yvnH eight; y++)
{
for(x=0; x<n Width; x++)
{
//指向DIB第y行，第x个象素的指针
lpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x;
un chValue = *lpSrc;
pCTData[y* nTran sWidth + x]=complex<double>(u nchValue,0);
}
}
FFT_2D(pCTData, nWidth, nH eight, pCFData) ; // 傅立叶正变换
for(x=0;x< nTran sWidth;x++) // 开始实施理想的低通滤波
{
for(y 二n Radius;y vnTran sHeight-nRadius;y++)
{
II把纵向所有大于nRadius2的高频分量设置为0
pCFData[x* nTran sHeight + y]=complex<double>(0,0);
}
}
for(x 二n Radius;x vnTran sWidth-nRadius;x++)
{
for(y=0 ;y<nTran sHeight;y++)
{
//把横向所有大于nRadiusl的高频分量设置为0
pCFData[x* nTran sHeight + y]=complex<double>(0,0);
}
}
IFFT_2D(pCFData, pCTData,nHeight, nWidth); // 经过低通滤波的图象进行反变换
for(y=0; y< nH eight; y++) // 反变换的数据传给lpDIBBits
{
for(x=0; x<n Width; x++)
{
//需要考虑信号的正负问题以及实际所用的图象数据是灰度值还是原始数据
dReal = pCTData[y* nTran sWidth + x].real() ; // 实部
dImag = pCTData[y* nTran sWidth + x].imag() ; // 虚部
un chValue 二dReal;
//指向DIB第y行，第x个象素的指针
lpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x;
*lpSrc 二un chValue ;
}
}
delete pCTData; // 释放内存
delete pCFData; // 释放内存
pCTData = NULL;
pCFData = NULL;
return (true); //返回结果
}
/*************************************************************************
* /函数名称：
* ButterWorthLowPass()
* /输入参数:
* LPSTR IpDIBBits -指向需要滤波的图像像素指针
* int nWidth - 数据宽度
* int nHeight - 数据高度
* int nRadius - ButterWorth 低通滤波的"半功率"点
* /返回值:
* BOOL //成功返回TRUE，否则返回FALSE。

* /说明:
* pDIBBits是指向需要滤波的图像像素指针。

*经过ButterWorth低通滤波的数据仍然存储在IpDIBBits当中。

**************************************************************************/
BOOL ButterWorthLowPass(LPSTR IpDIBBits, LONG nWidth, LONG nH eight, i nt n Radius)
{
un sig ned char* IpSrc; // 指向源图像的指针
int y ; //循环控制变量
int x ; //循环控制变量
double dTmpOne ; //存放临时变量
double dTmpTwo ; // 存放临时变量
double H ; // ButterWorth 滤波系数
int nTransWidth ; // 傅立叶变换的宽度(2的整数次幕)
int nTransHeight; // 傅立叶变换的高度(2的整数次幕)
double dReal ; // 傅立叶变换的实部
double dImag; // 傅立叶变换的虚部
double un chValue; // 存贮图像各像素灰度的临时变量
complex<double> * pCTData ; // 指向时域数据的指针
complex<double> * pCFData ; //dTmpO ne = log( nWidth)/log(2); dTmpTwo = ceil(dTmp One); dTmpTwo = pow(2,dTmpTwo); nTran sWidth = (int) dTmpTwo;
//计算进行傅立叶变换的点数—纵向 dTmpO ne = log( nH eight)/log(2); dTmpTwo = ceil(dTmp One); dTmpTwo = pow(2,dTmpTwo); nTran sHeight = (int) dTmpTwo;
//低通滤波的半径不能超过频域的最大半径
if(nRadius 〉(nTransWidth/2) || nRadius 〉(nTransHeight/2))
{
return (false); // 返回 FALSE
}
pCTData 二new complex<double>[ nTran sWidth * nTran sHeight]; // pCFData 二new complex<double>[ nTran sWidth * nTran
sHeight]; //
//图像数据的宽和高不一定是 2的整数次幕，所以pCTData 有
for(y=0; yvnTran sHeight; y++)
{
for(x=0; xvnTran sWidth; x++)
{
pCTData[y* nTra nsWidth + x]=complex<double>(0,0); // 补零 }
}
for(y=0; y<nHeight; y++) // 把图像数据传给 pCTData
for(x=0; x<n Width; x++)
{//计算进行傅立叶变换的点数—横向 （2的整数次幕）
指向频域数据的指针
（2的整数次幕）
分配内存 分配内存 部分数据需要补
//指向DIB第y行，第x个象素的指针
IpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x; un chValue = *lpSrc;
pCTData[y* nTran sWidth + x]=complex<double>(u nchValue,O); }
}
FFT_2D(pCTData, nWidth, nH eight, pCFData) ; // 傅立叶正变换for(y=0; y<nTransHeight; y++) // 开始实施ButterWorth 低通滤波{
for(x=0; xvnTran sWidth; x++)
{
H = (double)(y*y+x*x);
H = H / (nRadius * n Radius);
H = 1/(1+H); 〃求H 值
pCFData[y* nTran sWidth + x]=complex<double>(pCFData[y* nTra nsWidth + x].real()*H,
pCFData[y* nTran sWidth + x].imag()*H);
}
}
//经过ButterWorth低通滤波的图象进行反变换
IFFT_2D(pCFData, pCTData, nWidth, nHeight);
for(y=O; y<nHeight; y++) // 反变换的数据传给lpDIBBits
{
for(x=0; x<n Width; x++)
{
dReal = pCTData[y* nTra nsWidth + x].real();
dImag = pCTData[y* nTran sWidth + x].imag();
unchValue = max(0,min(255,sqrt(dReal*dReal+dlmag*dlmag)));
//指向DIB第y行，第x个象素的指针
lpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x;
*lpSrc 二un chValue ;
}
delete pCTData; // 释放内存
delete pCFData; // 释放内存
pCTData = NULL;
pCFData = NULL;
return (true); // 返回结果
}
/*************************************************************************
* /函数名称：
* GaussLowPass()
* /输入参数:
* LPSTR IpDIBBits -指向需要滤波的图像像素指针
* int nWidth - 数据宽度
* int nHeight - 数据高度
* int nRadius - Gauss 低通滤波的"半功率"点
* /返回值:
* BOOL //成功返回TRUE，否则返回FALSE。

* /说明:
* pDIBBits是指向需要滤波的图像像素指针。

*经过Gauss低通滤波的数据仍然存储在IpDIBBits当中。

**************************************************************************/
BOOL GaussLowPass(LPSTR IpDIBBits, LONG nWidth, LONG nH eight, i nt n Radius)
{
un sig ned char* IpSrc; // 指向源图像的指针
int y ; //循环控制变量
int x ; //循环控制变量
double dTmpOne ; //存放临时变量
double dTmpTwo ; // 存放临时变量
double H ; // ButterWorth 滤波系数
int nTransWidth ; // 傅立叶变换的宽度（2的整数次幕）
int nTransHeight;
// 傅立叶变换的高度（2的整数次幕） double dReal ;
// 傅立叶变换的实部 double dImag; // 傅立叶变换的虚部
double un chValue;
// 存贮图像各像素灰度的临时变量 complex<double> * pCTData ;
// complex<double> * pCFData ; // //计算进行傅立叶变换的点数—横向
dTmpO ne = log( nWidth)/log(2);
dTmpTwo = ceil(dTmp One);
dTmpTwo = pow(2,dTmpTwo); nTran
sWidth = (int) dTmpTwo;
//计算进行傅立叶变换的点数—纵向 dTmpO ne = log( nH eight)/log(2); dTmpTwo = ceil(dTmp One); dTmpTwo = pow(2,dTmpTwo); nTran sHeight = (int) dTmpTwo;
//低通滤波的半径不能超过频域的最大半径
if(n Radius> (nTran sWidth/2) || n Radius> (nTran sHeight/2))
{
return (false); // 返回 FALSE
}
pCTData 二new complex<double>[ nTran sWidth * nTran sHeight];
// pCFData 二new complex<double>[ nTran sWidth * nTran
sHeight]; //
//图像数据的宽和高不一定是 2的整数次幕，所以pCTData 有
for(y=0; yvnTran sHeight; y++)
{
for(x=0; xvnTran sWidth; x++)
{
pCTData[y* nTra nsWidth + x]=complex<double>(0,0); // 补零指向时域数据的指针 指向频域数据的指针
（2的整数次幕） （2的整数次幕）
分配内存 分配内存 部分数据需要补
}
for(y=0; yvnHeight; y++) // 把图像数据传给pCTData
{
for(x=0; x<n Width; x++)
{
//指向DIB第y行，第x个象素的指针
IpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x; un chValue = *IpSrc;
pCTData[y* nTran sWidth + x]=complex<double>(u nchVaIue,0);
}
}
FFT_2D(pCTData, nWidth, nH eight, pCFData) ; // 傅立叶正变换for(y=0; y<nTransHeight; y++) // 开始实施Gauss 低通滤波
{
for(x=0; xvnTran sWidth; x++)
{
H = (double)(y*y+x*x);
H = H / (2 * n Radius * n Radius);
H = exp (-H); 〃求H 值
pCFData[y* nTran sWidth + x]=complex<double>(pCFData[y* nTra nsWidth + x].real()*H,
pCFData[y* nTran sWidth + x].imag()*H);
}
}
//经过Gauss低通滤波的图象进行反变换
IFFT_2D(pCFData, pCTData, nWidth, nHeight);
for(y=0; y<nHeight; y++) // 反变换的数据传给lpDIBBits
{
for(x=0; x<n Width; x++)
{
dReal = pCTData[y* nTra nsWidth + x].real();
dImag = pCTData[y* nTran sWidth + x].imag();
unchValue = max(O,min(255,sqrt(dReal*dReal+dlmag*dlmag)));
//指向DIB第y行，第x个象素的指针
lpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x;
*lpSrc 二un chValue ;
}
}
delete pCTData; // 释放内存
delete pCFData; // 释放内存
pCTData = NULL;
pCFData = NULL;
return (true); // 返回结果
}
/*************************************************************************
* /函数名称：
* HighPassFilter()
* /输入参数:
* LPSTR lpDIBBits - 指向需要滤波的图像像素指针
* int nWidth - 数据宽度
* int nHeight - 数据高度
* int nRadius - 理想高通滤波的滤波半径
* /返回值:
* BOOL //成功返回TRUE，否则返回FALSE。

* /说明:
* lpDIBBits是指向需要滤波的图像像素指针。

*经过高通滤波的数据仍然存储在lpDIBBits当中。

**************************************************************************/
BOOL HighPassFilter(LPSTR lpDIBBits, LONG nWidth, LONG nHeight,int n Radius)
un sig ned char* lpSrc; // 指向源图像的指针int y ; //循环控制变量
int x ; //循环控制变量
double dTmp One ; double dTmpTwo ; int nTran sWidth ; int nTran sHeight; double un chValue;
//存放临时变量
//存放临时变量
//傅立叶变换的宽度（2的整数次幕）//傅立叶变换的高度（2的整数次幕）// 存贮图像各像素灰度的临时变量
complex<double> * pCTData ; // complex<double> * pCFData ; // //计算进行傅立叶变换的点数—横向dTmpO ne = log( nWidth)/log(2); dTmpTwo = ceil(dTmp One); dTmpTwo = pow(2,dTmpTwo); nTran sWidth = (int) dTmpTwo;
//计算进行傅立叶变换的点数—纵向指向时域数据的指针指向频域数据的指针（2的整数次幕）
（2的整数次
幕）
dTmpO ne = log( nH eight)/log(2);
dTmpTwo = ceil(dTmp One);
dTmpTwo = pow(2,dTmpTwo);
nTran sHeight = (int) dTmpTwo;
double dReal; // 傅立叶变换的实部
double dImag; // 傅立叶变换的虚部
//滤波的半径不能超过频域的最大半径
if(n Radius〉nTra nsWidth-1 || n Radius> nTra nsHeight-1) {
return (false); // 返回false
}
pCTData 二new complex<double>[ nTran sWidth * nTran sHeight]; // pCFData 二new complex<double>[ nTran sWidth * nTran sHeight]; //
//图像数据的宽和高不一定是2的整数次幕，所以pCTData有
分配内存
分配内存
部分数据需要补
for(y=0; yvnTran sHeight; y++)
{
pCTData[y* nTra nsWidth + x]=complex<double>(0,0); // 补零
}
}
for(y=0; y<nHeight; y++) // 把图像数据传给pCTData
{
for(x=0; x<n Width; x++)
{
//指向DIB第y行，第x个象素的指针
lpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x; un chValue = *lpSrc;
pCTData[y* nTran sWidth + x]=complex<double>(u nchValue,0); }
}
FFT_2D(pCTData, nWidth, nH eight, pCFData) ; // 傅立叶正变换for(x=0;x< nTra nsWidth;x++) // 开始实施理想的高通滤波
{
for(y=0 ;y<n Radius;y++)
{
//把纵向所有小于nRadius的低频分量设置为0
pCFData[x* nTran sHeight + y]=complex<double>(0,0);
}
for(y 二nTran sHeight-1- nRadius;y vnTran sHeight;y++)
{
pCFData[x* nTran sHeight + y]=complex<double>(0,0);
}
}
for(x=0;x <n Radius;x++)
{
//把横向所有小于nRadius的低频分量设置为0
pCFData[x* nTran sHeight + y]=complex<double>(0,0);
}
}
for(x 二nTran sWidth-nRadius;x vnTran sWidth;x++)
{
for(y=0 ;y<nTran sHeight;y++)
{
pCFData[x* nTran sHeight + y]=complex<double>(0,0);
}
}
IFFT_2D(pCFData, pCTData, nHeight, n Width); // 经过高通滤波的图象进行反变换
for(y=0; y<nHeight; y++) // 反变换的数据传给lpDIBBits
{
for(x=0; x<n Width; x++)
{
//需要考虑信号的正负问题以及实际所用的图象数据是灰度值还是原始数据
dReal = pCTData[y* nTra nsWidth + x].real();
dImag = pCTData[y* nTran sWidth + x].imag();
un chValue 二dReal;
//指向DIB第y行，第x个象素的指针
lpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x;
*lpSrc 二un chValue ;
}
}
delete pCTData; // 释放内存
delete pCFData; // 释放内存
pCTData = NULL;
pCFData = NULL;
return (true); // 返回结果
}
/*************************************************************************
* /函数名称：
* ButterWorthHighPass()
* /输入参数:
* LPSTR IpDIBBits -指向需要滤波的图像像素指针
* int nWidth - 数据宽度
* int nHeight - 数据高度
* int nRadius - ButterWorth 高通滤波的"半功率"点
* /返回值:
* BOOL //成功返回TRUE，否则返回FALSE。

* /说明:
* pDIBBits是指向需要滤波的图像像素指针。

*经过ButterWorth低通滤波的数据仍然存储在IpDIBBits当中。

**************************************************************************/
BOOL ButterWorthHighPass(LPSTR IpDIBBits, LONG nWidth, LONG nH eight, i nt n Radius)
{
un sig ned char* IpSrc; // 指向源图像的指针
int y ; //循环控制变量
int x ; //循环控制变量
double dTmpOne ; //存放临时变量
double dTmpTwo ; // 存放临时变量
double H ; // ButterWorth 滤波系数
int nTransWidth ; // 傅立叶变换的宽度(2的整数次幕)
int nTransHeight; // 傅立叶变换的高度(2的整数次幕)
double dReal ; // 傅立叶变换的实部
double dImag; // 傅立叶变换的虚部
double un chValue; // 存贮图像各像素灰度的临时变量
complex<double> * pCTData ; // 指向时域数据的指针complex<double> * pCFData ; // 指向频域数据的指针
//计算进行傅立叶变换的点数-横向(2的整数次幕) dTmpO ne = log( nWidth)/log(2);
dTmpTwo = ceil(dTmp On e);
dTmpTwo = pow(2,dTmpTwo);
nTran sWidth = (int) dTmpTwo;
//计算进行傅立叶变换的点数-纵向(2的整数次幕) dTmpO ne = log( nH eight)/log(2);
dTmpTwo = ceil(dTmp One);
dTmpTwo = pow(2,dTmpTwo);
nTran sHeight = (int) dTmpTwo;
pCTData 二new complex<double>[ nTran sWidth * nTran sHeight]; // pCFData 二new complex<double>[ nTran sWidth * nTran sHeight]; // //图像数据的宽和高不一定是2的整数次幕，所以pCTData有一
for(y=0; yvnTran sHeight; y++)
{
for(x=0; xvnTran sWidth; x++)
{
pCTData[y* nTra nsWidth + x]=complex<double>(0,0); // 补零
}
}
for(y=0; y<nHeight; y++) // 把图像数据传给pCTData
{
for(x=0; x<n Width; x++)
{
//指向DIB第y行，第x个象素的指针
lpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x; un chValue = *lpSrc;
pCTData[y* nTran sWidth + x]=complex<double>(u nchValue,0);
}
分配内存
分配内存
部分数据需要补
}
FFT_2D(pCTData, nWidth, nH eight, pCFData) ; // 傅立叶正变换
for(y=0; yvnTransHeight; y++) // 下面开始实施ButterWorth 高通滤波
{
for(x=0; xvnTran sWidth; x++)
{
H = (double)(y*y+x*x);
H = (n Radius * n Radius) / H ;
H = 1/(1+H); // 求H 值
pCFData[y* nTran sWidth + x]=complex<double>(H*(pCFData[y* nTra nsWidth + x].real()),
H*(pCFData[y* nTran sWidth + x].imag()) );
}
}
//经过ButterWorth高通滤波的图象进行反变换
IFFT_2D(pCFData, pCTData, nWidth, nHeight);
for(y=0; y<nHeight; y++) // 反变换的数据传给lpDIBBits
{
for(x=0; x<n Width; x++)
{
dReal = pCTData[y* nTra nsWidth + x].real();
dImag = pCTData[y* nTran sWidth + x].imag();
un chValue = max(0,mi n(255,sqrt(dReal*dReal+dlmag*dlmag)+100 ));
//指向DIB第y行，第x个象素的指针
lpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x;
*lpSrc 二un chValue ;
}
}
delete pCTData; // 释放内存
delete pCFData; // 释放内存
pCTData = NULL;
pCFData = NULL;
return (true); // 返回结果
* /函数名称：
* GaussHighPass（）
* /输入参数:
* LPSTR IpDIBBits -指向需要滤波的图像像素指针
* int nWidth -数据宽度
* int nHeight - 数据高度
* int nRadius - Gauss 高通滤波的"半功率"点
* /返回值:
* BOOL //成功返回TRUE，否则返回FALSE。

* /说明:
* pDIBBits是指向需要滤波的图像像素指针。

*经过Gauss高通滤波的数据仍然存储在IpDIBBits当中。

**************************************************************************/
BOOL GaussHighPass（LPSTR IpDIBBits, LONG nWidth, LONG nH eight, i nt n Radius）
{
un sig ned char* IpSrc; // 指向源图像的指针
int y ; //循环控制变量
int x ; //循环控制变量
double dTmpOne ; //存放临时变量
double dTmpTwo ; // 存放临时变量
double H ; // ButterWorth 滤波系数
int nTransWidth ; // 傅立叶变换的宽度（2的整数次幕）
int nTransHeight; // 傅立叶变换的高度（2的整数次幕）
double dReal ; // 傅立叶变换的实部
double dImag; // 傅立叶变换的虚部
double un chValue; // 存贮图像各像素灰度的临时变量complex<double> *
pCTData ; // 指向时域数据的指针complex<double> * pCFData ; // 指向频域数据的指针
//计算进行傅立叶变换的点数-横向（2的整数次幕）
dTmpO ne = log( nWidth)/log(2);
dTmpTwo = ceil(dTmp One);
dTmpTwo = pow(2,dTmpTwo);
nTran sWidth = (int) dTmpTwo;
//计算进行傅立叶变换的点数-纵向(2的整数次幕) dTmpO ne = log( nH eight)/log(2);
dTmpTwo = ceil(dTmp One);
dTmpTwo = pow(2,dTmpTwo);
nTran sHeight = (int) dTmpTwo;
//高通滤波的半径不能超过频域的最大半径
if(nRadius〉(nTransWidth/2) || nRadius〉(nTransHeight/2))
{
return (false); // 返回FALSE
}
pCTData 二new complex<double>[ nTran sWidth * nTran sHeight]; // pCFData 二new complex<double>[ nTran sWidth * nTran sHeight]; //
//图像数据的宽和高不一定是2的整数次幕，所以pCTData有
for(y=0; yvnTran sHeight; y++)
{
for(x=0; xvnTran sWidth; x++)
{
pCTData[y* nTra nsWidth + x]=complex<double>(0,0); // 补零
}
}
for(y=0; y<nHeight; y++) // 把图像数据传给pCTData
{
for(x=0; x<n Width; x++)
{
分配内存
分配内存
部分数据需要补
//指向DIB第y行，第x个象素的指针
lpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x;
un chValue = *lpSrc;
pCTData[y* nTran sWidth + x]=complex<double>(u nchValue,O);
}
}
FFT_2D(pCTData, nWidth, nH eight, pCFData) ; // 傅立叶正变换
for(y=0; y<nTransHeight; y++) // 开始实施Gauss 高通滤波
{
for(x=0; xvnTran sWidth; x++)
{
H = (double)(y*y+x*x);
H = H / (2 * n Radius * n Radius);
H = 1 - exp (-H); 〃求H 值
pCFData[y* nTran sWidth + x]=complex<double>(pCFData[y* nTra nsWidth + x].real()*H,
pCFData[y* nTran sWidth + x].imag()*H);
}
}
//经过Gauss高通滤波的图象进行反变换
IFFT_2D(pCFData, pCTData, nWidth, nHeight);
for(y=O; y<nHeight; y++) // 反变换的数据传给lpDIBBits
{
for(x=O; x<n Width; x++)
{
dReal = pCTData[y* nTra nsWidth + x].real();
dImag = pCTData[y* nTran sWidth + x].imag();
unchValue = max(O,min(255,sqrt(dReal*dReal+dlmag*dlmag)));
//指向DIB第y行，第x个象素的指针
lpSrc = (unsigned char*)lpDIBBits + nWidth * (nHeight - 1 - y) + x;
*lpSrc 二un chValue ;
}
}
delete pCTData; // 释放内存delete pCFData; // 释放内存pCTData = NULL;
pCFData = NULL;
return (true); // 返回结果
}。