空间滤波&频域滤波(2)

频域滤波

数字图像本质上是数字信号，图像的滤波就是对信号的不同频率进行筛选。体现在图像上，模糊操作实际是抑制高频，锐化操作实际是抑制低频。

因此，从频域上对图像信号做滤波更贴近信号处理的本质，同时也能实现空间域滤波很难达到的效果（如带通带阻）。

那么，如何将一幅数字图像在空间域和频率域中做转换呢？离散傅里叶变换（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)，以相乘，这样转换得到的傅里叶谱就是以中心点为原点的了；或者也可以省略这步，将得到的傅里叶谱（此时四个角的位置频率最低）左上角与右下角互换，右上角与左下角互换（OpenCV的例程就是这么处理的）。

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）拉普拉斯增强

反傅里叶变换后即为原图像的拉普拉斯变换，再加上原图即实现图像增强

公式是个更一般的高频强调滤波方案，k2控制高频贡献。

类似于高频强调滤波器，同态滤波器也能衰减低频增强高频

4）选择性滤波器

左边为布特沃斯带阻滤波器，右边是高斯带阻滤波器，W是带宽。用1减去上述公式即得相应的带通滤波器。

5）陷波滤波器

为一个包含了3个陷波对（考虑频域的对称关系，一个中心位于(u,v)的陷波在(-u,-v)上有一个对应的陷波）的n 阶布特沃斯陷波带阻滤波器。其中为每组陷波对的截止频率，为陷波中心位置。

类似的，为对应的带通滤波器。

图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 sizeint m = getOptimalDFTSize(in.rows);int n = getOptimalDFTSize(in.cols); // on the border add zero valuescopyMakeBorder(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 zerosdft(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 scalelog(magI, magI);// crop the spectrum, if it has an odd number of rows or columnsmagI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));// rearrange the quadrants of Fourier image  so that the origin is at the image centerint cx = magI.cols / 2;int cy = magI.rows / 2;Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrantMat q1(magI, Rect(cx, 0, cx, cy));  // Top-RightMat q2(magI, Rect(0, cy, cx, cy));  // Bottom-LeftMat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-RightMat 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