频域滤波
数字图像本质上是数字信号，图像的滤波就是对信号的不同频率进行筛选。体现在图像上，模糊操作实际是抑制高频，锐化操作实际是抑制低频。
因此，从频域上对图像信号做滤波更贴近信号处理的本质，同时也能实现空间域滤波很难达到的效果（如带通带阻）。
那么，如何将一幅数字图像在空间域和频率域中做转换呢？离散傅里叶变换（DFT）
1. 离散傅里叶变换
傅里叶提出：任何周期函数都可以用一系列正余弦函数的累加来逼近。
因此对于一个空间域的信号（比如图像），可以用频率、幅值、相位来描述
上面两个公式即为二维离散傅里叶变换（具体的理论推导是数字信号处理中的内容，网上内容很多，冈萨雷斯的数字图像处理也有详述），u,v表示频率采样位置（在这里将频率范围从0-2PI换到图像的高/宽），x,y表示空间采样位置。
由复变函数的知识可以知道，e^jsita可以表示成实部+虚部，这样上述傅里叶正变换可表示为
图1-1. DFT（左：原图，中：幅值，右：相位）
上图使用不同的素材做DFT，相应的傅里叶谱和相位图都做了中心化转换（频率为0的点由左上角移至中心，具体做法是将得到的傅里叶谱图切分成2*2的子窗口，左上角的窗口和右下角的窗口互换，右上角的窗口和左下角的窗口互换；当然也可以通过将原图 * (-1)^(x+y)实现中心化转换）。
可以看到，在傅里叶谱中中心的点（代表频率为0）灰度值最高，其表示的是图像的平均灰度值；从整体形状上看，傅里叶谱指示着原图梯度的方向，并且这种指示不随图像的平移发生变化，然而此时相位图会发生很大变化。
从公式可以看到，DFT的每个点都需要遍历所有像素，时间复杂度为O(MN*MN)，很显然无法实用，好在还有快速傅里叶变换（FFT），其时间复杂度为O(MN*log(MN))。FFT的具体算法可以参考FFT算法的完整DSP实现。
2. 频域滤波器
在空间域上实现的平滑、锐化滤波器对应频域上的低通、高通滤波器，两者的效果相似，但应用了FFT的频域操作较空间域操作更便捷（空间域上的卷积对应频域上的点乘）
此外频域操作还能实现空间域不能完成的带通、带阻滤波（这是频域滤波的意义所在）
基本步骤
1）对原图的高宽进行扩充，一般是原始高宽均扩大2倍，多出来的部分补0。这么做的原因是为了将数字图像模拟成周期信号时避免信号混淆。
2）对扩充后的图像 f(x,y)，以
3）设计一个频域滤波函数 H(u,v)，其尺寸与扩充后的图像尺寸一致，中心为低频。将FFT后的图像 F(u,v) 与 H(u,v)点乘得滤波后的结果G(u,v)。
4）对G(u,v)做反傅里叶变换以转换到空间域，取其实部的值做归一化等操作为图像的值（不使用虚部值是为了避免由于计算误差导致的寄生复变量；如果之前是用
5）对图像做裁剪，取自左上角开始的原始图像高宽的矩形区域为最终图像。
常用滤波器
1）低通滤波器
左边为布特沃斯低通滤波器，右边是高斯低通滤波器。D(u,v)为点(u,v)到中心的欧式距离，D0为截止频率
2）高通滤波器
左边为布特沃斯高通滤波器，右边是高斯高通滤波器。
3）拉普拉斯增强
反傅里叶变换后即为原图像的拉普拉斯变换，再加上原图即实现图像增强
公式
类似于高频强调滤波器，同态滤波器也能衰减低频增强高频
4）选择性滤波器
左边为布特沃斯带阻滤波器，右边是高斯带阻滤波器，W是带宽。用1减去上述公式即得相应的带通滤波器。
5）陷波滤波器
类似的，
图2-1. 1个陷波对的2阶布特沃斯陷波带阻滤波
上图是使用1个陷波对的2阶布特沃斯陷波带阻滤波的部分效果。左上是原图，由于莫尔效应，上面包含了具有周期性的网状噪声（在空间域上很难消除），左下是原图的频谱，可以看到网状噪声反映在频谱上是关于中心对称的多组亮点（每组包含4个亮点，对应1个陷波对，这里仅对最亮的一组做处理）。显然，只需我们的频域滤波器在这些亮点的部分取低值就可压制噪声。上图中间为当陷波中心取15，截止频率取20的效果，可以看到有4个代表陷波的黑点。上图右边继续调整陷波中心（29），直至其覆盖4个亮点，增大截止频率（149）以扩大黑点半径，可以看到，经过滤波后图像的网状噪声得到有效抑制。
3. 代码
FFT的实现网上有很多，这里参考了Paul Bourke的实现
顺便贴个OpenCV的例程（opencv/samples/cpp/dft.cpp）
<span style="white-space:pre"> </span>Mat padded;
<span style="white-space:pre"> </span>//expand input image to optimal size
int m = getOptimalDFTSize(in.rows);
int n = getOptimalDFTSize(in.cols);
// on the border add zero values
copyMakeBorder(in, padded, 0, m - in.rows, 0, n - in.cols, BORDER_CONSTANT, Scalar::all(0));
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
Mat complexI;
merge(planes, 2, complexI); // Add to the expanded another plane with zeros
dft(complexI, complexI); // this way the result may fit in the source matrix
// compute the magnitude and switch to logarithmic scale
// => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
split(complexI, planes); // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
magnitude(planes[0], planes[1], planes[0]);// planes[0] = sqrt(planes[0]^2 + planes[1]^2)
Mat magI = planes[0];
magI += Scalar::all(1); // switch to logarithmic scale
log(magI, magI);
// crop the spectrum, if it has an odd number of rows or columns
magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));
// rearrange the quadrants of Fourier image so that the origin is at the image center
int cx = magI.cols / 2;
int cy = magI.rows / 2;
Mat q0(magI, Rect(0, 0, cx, cy)); // Top-Left - Create a ROI per quadrant
Mat q1(magI, Rect(cx, 0, cx, cy)); // Top-Right
Mat q2(magI, Rect(0, cy, cx, cy)); // Bottom-Left
Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right
Mat tmp; // swap quadrants (Top-Left with Bottom-Right)
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp); // swap quadrant (Top-Right with Bottom-Left)
q2.copyTo(q1);
tmp.copyTo(q2);
normalize(magI, magI, 0, 1, CV_MINMAX); // Transform the matrix with float values into a
// viewable image form (float between values 0 and 1).
magI.copyTo(out);4. 参考
[1] 冈萨雷斯 《数字图像处理》
[2] Paul Bourke: Fast Fourier Transform
[3] http://blog.csdn.net/abcjennifer/article/details/7622228
1831
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
