10. 频率域滤波
10.1 概述及原理详解
频率域滤波的步骤如下:
10.2 低通滤波
低通滤波器越靠近中心点位置的值越接近于 1,越远离中心位置的值就越小于 1,与傅里叶变换相乘后,相当于保留了低频信息,消弱或者移除了高频信息。
10.2.1 理想低通滤波器
假设 𝐻、𝑊 分别代表图像快速傅里叶变换的高、宽。傅里叶谱的最大值(中心点)的位置在 (maxR, maxC)。radius 代表截断频率。𝐷(𝑟, 𝑐) 代表 (𝑟, 𝑐) 到中心位置的距离。
i
l
p
F
i
l
t
e
r
(
r
,
c
)
=
{
1
D
(
r
,
c
)
≤
radius
0
D
(
r
,
c
)
>
radius
\mathbf{ilpFilter}(r,c)=\left\{ \begin{matrix} 1 & D(r,c)\le \text{radius} \\ 0 & D(r,c)>\text{radius} \\ \end{matrix} \right.
ilpFilter(r,c)={10D(r,c)≤radiusD(r,c)>radius
10.2.2 巴特沃斯低通滤波器
b l p F i l t e r ( r , c ) = 1 1 + [ D ( r , c ) / radius ] 2 n \mathbf{blpFilter}(r,c)=\frac{1}{1+{{\left[ D(r,c)/\text{radius} \right]}^{2n}}} blpFilter(r,c)=1+[D(r,c)/radius]2n1
其中,n 代表阶数。
10.2.3 高斯低通滤波器
g l p F i l t e r ( r , c ) = e − D 2 ( r , c ) 2 radiu s 2 \mathbf{glpFilter}(r,c)={{e}^{-\frac{{{D}^{2}}(r,c)}{2\text{radiu}{{\text{s}}^{2}}}}} glpFilter(r,c)=e−2radius2D2(r,c)
通过定义函数 createLPFilter
构建低通滤波器,其中参数 size 代表滤波器的尺寸,即快速傅里叶变换的尺寸;center 代表傅里叶谱的中心位置(即:最大值的位置);radius 代表截断频率;type 代表所定义的枚举类型,enum LPFILTER_TYPE { ILP_FILTER = 0, BLP_FILTER= 1, GLP_FILTER = 2 }
代表三种不同的低通滤波器;n 是只有在构建巴特沃斯滤波器时才用到的参数。
低通滤波器的具体实现如下:
Mat createLPFilter(Size size, Point center, float radius, int type, int n=2)
{
Mat lpFilter = Mat::zeros(size, CV_32FC1);
int rows = size.height;
int cols = size.width;
if (radius <= 0)
return lpFilter;
//构建理想低通滤波器
if (type == ILP_FILTER)
{
for (int r = 0; r < rows; r++)
{
for (int c = 0; c < cols; c++)
{
float norm2 = pow(abs(float(r - center.y)), 2) + pow(abs(float(c - center.x)), 2);
if (sqrt(norm2) < radius)
lpFilter.at<float>(r, c) = 1;
else
lpFilter.at<float>(r, c) = 0;
}
}
}
//构建巴特沃斯低通滤波器
if (type == BLP_FILTER)
{
for (int r = 0; r < rows; r++)
{
for (int c = 0; c<cols; c++)
{
lpFilter.at<float>(r, c) = float(1.0 / (1.0 + pow(sqrt(pow(r - center.y, 2.0) + pow(c - center.x, 2.0)) /radius, 2.0*n)));
}
}
}
//构建高斯低通滤波器
if (type == GLP_FILTER)
{
for (int r = 0; r< rows; r++)
{
for (int c = 0; c < cols; c++)
{
lpFilter.at<float>(r, c) = float(exp(-(pow(c - center.x, 2.0) + pow(r - center.y, 2.0)) / (2 * pow(radius, 2.0))));
}
}
}
return lpFilter;
}
对于低通滤波的实现,将使用在第 10 章中定义的求傅里叶谱的函数 amplitudeSpectrum
和傅里叶谱的灰度级显示函数 graySpectrum
。注意,傅里叶谱在低通滤波中并没有起到任何作用,只是通过傅里叶谱的灰度级显示来观察低通滤波器与傅里叶变换点乘后的灰度级显示是怎样的效果。为了同时观察三种滤波器和截断频率对低通滤波效果的影响,在下面实现中加入了两个进度条来实时调整这两个参数。
低通滤波的实现如下:
#include<opencv2/opencv.hpp>
#include<vector>
#include<string>
#include<iostream>
using namespace cv;
using namespace std;
Mat I;//输入的图像矩阵
Mat F;//图像的快速傅里叶变换
Point maxLoc;//傅里叶谱的最大值的坐标
int radius = 20;//截断频率
const int Max_RADIUS = 100;//设置最大的截断频率
Mat lpFilter;//低通滤波器
int lpType = 0;//低通滤波器的类型
const int MAX_LPTYPE = 2;
Mat F_lpFilter;//低通傅里叶变换
Mat FlpSpectrum;//低通傅里叶变换的傅里叶谱的灰度级
Mat result;//低通滤波后的效果
string lpFilterspectrum = "低通傅里叶谱";//显示窗口的名称
//快速傅里叶变换
void fft2Image(Mat I, Mat &F);
//幅度谱
void amplitudeSpectrum(InputArray _srcFFT, OutputArray _dstSpectrum);
//幅度谱的灰度级显示
Mat graySpectrum(Mat spectrum);
void callback_lpFilter(int, void*);
//低通滤波器类型:理想低通滤波器、巴特沃斯低通滤波器、高斯低通滤波器
enum LPFILTER_TYPE { ILP_FILTER = 0, BLP_FILTER = 1, GLP_FILTER = 2 };
//构建低通滤波器
Mat createLPFilter(Size size, Point center, float radius, int type, int n);
//主函数
int main(int argc, char*argv[])
{
/* -- 第一步:读入图像矩阵 -- */
I = imread("../eye.png", 0);
if (!I.data)
return -1;
imwrite("I1.jpg", I);
//数据类型转换,转换为浮点型
Mat fI;
I.convertTo(fI, CV_32FC1, 1.0, 0.0);
/* -- 第二步:每一个数乘以(-1)^(r+c) -- */
for (int r = 0; r < fI.rows; r++)
{
for (int c = 0; c < fI.cols; c++)
{
if ((r + c) % 2)
{
fI.at<float>(r, c) *= -1;
}
}
}
/* -- 第三、四步:补0和快速傅里叶变换 -- */
fft2Image(fI, F);
//傅里叶谱
Mat amplSpec;
amplitudeSpectrum(F, amplSpec);
//傅里叶谱的灰度级显示
Mat spectrum = graySpectrum(amplSpec);
imshow("原傅里叶谱的灰度级显示", spectrum);
imwrite("spectrum.jpg", spectrum);
//找到傅里叶谱的最大值的坐标
minMaxLoc(spectrum, NULL, NULL, NULL, &maxLoc);
/* -- 低通滤波 -- */
namedWindow(lpFilterspectrum, WINDOW_AUTOSIZE);
createTrackbar("低通类型:", lpFilterspectrum, &lpType, MAX_LPTYPE,
callback_lpFilter);
createTrackbar("半径:", lpFilterspectrum, &radius, Max_RADIUS,
callback_lpFilter);
callback_lpFilter(0, 0);
waitKey(0);
return 0;
}
void fft2Image(Mat I, Mat &F)
{
//得到I的行数和列数
int rows = I.rows;
int cols = I.cols;
//满足快速傅里叶变换的最优行数和列数
int rPadded = getOptimalDFTSize(rows);
int cPadded = getOptimalDFTSize(cols);
//左侧和下侧补0
Mat f;
copyMakeBorder(I, f, 0, rPadded - rows, 0, cPadded - cols, BORDER_CONSTANT, Scalar::all(0));
//单通道转为双通道
Mat planes[] = {f, Mat::zeros(f.size(), CV_32F)};
merge(planes, 2, f);
//快速傅里叶变换(双通道,用于存储实部和虚部)
dft(f, F, DFT_COMPLEX_OUTPUT);
}
void amplitudeSpectrum(InputArray _srcFFT, OutputArray _dstSpectrum)
{
//判断傅里叶变换有两个通道
CV_Assert(_srcFFT.channels() == 2);
//分离通道
vector<Mat> FFT2Channel;
split(_srcFFT, FFT2Channel);
//计算傅里叶变换的幅度谱 sqrt(pow(R,2)+pow(I,2))
magnitude(FFT2Channel[0], FFT2Channel[1], _dstSpectrum);
}
Mat graySpectrum(Mat spectrum)
{
Mat dst;
log(spectrum + 1, dst);
//归一化
normalize(dst, dst, 0, 1, NORM_MINMAX);
//为了进行灰度级显示,做类型转换
dst.convertTo(dst, CV_8UC1, 255, 0);
return dst;
}
//回调函数:调整低通滤波器的类型及截断频率
void callback_lpFilter(int, void*)
{
/* -- 第五步:构建低通滤波器 -- */
lpFilter = createLPFilter(F.size(), maxLoc, radius, lpType, 2);
/*-- 第六步:低通滤波器和图像的快速傅里叶变换点乘 --*/
F_lpFilter.create(F.size(), F.type());
for (int r = 0; r < F_lpFilter.rows; r++)
{
for (int c = 0; c < F_lpFilter.cols; c++)
{
//分别取出当前位置的快速傅里叶变换和理想低通滤波器的值
Vec2f F_rc = F.at<Vec2f>(r, c);
float lpFilter_rc = lpFilter.at<float>(r, c);
//低通滤波器和图像的快速傅里叶变换的对应位置相乘
F_lpFilter.at<Vec2f>(r, c) = F_rc*lpFilter_rc;
}
}
//低通傅里叶变换的傅里叶谱
amplitudeSpectrum(F_lpFilter,FlpSpectrum);
//低通傅里叶谱的灰度级显示
FlpSpectrum = graySpectrum(FlpSpectrum);
imshow(lpFilterspectrum, FlpSpectrum);
imwrite("FlpSpectrum.jpg", FlpSpectrum);
/* -- 第七、八步:对低通傅里叶变换执行傅里叶逆变换,并只取实部 -- */
dft(F_lpFilter, result, DFT_SCALE + DFT_INVERSE + DFT_REAL_OUTPUT);
/* -- 第九步:同乘以(-1)^(x+y) -- */
for (int r = 0; r < result.rows; r++)
{
for (int c = 0; c < result.cols; c++)
{
if ((r + c) % 2)
result.at<float>(r, c) *= -1;
}
}
//注意将结果转换为 CV_8U 类型
result.convertTo(result, CV_8UC1, 1.0, 0);
/* -- 第十步:截取左上部分,其大小与输入图像的大小相同--*/
result = result(Rect(0, 0, I.cols, I.rows)).clone();
imshow("经过低通滤波后的图片", result);
imwrite("lF.jpg", result);
}
Mat createLPFilter(Size size, Point center, float radius, int type, int n=2)
{
Mat lpFilter = Mat::zeros(size, CV_32FC1);
int rows = size.height;
int cols = size.width;
if (radius <= 0)
return lpFilter;
//构建理想低通滤波器
if (type == ILP_FILTER)
{
for (int r = 0; r < rows; r++)
{
for (int c = 0; c < cols; c++)
{
float norm2 = pow(abs(float(r - center.y)), 2) + pow(abs(float(c - center.x)), 2);
if (sqrt(norm2) < radius)
lpFilter.at<float>(r, c) = 1;
else
lpFilter.at<float>(r, c) = 0;
}
}
}
//构建巴特沃斯低通滤波器
if (type == BLP_FILTER)
{
for (int r = 0; r < rows; r++)
{
for (int c = 0; c<cols; c++)
{
lpFilter.at<float>(r, c) = float(1.0 / (1.0 + pow(sqrt(pow(r - center.y, 2.0) + pow(c - center.x, 2.0)) /radius, 2.0*n)));
}
}
}
//构建高斯低通滤波器
if (type == GLP_FILTER)
{
for (int r = 0; r< rows; r++)
{
for (int c = 0; c < cols; c++)
{
lpFilter.at<float>(r, c) = float(exp(-(pow(c - center.x, 2.0) + pow(r - center.y, 2.0)) / (2 * pow(radius, 2.0))));
}
}
}
return lpFilter;
}
注意:高斯低通滤波的效果比其他两者更好,显示更平滑。
10.3 高通滤波
频谱中的高频部分代表的是灰度值迅速变化的区域,通常是图像的边缘。
10.3.1 理想高通滤波器
假设 𝐻、𝑊 分别代表图像快速傅里叶变换的高、宽。傅里叶谱的最大值(中心点)的位置在 (maxR, maxC)。radius 代表截断频率。𝐷(𝑟, 𝑐) 代表 (𝑟, 𝑐) 到中心位置的距离。
i
h
p
F
i
l
t
e
r
(
r
,
c
)
=
{
0
D
(
r
,
c
)
≤
radius
1
D
(
r
,
c
)
>
radius
\mathbf{ihpFilter}(r,c)=\left\{ \begin{matrix} 0 & D(r,c)\le \text{radius} \\ 1 & D(r,c)>\text{radius} \\ \end{matrix} \right.
ihpFilter(r,c)={01D(r,c)≤radiusD(r,c)>radius
10.3.2 巴特沃斯高通滤波器
b h p F i l t e r ( r , c ) = 1 − 1 1 + [ D ( r , c ) / radius ] 2 n \mathbf{bhpFilter}(r,c)=1-\frac{1}{1+{{\left[ D(r,c)/\text{radius} \right]}^{2n}}} bhpFilter(r,c)=1−1+[D(r,c)/radius]2n1
其中,n 代表阶数。
10.3.3 高斯高通滤波器
g h p F i l t e r ( r , c ) = 1 − e − D 2 ( r , c ) 2 radiu s 2 \mathbf{ghpFilter}(r,c)=1-{{e}^{-\frac{{{D}^{2}}(r,c)}{2\text{radiu}{{\text{s}}^{2}}}}} ghpFilter(r,c)=1−e−2radius2D2(r,c)
对于图像的高通滤波,只需要将 10.2.3 中低通滤波的 C++ 实现中的第五步替换成高通滤波器就可以了。
10.4 带通滤波
假设 BW \text{BW} BW 代表带宽, D 0 D_0 D0 代表带宽的径向中心,其他符号与低通、高通滤波相同。
10.4.1 理想带通滤波器
假设 𝐻、𝑊 分别代表图像快速傅里叶变换的高、宽。傅里叶谱的最大值(中心点)的位置在 (maxR, maxC)。radius 代表截断频率。𝐷(𝑟, 𝑐) 代表 (𝑟, 𝑐) 到中心位置的距离。
KaTeX parse error: Unknown column alignment: * at position 48: … \begin{array}{*̲{35}{l}} 1 &…
10.4.2 巴特沃斯带通滤波器
b b p F i l t e r ( r , c ) = 1 − 1 1 + [ D ( r , c ) × BW D ( r , c ) 2 − D 0 2 ] 2 n \mathbf{bbpFilter}(r,c)=1-\frac{1}{1+{{\left[ \frac{D(r,c)\times \text{BW}}{D{{(r,c)}^{2}}-{{D}_{0}}^{2}} \right]}^{2n}}} bbpFilter(r,c)=1−1+[D(r,c)2−D02D(r,c)×BW]2n1
其中,n 代表阶数。
10.4.3 高斯带通滤波器
g b p F i l t e r ( r , c ) = e − ( D ( r , c ) 2 − D 0 2 D ( r , c ) × BW ) 2 \mathbf{gbpFilter}(r,c)={{e}^{-{{(\frac{D{{(r,c)}^{2}}-{{D}_{0}}^{2}}{D(r,c)\times \text{BW}})}^{2}}}} gbpFilter(r,c)=e−(D(r,c)×BWD(r,c)2−D02)2
10.5 带阻滤波
假设 BW \text{BW} BW 代表带宽, D 0 D_0 D0 代表带宽的径向中心,其他符号与低通、高通滤波相同。
10.5.1 理想带阻滤波器
假设 𝐻、𝑊 分别代表图像快速傅里叶变换的高、宽。傅里叶谱的最大值(中心点)的位置在 (maxR, maxC)。radius 代表截断频率。𝐷(𝑟, 𝑐) 代表 (𝑟, 𝑐) 到中心位置的距离。
KaTeX parse error: Unknown column alignment: * at position 48: … \begin{array}{*̲{35}{l}} 0 &…
10.5.2 巴特沃斯带阻滤波器
b b r F i l t e r ( r , c ) = 1 1 + [ D ( r , c ) × BW D ( r , c ) 2 − D 0 2 ] 2 n \mathbf{bbrFilter}(r,c)=\frac{1}{1+{{\left[ \frac{D(r,c)\times \text{BW}}{D{{(r,c)}^{2}}-{{D}_{0}}^{2}} \right]}^{2n}}} bbrFilter(r,c)=1+[D(r,c)2−D02D(r,c)×BW]2n1
其中,n 代表阶数。
10.5.3 高斯带阻滤波器
g b r F i l t e r ( r , c ) = 1 − e − ( D ( r , c ) 2 − D 0 2 D ( r , c ) × BW ) 2 \mathbf{gbrFilter}(r,c)=1-{{e}^{-{{(\frac{D{{(r,c)}^{2}}-{{D}_{0}}^{2}}{D(r,c)\times \text{BW}})}^{2}}}} gbrFilter(r,c)=1−e−(D(r,c)×BWD(r,c)2−D02)2
10.6 自定义滤波器
低通、高通、带通、带阻滤波均是按照某一特定规则构建滤波器的。自定义滤波器通常用于消除结构化噪声或者目标。(具体略)
10.7 同态滤波
同态滤波与频率域滤波的不同之处是,它在最开始对输入的图像矩阵进行对数运算,在最后一步进行对数运算的逆运算,即指数运算,其中间步骤就是频率域滤波的步骤。
使用同态滤波处理后可以看到原图中更多的信息。