【图像处理】-013 频域滤波处理
在上一篇中,我们讲到了进行同态滤波对图像不均匀光照进行处理,其中用到了频域中对图像的高通滤波。这一篇文章中,我将对图像频域滤波进行进一步讨论,尝试使用各种低通和高通滤波器对图像进行滤波处理。
文章目录
1 理论依据
1.1 频率域滤波基础
在二维离散傅立叶变换DFT中,
(1)
F
(
u
,
v
)
=
∑
x
=
0
M
−
1
∑
y
=
0
N
−
1
f
(
x
,
y
)
e
−
j
2
π
(
u
x
/
M
+
v
y
/
N
)
F(u,v)=\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)e^{-j2\pi (ux/M+ vy/N)} \tag{1}
F(u,v)=x=0∑M−1y=0∑N−1f(x,y)e−j2π(ux/M+vy/N)(1)
可以看出,
F
(
u
,
v
)
F(u,v)
F(u,v)的每一项都包含有用指数形式修改过的
f
(
x
,
y
)
f(x,y)
f(x,y)的所有值,因此,将一幅图像的特定分量与其变换直接联系起来通常是不可能的。然而,关于傅立叶变换的频率分量和一幅图像的空间特性间的关系可以做出某些一般的叙述。
频率域的滤波技术是以如下处理为基础的:修改傅立叶变换以达到特殊目的,然后计算IDFT返回到图像域。遵循 F ( u , v ) = ∣ F ( u , v ) ∣ e j ϕ ( u , v ) F(u,v)=|F(u,v)|e^{j\phi(u,v)} F(u,v)=∣F(u,v)∣ejϕ(u,v),我们选择对傅立叶变换的幅度(谱)和相角进行处理。
频率域滤波由修改一幅图像的傅立叶变换然后计算其反变换得到处理后的结果组成。
若给定一幅大小为
M
×
N
M\times N
M×N的数字图像
f
(
x
,
y
)
f(x,y)
f(x,y),则基本滤波公式如下:
(2)
g
(
x
,
y
)
=
ℑ
−
1
[
H
(
u
,
v
)
F
(
u
,
v
)
]
g(x,y)=\Im^{-1}[H(u,v)F(u,v)] \tag{2}
g(x,y)=ℑ−1[H(u,v)F(u,v)](2)
其中,
ℑ
−
1
\Im^{-1}
ℑ−1是IDFT,
F
(
u
,
v
)
F(u,v)
F(u,v)是输入图像
f
(
x
,
y
)
f(x,y)
f(x,y)的DFT,
H
(
u
,
v
)
H(u,v)
H(u,v)是滤波函数(也称为滤波器,或者滤波传递函数),
g
(
x
,
y
)
g(x,y)
g(x,y)是滤波后的输出图像。函数
F
,
H
,
g
F,H,g
F,H,g是大小和输入图像相同的
M
×
N
M\times N
M×N阵列。乘积
[
H
(
u
,
v
)
F
(
u
,
v
)
]
[H(u,v)F(u,v)]
[H(u,v)F(u,v)]是由阵列相乘形成的。滤波函数修改输入图像的变换来得到处理后的输出
g
(
x
,
y
)
g(x,y)
g(x,y)。使用关于中心对称的函数可以显著的简化
H
(
u
,
v
)
H(u,v)
H(u,v)的技术条件,它要求
F
(
u
,
v
)
F(u,v)
F(u,v)也被中心化。这是通过在变换前使用
(
−
1
)
x
+
y
(-1)^{x+y}
(−1)x+y乘以输入图像来完成的。
如果 H H H是实对称函数而 f f f是实函数(通常如此),IDFT理论上应生成实数量。实际上,该反变换通常包含由舍入误差和其他计算错误引起的寄生复数项。通常取IDFT得实部来形成函数 g g g。
傅立叶变换中的低频与图像中缓慢变化的灰度分量有关,例如室内的墙面和室外少云的天空等。另一方面,高频由灰尘的尖锐过渡造成,如边缘和噪声。因此,衰减高频而通过低频的滤波器 H ( u , v ) H(u,v) H(u,v)(近似低通滤波器)将模糊一幅图像,相反特性的滤波器称为高通滤波器将增强尖锐的细节,但将导致图像对比度降低。
在式 ( 2 ) (2) (2)中,涉及到频率域中两个函数的乘积,由卷积定理可知,频率域中的乘积等于空间域内的卷积。
在DFT复数阵列中,傅立叶变换的复数形式如下:
(3)
F
(
u
,
v
)
=
R
(
u
,
v
)
+
j
I
(
u
,
v
)
F(u,v)=R(u,v)+jI(u,v) \tag{3}
F(u,v)=R(u,v)+jI(u,v)(3)
然后
(4)
g
(
x
,
y
)
=
ℑ
−
1
[
H
(
u
,
v
)
R
(
u
,
v
)
+
j
H
(
u
,
v
)
I
(
i
,
v
)
]
g(x,y)=\Im^{-1}[H(u,v)R(u,v)+jH(u,v)I(i,v)] \tag{4}
g(x,y)=ℑ−1[H(u,v)R(u,v)+jH(u,v)I(i,v)](4)
由于相角计算时,
(5)
ϕ
(
u
,
v
)
=
a
r
c
t
a
n
[
I
m
g
(
u
,
v
)
R
e
a
l
(
u
,
v
)
]
\phi (u,v) = arctan[\frac{Img(u,v)}{Real(u,v)}] \tag{5}
ϕ(u,v)=arctan[Real(u,v)Img(u,v)](5)
从上式可以看出,经过滤波之后,相位角并不会改变,因为计算相位角时,
H
(
u
,
v
)
H(u,v)
H(u,v)被相除抵消了。等同的影响实部和虚部而不影响相位的滤波器成为零相移滤波器。
1.2 阵列相乘
阵列相乘表示两个矩阵的对应元素分别相乘,需要与矩阵乘法相区别。在MATLAB中,阵列相乘通过矩阵的
.
∗
.*
.∗实现。在OpenCV中,阵列相乘通过cv::Mat::mul()
函数实现,该函数执行了矩阵的阵列相乘,而cv::Mat::dot()
函数执行的是向量的点积,对于矩阵调用该接口,会将矩阵看成1行N列的向量,执行点积操作,返回一个double值。
1.3 频率域滤波步骤小结
- (1) 给定一幅大小为 M × N M\times N M×N的输入图像 f ( x , y ) f(x,y) f(x,y),按照 P ≥ 2 M − 1 P \geq 2M-1 P≥2M−1和 Q ≥ 2 N − 1 Q \geq 2N-1 Q≥2N−1得到填充参数 P P P和 Q Q Q,典型地,可以选择 P = 2 M P=2M P=2M和 Q = 2 N Q=2N Q=2N;
- (2) 对 f ( x , y ) f(x,y) f(x,y)添加必要数量的0,形成大小为 P × Q P\times Q P×Q的填充图像 f p ( x , y ) f_p(x,y) fp(x,y);
- (3) 用 ( − 1 ) x + y (-1)^{x+y} (−1)x+y乘以 f p ( x , y ) f_p(x,y) fp(x,y)移到其变换的中心;
- (4) 计算来自步骤(3)的图像的DFT,得到 F ( u , v ) F(u,v) F(u,v);
- (5) 生成一个实的、对称的滤波函数 H ( u , v ) H(u,v) H(u,v),其大小为 P × Q P\times Q P×Q,中心在 ( P / 2 , Q / 2 ) (P/2,Q/2) (P/2,Q/2)处。用阵列相乘形成乘积 G ( u , v ) = H ( u , v ) F ( u , v ) G(u,v)=H(u,v)F(u,v) G(u,v)=H(u,v)F(u,v),即 G ( i , k ) = H ( i , k ) F ( i , k ) G(i,k)=H(i,k)F(i,k) G(i,k)=H(i,k)F(i,k);
- (6) 得到处理后的图像:
g p ( x , y ) = { r e a l [ ℑ − 1 [ G ( u , v ) ] ] } ( − 1 ) ( x + y ) g_p(x,y)=\{real[\Im^{-1}[G(u,v)]]\}(-1)^(x+y) gp(x,y)={real[ℑ−1[G(u,v)]]}(−1)(x+y)
为了忽略由于计算不准确导致的寄生分量,选择了实部,并且处理的是填充后的阵列; - (7) 通过从 g ( x , y ) g(x,y) g(x,y)的左上象限提取 M × N M\times N M×N区域,得到最终处理结果 g ( x , y ) g(x,y) g(x,y)。
1.4 空间与频率滤波间的对应
空间域和频率域滤波间的纽带时卷积定理。频率域滤波的定义中,滤波的结果等于滤波函数 H ( u , v ) H(u,v) H(u,v)和输入图像的傅立叶变换 F ( u , v ) F(u,v) F(u,v)的乘积,而频域中的乘积等于空域的卷积,因此,空间域滤波是滤波函数与输入图像的卷积。
p 对于给定的频率域滤波器 H ( u , v ) H(u,v) H(u,v),如果我们要给出其空间域的等价形式,可以令图像 f ( x , y ) = δ ( x , y ) f(x,y)=\delta (x,y) f(x,y)=δ(x,y),那么 F ( u , v ) = 1 F(u,v)=1 F(u,v)=1,此时,由 g ( x , y ) = ℑ − 1 [ H ( u , v ) F ( u , v ) ] g(x,y)=\Im^{-1}[H(u,v)F(u,v)] g(x,y)=ℑ−1[H(u,v)F(u,v)]可知,滤波器的输出为 ℑ − 1 [ H ( u , v ) ] \Im^{-1}[H(u,v)] ℑ−1[H(u,v)],这就是频域滤波器的傅立叶反变换,此时,它对应的是空间域的滤波器。
因此,两个滤波器形成了傅立叶变换对:
(6)
h
(
x
,
y
)
⇔
H
(
u
,
v
)
h(x,y)\Leftrightarrow H(u,v) \tag{6}
h(x,y)⇔H(u,v)(6)
其中
h
(
x
,
y
)
h(x,y)
h(x,y)是一个空间滤波器,因为该滤波器可以由频率域滤波器对一个冲激的响应得到,所以
h
(
x
,
y
)
h(x,y)
h(x,y)也称为
H
(
u
,
v
)
H(u,v)
H(u,v)的脉冲响应,因为式6的离散实现中的所有数值都是有限的,这样的滤波器称为有限冲激响应(FIR)滤波器。
2 频域低通滤波
2.1 理想低通滤波器
2.1.1 数学表示
在以原点为圆心,以
D
0
D_0
D0为半径的圆内,无衰减的通过所有频率,而在该圆外“切断”所有频率的二维低通滤波器,称为理想低通滤波器(ILPF),它的函数表示如下:
(7)
H
(
u
,
v
)
=
{
1
,
D
(
u
,
v
)
≤
D
0
0
,
D
(
u
,
v
)
>
D
0
H(u,v)=\begin{cases} 1,&\qquad D(u,v) \leq D_0 & \\ 0,&\qquad D(u,v) > D_0 \end{cases} \tag{7}
H(u,v)={1,0,D(u,v)≤D0D(u,v)>D0(7)
其中,
D
0
D_0
D0是一个正常数,
D
(
u
,
v
)
D(u,v)
D(u,v)是频率域中点
(
u
,
v
)
(u,v)
(u,v)与频率矩形中心的距离,即:
(8)
D
(
u
,
v
)
=
[
(
u
−
P
/
2
)
2
+
(
v
−
Q
/
2
)
2
]
1
/
2
D(u,v)=[(u-P/2)^2 + (v-Q/2)^2]^{1/2} \tag{8}
D(u,v)=[(u−P/2)2+(v−Q/2)2]1/2(8)
其中,
P
P
P和
Q
Q
Q是按照
P
≥
2
M
−
1
P \geq 2M-1
P≥2M−1和
Q
≥
2
N
−
1
Q \geq 2N-1
Q≥2N−1填充后的尺寸。
对于一个理想低通滤波器(ILPF),在 H ( u , v ) = 1 H(u,v)=1 H(u,v)=1和 H ( u , v ) = 0 H(u,v)=0 H(u,v)=0之间的过渡点称为截止频率。
2.1.2 图示
这里,我选择使用MATLAB来进行理想低通滤波器本身的创建和图示。
a=100;
b=100;
U=0:a;
V=0:b;
M=length(U);N=length(V);
D0=10;
m=fix(M/2); n=fix(N/2);
H=zeros(M,N);
n=8;
for u=1:M
for v=1:N
a=sqrt((U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v) - 50));%D(u,v)的值
if(a<=D0)%理想低通滤波器,在圆内的频率完全通过,圆外的频率完全截止
H(u,v)=1;
else
H(u,v)=0;
end
end
end
%在绘制高斯曲面的时候,加上下述代码,显示得美观
figure;
surf(U,V,H)
![](https://i-blog.csdnimg.cn/blog_migrate/e953f8c516d126fdc271bec8ecfb08b8.png)
Figure 1. 理想低通滤波器
2.2 巴特沃斯低通滤波器
2.2.1 数学表示
截止频率位于距离原点
D
0
D_0
D0处的
n
n
n阶巴特沃斯低通滤波器(BLPF)的传递函数定义为:
H
(
u
,
v
)
=
1
1
+
[
D
(
u
,
v
)
/
D
0
]
2
n
H(u,v)=\frac{1}{1+[D(u,v)/D_0]^{2n}}
H(u,v)=1+[D(u,v)/D0]2n1
2.2.2 图示
a=100;
b=100;
U=0:a;
V=0:b;
M=length(U);N=length(V);
D0=10;
n = 1;
H=zeros(M,N);
for u=1:M
for v=1:N
a=sqrt((U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v) - 50));
H(u,v)=1/(1+(a/D0)^(2*n));
end
end
%在绘制高斯曲面的时候,加上下述代码,显示得美观
figure;
surf(U,V,H),title('1阶巴特沃斯低通滤波器');
2.3 高斯低通滤波器
2.3.1 数学表示
高斯低通滤波器(GLPF)的数学表达式如下:
H
(
u
,
v
)
=
e
−
D
2
(
u
,
v
)
/
2
σ
2
H(u,v)=e^{-D^{2}(u,v)/2\sigma ^2}
H(u,v)=e−D2(u,v)/2σ2
通常讨论时,可以去截止频率
D
0
D_0
D0,表示形式如下:
H
(
u
,
v
)
=
e
−
D
2
(
u
,
v
)
/
2
D
0
2
H(u,v)=e^{-D^{2}(u,v)/2D_0^2}
H(u,v)=e−D2(u,v)/2D02
2.3.2 图示
a=100;
b=100;
U=0:a;
V=0:b;
M=length(U);N=length(V);
D0=10;
H=zeros(M,N);
e = exp(1);
for u=1:M
for v=1:N
a=sqrt((U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v) - 50));
H(u,v)=e^(-a*a/(2*D0*D0));
end
end
%在绘制高斯曲面的时候,加上下述代码,显示得美观
figure;
surf(U,V,H),title('高斯低通滤波器');
3 频域高通滤波
3.1 理想高通滤波器
3.1.1 数学表示
在以原点为圆心,以
D
0
D_0
D0为半径的圆外,无衰减的通过所有频率,而在该圆内“切断”所有频率的二维低通滤波器,称为理想高通滤波器(IHPF),它的函数表示如下:
(7)
H
(
u
,
v
)
=
{
0
,
D
(
u
,
v
)
≤
D
0
1
,
D
(
u
,
v
)
>
D
0
H(u,v)=\begin{cases} 0,&\qquad D(u,v) \leq D_0 & \\ 1,&\qquad D(u,v) > D_0 \end{cases} \tag{7}
H(u,v)={0,1,D(u,v)≤D0D(u,v)>D0(7)
其中,
D
0
D_0
D0是一个正常数,
D
(
u
,
v
)
D(u,v)
D(u,v)是频率域中点
(
u
,
v
)
(u,v)
(u,v)与频率矩形中心的距离,即:
(8)
D
(
u
,
v
)
=
[
(
u
−
P
/
2
)
2
+
(
v
−
Q
/
2
)
2
]
1
/
2
D(u,v)=[(u-P/2)^2 + (v-Q/2)^2]^{1/2} \tag{8}
D(u,v)=[(u−P/2)2+(v−Q/2)2]1/2(8)
其中,
P
P
P和
Q
Q
Q是按照
P
≥
2
M
−
1
P \geq 2M-1
P≥2M−1和
Q
≥
2
N
−
1
Q \geq 2N-1
Q≥2N−1填充后的尺寸。
3.1.2 图示
a=100;
b=100;
U=0:a;
V=0:b;
M=length(U);N=length(V);
D0=10;
H=zeros(M,N);
for u=1:M
for v=1:N
a=sqrt((U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v) - 50));%D(u,v)的值
if(a<=D0)%理想高通滤波器,在圆外的频率完全通过,圆内的频率完全截止
H(u,v)=0;
else
H(u,v)=1;
end
end
end
%在绘制高斯曲面的时候,加上下述代码,显示得美观
figure;
surf(U,V,H)
3.2 巴特沃斯高通滤波器
3.2.1 数学表示
截止频率位于距离原点
D
0
D_0
D0处的
n
n
n阶巴特沃斯低通滤波器(BLPF)的传递函数定义为:
H
(
u
,
v
)
=
1
1
+
[
D
0
/
D
(
u
,
v
)
]
2
n
H(u,v)=\frac{1}{1+[D_{0}/D(u,v)]^{2n}}
H(u,v)=1+[D0/D(u,v)]2n1
3.2.2 图示
a=100;
b=100;
U=0:a;
V=0:b;
M=length(U);N=length(V);
D0=10;
n = 1;
H=zeros(M,N);
for u=1:M
for v=1:N
a=sqrt((U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v) - 50));
H(u,v)=1/(1+(D0/a)^(2*n));
end
end
%在绘制高斯曲面的时候,加上下述代码,显示得美观
figure;
surf(U,V,H),title('1阶巴特沃斯高通滤波器');
3.3 高斯高通滤波器
3.3.1 数学表示
高斯低通滤波器(GLPF)的数学表达式如下:
H
(
u
,
v
)
=
1
−
e
−
D
2
(
u
,
v
)
/
2
σ
2
H(u,v)=1-e^{-D^{2}(u,v)/2\sigma ^2}
H(u,v)=1−e−D2(u,v)/2σ2
通常讨论时,可以去截止频率
D
0
D_0
D0,表示形式如下:
H
(
u
,
v
)
=
1
−
e
−
D
2
(
u
,
v
)
/
2
D
0
2
H(u,v)=1-e^{-D^{2}(u,v)/2D_0^2}
H(u,v)=1−e−D2(u,v)/2D02
3.3.2 图示
a=100;
b=100;
U=0:a;
V=0:b;
M=length(U);N=length(V);
D0=10;
H=zeros(M,N);
e = exp(1);
for u=1:M
for v=1:N
a=sqrt((U(u) - 50) .* (U(u)-50) + (V(v) - 50) .* (V(v) - 50));
H(u,v)=1-e^(-a*a/(2*D0*D0));
end
end
%在绘制高斯曲面的时候,加上下述代码,显示得美观
figure;
surf(U,V,H),title('高斯高通滤波器');
4 OpenCV实现效果
struct SILPFParam
{
cv::Size sz;
float r;
};
typedef SILPFParam SIHPFParam;
struct SBLPFParam
{
cv::Size sz;
float r;
int n;
};
typedef SBLPFParam SBHPFParam;
void CreateIdealLowpassFilter(SILPFParam& param, cv::Mat& ilpf)
{
cv::Mat single(param.sz.height, param.sz.width, CV_32F);
cv::Point centre = cv::Point(param.sz.height / 2, param.sz.width / 2);
double radius;
for (int i = 0; i < param.sz.height; i++)
{
for (int j = 0; j < param.sz.width; j++)
{
radius = sqrt(pow((float)(i - centre.x), 2) + pow((float)(j - centre.y), 2));
if (radius <= param.r)
single.at<float>(i, j) = 1;
else
single.at<float>(i, j) = 0;
}
}
cv::Mat IDLP_channels[] = { cv::Mat_<float>(single),cv::Mat_<float>(single) };
cv::merge(IDLP_channels, 2, ilpf);
}
void CreateIdealHighpassFilter(SIHPFParam& param, cv::Mat& ihpf)
{
cv::Mat single(param.sz.height, param.sz.width, CV_32F);
cv::Point centre = cv::Point(param.sz.height / 2, param.sz.width / 2);
double radius;
for (int i = 0; i < param.sz.height; i++)
{
for (int j = 0; j < param.sz.width; j++)
{
radius = sqrt(pow((float)(i - centre.x), 2) + pow((float)(j - centre.y), 2));
if (radius <= param.r)
single.at<float>(i, j) = 0;
else
single.at<float>(i, j) = 1;
}
}
cv::Mat IDLP_channels[] = { cv::Mat_<float>(single),cv::Mat_<float>(single) };
cv::merge(IDLP_channels, 2, ihpf);
}
//创建巴特沃斯低通滤波器
void CreateButterworthLowpassFilter(SBLPFParam& param, cv::Mat& blpf)
{
cv::Mat single(param.sz.height, param.sz.width, CV_32F);
cv::Point centre = cv::Point(param.sz.height / 2, param.sz.width / 2);
double radius;
for (int i = 0; i < param.sz.height; i++)
{
for (int j = 0; j < param.sz.width; j++)
{
radius = sqrt(pow((float)(i - centre.x), 2) + pow((float)(j - centre.y), 2));
single.at<float>(i, j) = 1.0 / (1.0 + pow(radius / param.r, 2 * param.n));
}
}
cv::Mat IDLP_channels[] = { cv::Mat_<float>(single),cv::Mat_<float>(single) };
cv::merge(IDLP_channels, 2, blpf);
}
//创建巴特沃斯高通滤波器
void CreateButterworthHighpassFilter(SBHPFParam& param, cv::Mat& bhpf)
{
cv::Mat single(param.sz.height, param.sz.width, CV_32F);
cv::Point centre = cv::Point(param.sz.height / 2, param.sz.width / 2);
double radius;
for (int i = 0; i < param.sz.height; i++)
{
for (int j = 0; j < param.sz.width; j++)
{
radius = sqrt(pow((float)(i - centre.x), 2) + pow((float)(j - centre.y), 2));
single.at<float>(i, j) = 1.0 / (1.0 + pow( param.r/radius, 2 * param.n));
}
}
cv::Mat IDLP_channels[] = { cv::Mat_<float>(single),cv::Mat_<float>(single) };
cv::merge(IDLP_channels, 2, bhpf);
}
//创建高斯低通滤波器
void CreateGaussianLowpassFilter(SGLPFParam& param, cv::Mat& glpf)
{
cv::Mat single(param.sz.height, param.sz.width, CV_32F);
cv::Point centre = cv::Point(param.sz.height / 2, param.sz.width / 2);
double radius;
for (int i = 0; i < param.sz.height; i++)
{
for (int j = 0; j < param.sz.width; j++)
{
radius = sqrt(pow((float)(i - centre.x), 2) + pow((float)(j - centre.y), 2));
single.at<float>(i, j) = exp(-radius*radius/(2*param.r*param.r));
}
}
cv::Mat IDLP_channels[] = { cv::Mat_<float>(single),cv::Mat_<float>(single) };
cv::merge(IDLP_channels, 2, glpf);
}
//创建高斯高通滤波器
void CreateGaussianHighpassFilter(SGHPFParam& param, cv::Mat& ghpf)
{
cv::Mat single(param.sz.height, param.sz.width, CV_32F);
cv::Point centre = cv::Point(param.sz.height / 2, param.sz.width / 2);
double radius;
for (int i = 0; i < param.sz.height; i++)
{
for (int j = 0; j < param.sz.width; j++)
{
radius = sqrt(pow((float)(i - centre.x), 2) + pow((float)(j - centre.y), 2));
single.at<float>(i, j) = 1-exp(-radius*radius / (2 * param.r*param.r));
}
}
cv::Mat IDLP_channels[] = { cv::Mat_<float>(single),cv::Mat_<float>(single) };
cv::merge(IDLP_channels, 2, ghpf);
}
#include "../include/importOpenCV.h"
#include "../include/baseOps.h"
#include "../include/opencv400/opencv2/core.hpp"
#include <iostream>
int main()
{
//将工作目录设置到EXE所在的目录。
SetCurrentDirectoryToExePath();
cv::Mat src = cv::imread("../images/18.jpg");
cv::imshow("原图", src);
cv::Mat input;
cv::Mat imgHls;
std::vector<cv::Mat> vHls;
if (src.channels() == 3)
{
cvtColor(src, imgHls, cv::COLOR_BGR2HSV);
split(imgHls, vHls);
vHls[2].copyTo(input);
}
else
src.copyTo(input);
int w = cv::getOptimalDFTSize(input.cols);
int h = cv::getOptimalDFTSize(input.rows);
cv::Mat padded;
cv::copyMakeBorder(input, padded, 0, h - input.rows, 0, w - input.cols, cv::BORDER_CONSTANT, cv::Scalar::all(0));
padded.convertTo(padded, CV_32FC1);
for (int i = 0; i<padded.rows; i++)
{
float *ptr = padded.ptr<float>(i);
for (int j = 0; j<padded.cols; j++) ptr[j] *= pow(-1, i + j);
}
cv::Mat plane[] = { padded,cv::Mat::zeros(padded.size(),CV_32F) };
cv::Mat complexImg;
cv::merge(plane, 2, complexImg);
cv::dft(complexImg, complexImg);
//************************filter****************************
cv::Mat ILPF;
cv::Mat IHPF;
SILPFParam param;
param.sz = padded.size();
param.r = 50;
CreateIdealLowpassFilter(param, ILPF);
CreateIdealHighpassFilter(param, IHPF);
cv::multiply(complexImg, ILPF, ILPF);
cv::multiply(complexImg, IHPF, IHPF);
cv::Mat BLPF;
cv::Mat BHPF;
SBLPFParam bfparam;
bfparam.sz = padded.size();
bfparam.r = 50;
bfparam.n = 2;
CreateButterworthLowpassFilter(bfparam, BLPF);
CreateButterworthHighpassFilter(bfparam, BHPF);
cv::multiply(complexImg, BLPF, BLPF);
cv::multiply(complexImg, BHPF, BHPF);
cv::Mat GLPF;
cv::Mat GHPF;
SGLPFParam gfparam;
gfparam.sz = padded.size();
gfparam.r = 50;
CreateGaussianLowpassFilter(gfparam, GLPF);
CreateGaussianHighpassFilter(gfparam, GHPF);
cv::multiply(complexImg, GLPF, GLPF);
cv::multiply(complexImg, GHPF, GHPF);
ShowComplexImg("dft", complexImg);
ShowComplexImg("ILPF", ILPF);
ShowComplexImg("IHPF", IHPF);
ShowComplexImg("ILPF", BLPF);
ShowComplexImg("IHPF", BHPF);
ShowComplexImg("GLPF", GLPF);
ShowComplexImg("GHPF", GHPF);
//******************************************************************
//*************************idft*************************************
cv::idft(ILPF, ILPF);
split(ILPF, plane);
cv::magnitude(plane[0], plane[1], plane[0]);
cv::normalize(plane[0], plane[0], 1, 0, cv::NORM_MINMAX);
cv::Mat dst;
if (src.channels() == 3)
{
int nWidth = vHls[0].cols > plane[0].cols ? plane[0].cols : vHls[0].cols;
int nHeight = vHls[0].rows > plane[0].rows ? plane[0].rows : vHls[0].rows;
plane[0] = plane[0] * 255;
plane[0].convertTo(vHls[2], CV_8UC1);
cv::Rect rt(0, 0, nWidth, nHeight);
vHls[0] = vHls[0](rt);
vHls[1] = vHls[1](rt);
vHls[2] = vHls[2](rt);
merge(vHls, imgHls);
cvtColor(imgHls, dst, cv::COLOR_HSV2BGR);
}
else
plane[0].copyTo(dst);
cv::imshow("idft-ILPF",dst);
cv::idft(IHPF, IHPF);
cv::split(IHPF, plane);
cv::magnitude(plane[0], plane[1], plane[0]);
cv::normalize(plane[0], plane[0], 1, 0, cv::NORM_MINMAX);
if (src.channels() == 3)
{
plane[0] = plane[0] * 255;
plane[0].convertTo(vHls[2], CV_8UC1);
int nWidth = vHls[0].cols > plane[0].cols ? plane[0].cols : vHls[0].cols;
int nHeight = vHls[0].rows > plane[0].rows ? plane[0].rows : vHls[0].rows;
cv::Rect rt(0, 0, nWidth, nHeight);
vHls[0] = vHls[0](rt);
vHls[1] = vHls[1](rt);
vHls[2] = vHls[2](rt);
merge(vHls, imgHls);
cvtColor(imgHls, dst, cv::COLOR_HSV2BGR);
}
else
plane[0].copyTo(dst);
cv::imshow("idft-IHPF", dst);
cv::idft(BLPF, BLPF);
cv::split(BLPF, plane);
cv::magnitude(plane[0], plane[1], plane[0]);
cv::normalize(plane[0], plane[0], 1, 0, cv::NORM_MINMAX);
if (src.channels() == 3)
{
plane[0] = plane[0] * 255;
plane[0].convertTo(vHls[2], CV_8UC1);
int nWidth = vHls[0].cols > plane[0].cols ? plane[0].cols : vHls[0].cols;
int nHeight = vHls[0].rows > plane[0].rows ? plane[0].rows : vHls[0].rows;
cv::Rect rt(0, 0, nWidth, nHeight);
vHls[0] = vHls[0](rt);
vHls[1] = vHls[1](rt);
vHls[2] = vHls[2](rt);
merge(vHls, imgHls);
cvtColor(imgHls, dst, cv::COLOR_HSV2BGR);
}
else
plane[0].copyTo(dst);
cv::imshow("idft-BLPF", dst);
cv::idft(BHPF, BHPF);
cv::split(BHPF, plane);
cv::magnitude(plane[0], plane[1], plane[0]);
cv::normalize(plane[0], plane[0], 1, 0, cv::NORM_MINMAX);
if (src.channels() == 3)
{
plane[0] = plane[0] * 255;
plane[0].convertTo(vHls[2], CV_8UC1);
int nWidth = vHls[0].cols > plane[0].cols ? plane[0].cols : vHls[0].cols;
int nHeight = vHls[0].rows > plane[0].rows ? plane[0].rows : vHls[0].rows;
cv::Rect rt(0, 0, nWidth, nHeight);
vHls[0] = vHls[0](rt);
vHls[1] = vHls[1](rt);
vHls[2] = vHls[2](rt);
merge(vHls, imgHls);
cvtColor(imgHls, dst, cv::COLOR_HSV2BGR);
}
else
plane[0].copyTo(dst);
cv::imshow("idft-BHPF", dst);
cv::idft(GLPF, GLPF);
cv::split(GLPF, plane);
cv::magnitude(plane[0], plane[1], plane[0]);
cv::normalize(plane[0], plane[0], 1, 0, cv::NORM_MINMAX);
if (src.channels() == 3)
{
plane[0] = plane[0] * 255;
plane[0].convertTo(vHls[2], CV_8UC1);
int nWidth = vHls[0].cols > plane[0].cols ? plane[0].cols : vHls[0].cols;
int nHeight = vHls[0].rows > plane[0].rows ? plane[0].rows : vHls[0].rows;
cv::Rect rt(0, 0, nWidth, nHeight);
vHls[0] = vHls[0](rt);
vHls[1] = vHls[1](rt);
vHls[2] = vHls[2](rt);
merge(vHls, imgHls);
cvtColor(imgHls, dst, cv::COLOR_HSV2BGR);
}
else
plane[0].copyTo(dst);
cv::imshow("idft-GLPF", dst);
cv::idft(GHPF, GHPF);
cv::split(GHPF, plane);
cv::magnitude(plane[0], plane[1], plane[0]);
cv::normalize(plane[0], plane[0], 1, 0, cv::NORM_MINMAX);
if (src.channels() == 3)
{
plane[0] = plane[0] * 255;
plane[0].convertTo(vHls[2], CV_8UC1);
int nWidth = vHls[0].cols > plane[0].cols ? plane[0].cols : vHls[0].cols;
int nHeight = vHls[0].rows > plane[0].rows ? plane[0].rows : vHls[0].rows;
cv::Rect rt(0, 0, nWidth, nHeight);
vHls[0] = vHls[0](rt);
vHls[1] = vHls[1](rt);
vHls[2] = vHls[2](rt);
merge(vHls, imgHls);
cvtColor(imgHls, dst, cv::COLOR_HSV2BGR);
}
else
plane[0].copyTo(dst);
cv::imshow("idft-GHPF", dst);
cv::waitKey();
return 0;
}
原图及幅度谱
理想低通滤波器及理想高通滤波器
巴特沃斯低通滤波器及巴特沃斯高通滤波器
高斯低通滤波器及高斯高通滤波器