- 数字图像处理（3）——频率域图像增强
- 1 背景知识
- 1.1 图像的频率
- 1.2 傅立叶变换
- 1.2.1 一维连续傅立叶变换
- 1.2.2 二维连续傅立叶变换
- 1.2.3 一维离散傅立叶变换
- 1.2.4 二维离散傅立叶变换
- 1.2.5 频谱图
- 2 频率域滤波
- 2.1 频率域平滑（低通）滤波器
- 2.1.1 理想低通滤波器（ILPF）
- 2.1.2 巴特沃思低通滤波器(BLPF)
- 2.1.3 高斯低通滤波器（GLPF）
- 2.2 频率域锐化（高通）滤波器
- 2.2.1 理想高通滤波器（IHPF）
- 2.2.2 巴特沃思高通滤波器（BHPF）
- 2.2.3 高斯高通滤波器（GHPF）
- 2.3 同态滤波
- 3 Matlab命令
- 4 嘘。TM
图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：风景照中大面积的陆地在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。
所以利用频率域研究图像，一些在空间域表述困难的增强任务，在频率域中变得非常普通；且滤波在频率域更直观，它可解释空间域滤波的某些性质。
1.2 傅立叶变换首先我们知道，周期函数可表示为不同频率的正弦和/或余弦和的形式，非周期函数可以用正弦和/或余弦乘以加权函数的积分来表示。此情况下的公式就是傅立叶变换。
傅里叶反变换：函数特征可以通过反变换来重建，不丢失任何信息。
Remark :记得宇宙第一欧拉公式。
1.2.1 一维连续傅立叶变换Lemma1.2.1 若一维实空间上的函数
f
(
x
)
f(x)
f(x)绝对可积，则函数
F
(
u
)
=
∫
−
∞
∞
f
(
x
)
e
−
i
2
π
u
x
d
x
F(u)=int_{-infty}^{infty}f(x)e^{-i2pi ux}dx
F(u)=∫−∞∞f(x)e−i2πuxdx
存在，若
F
(
u
)
F(u)
F(u)绝对可积，则
f
(
x
)
=
∫
−
∞
∞
F
(
u
)
e
i
2
π
u
x
d
u
f(x)=int_{-infty}^{infty}F(u)e^{i2pi ux}du
f(x)=∫−∞∞F(u)ei2πuxdu
**Def1.2.1 **记函数
f
(
x
)
f(x)
f(x)，其Fourier变换为
F
(
u
)
F(u)
F(u)，
f
(
x
)
f(x)
f(x)为
F
(
u
)
F(u)
F(u)的Fourier逆变换，记
F
(
u
)
=
F
(
f
(
x
)
)
,
f
(
x
)
=
F
−
1
(
F
(
u
)
)
F(u)=mathcal F(f(x)) , f(x)=mathcal F^{-1}(F(u))
F(u)=F(f(x)) , f(x)=F−1(F(u))
Lemma1.2.2 若二维实空间上的函数
f
(
x
,
y
)
f(x,y)
f(x,y)绝对可积，则函数
F
(
u
,
v
)
=
∫
R
2
f
(
x
,
y
)
e
−
i
2
π
(
u
x
+
v
y
)
d
x
d
y
F(u,v)=int_{mathbb R^2}f(x,y)e^{-i2pi (ux+vy)}dxdy
F(u,v)=∫R2f(x,y)e−i2π(ux+vy)dxdy
存在，若
F
(
u
,
v
)
F(u,v)
F(u,v)绝对可积，则
f
(
x
,
y
)
=
∫
R
2
F
(
u
,
v
)
e
i
2
π
(
u
x
+
v
y
)
d
u
d
v
f(x,y)=int_{mathbb R^2}F(u,v)e^{i2pi (ux+vy)}dudv
f(x,y)=∫R2F(u,v)ei2π(ux+vy)dudv
Def1.2.3 单变量离散函数
f
(
x
)
,
(
x
=
0
,
1
,
.
.
.
,
M
−
1
)
f(x),(x=0,1,...,M-1)
f(x),(x=0,1,...,M−1)的傅立叶变换
F
(
u
)
F(u)
F(u)定义为
F
(
u
)
=
1
M
∑
x
=
0
M
−
1
f
(
x
)
e
−
i
2
π
u
x
M
,
u
=
0
,
1
,
.
.
.
,
M
−
1
F(u)=frac1Msum_{x=0}^{M-1}f(x)e^frac{-i2pi ux}{M} , u=0,1,...,M-1
F(u)=M1x=0∑M−1f(x)eM−i2πux , u=0,1,...,M−1
给定
F
(
u
)
F(u)
F(u)，通过傅立叶逆变换可得到
f
(
x
)
f(x)
f(x)
f
(
x
)
=
∑
u
=
0
M
−
1
F
(
u
)
e
i
2
π
u
x
M
,
x
=
0
,
1
,
.
.
.
,
M
−
1
f(x)=sum_{u=0}^{M-1}F(u)e^frac{i2pi ux}M , x=0,1,...,M-1
f(x)=u=0∑M−1F(u)eMi2πux , x=0,1,...,M−1
Def1.2.4 二维图像尺寸为
M
×
N
Mtimes N
M×N的函数
f
(
x
,
y
)
f(x,y)
f(x,y)的离散傅立叶变换
F
(
u
,
v
)
F(u,v)
F(u,v)为
F
(
u
,
v
)
=
1
M
N
∑
x
=
0
M
−
1
∑
y
=
0
N
−
1
f
(
x
,
y
)
e
−
i
2
π
(
u
x
/
M
+
v
y
/
N
)
,
u
=
0
,
1
,
.
.
.
,
M
−
1
,
v
=
0
,
1
,
.
.
.
,
N
−
1
F(u,v)=frac1{MN}sum_{x=0}^{M-1}sum_{y=0}^{N-1}f(x,y)e^{-i2pi (ux/M+vy/N)} , u=0,1,...,M-1 , v=0,1,...,N-1
F(u,v)=MN1x=0∑M−1y=0∑N−1f(x,y)e−i2π(ux/M+vy/N) , u=0,1,...,M−1 , v=0,1,...,N−1
给定
F
(
u
,
v
)
F(u,v)
F(u,v)，通过傅立叶逆变换可得到
f
(
x
,
y
)
f(x,y)
f(x,y)
f
(
x
,
y
)
=
∑
u
=
0
M
−
1
∑
v
=
0
N
−
1
F
(
u
,
v
)
e
i
2
π
(
u
x
/
M
+
v
y
/
N
)
,
x
=
0
,
1
,
.
.
.
,
M
−
1
,
y
=
0
,
1
,
.
.
.
,
N
−
1
f(x,y)=sum_{u=0}^{M-1}sum_{v=0}^{N-1}F(u,v)e^{i2pi (ux/M+vy/N)} , x=0,1,...,M-1 , y=0,1,...,N-1
f(x,y)=u=0∑M−1v=0∑N−1F(u,v)ei2π(ux/M+vy/N) , x=0,1,...,M−1 , y=0,1,...,N−1
remark：u和v是频率变量，x和y是空间图像变量。
由Def1.2.4，通过宇宙第一欧拉公式转化为极坐标形式，由此可以得出幅度谱（频率谱），相位谱和功率谱。
二维图像的频谱图像就是将傅立叶变换后的二维图像作为亮度给显示出来，但由于图像的频谱幅度随频率增大而迅速衰减，故引入对数化的操作 D ( u , v ) = l o g ( 1 + K ∣ F ( u , v ) ∣ ) D(u,v)=log(1+K|F(u,v)|) D(u,v)=log(1+K∣F(u,v)∣)使其便于观察。
相关代码：
clear all
clc;
c=imread('c.tif');
imshow(c);
fftc=fft2(c); %二维离散傅立叶变换
%figure;
%imshow(fftc);
sfftc=fftshift(fftc); %直流分量移到频谱中心
R=real(sfftc); %取傅立叶变换的实部
I=imag(sfftc); %取傅立叶变换的虚部
A=sqrt(R.^2+I.^2); %计算频谱幅值
A=(A-min(min(A)))/(max(max(A-min(min(A)))))*225;%归一化
figure;
imshow(A);
在这里插入图片描述
那么从频谱图我们可以得出图像的能量分布。如果频谱图中暗点多，那么实际图像是比较柔和的（因为各点与邻域差异都不大，梯度相对较小），反之，如果频谱图中亮点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大的。还有就是明亮线和原始图像中对应的轮廓线是垂直的。如果原始图像中有圆形区域那么幅度谱中也呈圆形分布。
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)，H为滤波函数。
频率域的滤波步骤：
（1）图像预处理（中心变换）；
（2）计算（1）中图像DFT：F（u，v）；
（3）用滤波函数H（u，v）乘以F（u，v）；
（4）计算（3）中结果的IDFT；
（5）得到（4）中结果的实部；
（6）去中心化；
2.1 频率域平滑（低通）滤波器已知边缘和噪声等尖锐变化处于傅里叶变换的高频部分，那么平滑可以通过衰减高频成分的范围来实现。频率域平滑（低通）滤波器主要类别有：理想低通滤波器，巴特沃思低通滤波器，高斯低通滤波器。
2.1.1 理想低通滤波器（ILPF）截断傅里叶变换中的所有高频成分，这些高频成分处于指定距离D0（截止频率）之外
H
(
u
,
v
)
=
{
1
,
D
(
u
,
v
)
≤
D
0
0
,
D
(
u
,
v
)
>
D
0
H(u,v)=begin {cases} 1 , D(u,v)le D_0\ 0 , D(u,v)> D_0 end{cases}
H(u,v)={1 , D(u,v)≤D00 , D(u,v)>D0
频率矩形的中心在(u,v)=(M/2,N/2)，从点(u,v)到中心（原点）的距离如下
D
(
u
,
v
)
=
[
(
u
−
M
2
)
2
+
(
v
−
N
2
)
2
]
1
2
D(u,v)=[(u-frac M2)^2+(v-frac N2)^2]^{frac12}
D(u,v)=[(u−2M)2+(v−2N)2]21
说明：在半径为D0的圆内，所有频率没有衰减地通过滤波器，而在此半径的圆之外的所有频率完全被衰减掉。半径D0越小，模糊越大；半径D0越大，模糊越小。
相关代码
a=imread('1.jpg');
b=rgb2gray(a);
b=double(b);
f=fft2(b);
g=fftshift(f);
[M,N]=size(f);
d0=15;
n1=floor(M/2);
n2=floor(N/2);
for i=1:M
for j=1:N
d=sqrt((i-n1)^2+(j-n2)^2);
if d<=d0 h=1;
else h=0;
end
g(i,j)=h*g(i,j);
end
end
g=ifftshift(g);
g=uint8(real(ifft2(g)));
imshow(g,[0,255])
2.1.2 巴特沃思低通滤波器(BLPF)
n级巴特沃思低通滤波器定义如下
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
D0为截至频率距原点的距离，D(u,v)是点(u,v)距原点的距离。
不同于ILPF，BLPF变换函数在通带与被滤除的频率之间没有明显的截断。
应用：可用于平滑处理，如图像由于量化不足产生虚假轮廓时，常可用低通滤波进行平滑以改进图像质量。通常，BLPF的平滑效果好于ILPF（振铃现象）。
2.1.3 高斯低通滤波器（GLPF）二维高斯低通滤波器定义如下
H
(
u
,
v
)
=
e
−
D
(
u
,
v
)
2
/
2
σ
2
H(u,v)=e^{-D(u,v)^2/2sigma^2}
H(u,v)=e−D(u,v)2/2σ2
D(u,v)是点(u,v)距原点的距离。
结论:
• GLPF不能达到有相同截止频率的二阶BLPF的平滑效果
• GLPF没有振铃
• 如果需要严格控制低频和高频之间截至频率的过渡，选用BLPF，代价是可能产生振铃
2.2 频率域锐化（高通）滤波器对比频率域平滑（低通）滤波器，不难发现高通滤波器与其三还是有很大联系，在这就不过多介绍了。
截断傅里叶变换中的所有低频成分，这些低频成分处于指定距离D0（截止频率）之内
H
(
u
,
v
)
=
{
0
,
D
(
u
,
v
)
≤
D
0
1
,
D
(
u
,
v
)
>
D
0
H(u,v)=begin{cases} 0 , D(u,v)le D_0\ 1 , D(u,v)>D_0 end{cases}
H(u,v)={0 , D(u,v)≤D01 , D(u,v)>D0
频率矩形的中心在(u,v)=(M/2,N/2)，从点(u,v)到中心（原点）的距离如下
D
(
u
,
v
)
=
[
(
u
−
M
2
)
2
+
(
v
−
N
2
)
2
]
1
2
D(u,v)=[(u-frac M2)^2+(v-frac N2)^2]^{frac12}
D(u,v)=[(u−2M)2+(v−2N)2]21
n级且截止频率为距原点距离为D0的巴特沃思高通滤波器定义如下
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
截距距原点为D0的高斯高通滤波器定义为
H
(
u
,
v
)
=
1
−
e
−
D
(
u
,
v
)
2
/
2
D
0
2
H(u,v)=1-e^{-D(u,v)^2/2D_0^2}
H(u,v)=1−e−D(u,v)2/2D02
一般图像中的大部分能量集中在低频分量里,高通滤波会将很多低频分量(特别是直流分量)滤除, 导致增强图像中边缘得到加强, 但光滑区域灰度减弱变暗甚至接近黑色。
频域滤波可以灵活地解决加性噪声问题，但无法消减乘性或卷积性噪声。
同态滤波是一种在频域中同时将图像亮度范围进行压缩和将图像对比度进行增强的方法，是基于图像成像模型进行的。一幅图 f(x,y)可以表示成照度分量 i(x,y)与反射分量 r(x,y)的乘积， f ( x , y ) = i ( x , y ) r ( x , y ) f(x,y)=i(x,y)r(x,y) f(x,y)=i(x,y)r(x,y)。同态滤波的基本思想是将非线性问题转化成线性问题处理,即先对非线性混杂信号作某种数学运算,变换成加性的,然后用线性滤波方法处理,最后作反变换运算,恢复处理后图像。
同态滤波特点：能消除乘性噪声，能同时压缩图像的整体动态范围，并增加图像中相邻区域间的对比度。
空域滤波与频域滤波的区别：
◆空域基于模板操作，每次只是基于部分像素的性质
◆频域利用图像中所有像素的数据，具有全局性质，更好地体现图像的整体特征
3 Matlab命令| 函数名 | 参数说明 | 作用 |
|---|---|---|
| Y = fft2(X) | X - Input array(matrix) Y - Frequency domain representation(vector | matrix ) | 二维快速傅立叶变换 |
| Y = fftshift(X) | X - Input array（ vector｜matrix） | 中心化 |
| X = ifftshift(Y) | Y - Input array（ vector｜matrix） | 逆中心化变换（重排数据序列，使得原点位于频率矩阵的左上角） |
| X = ifft2(Y) | X - Input array（matrix） | 二维快速傅立叶逆变换 |
滤波函数的程序就自己写一写Function咯。
Img = imread('1.jpg');
Img = fftshift(fft2(im2double(rgb2gray(Img))));
[row,col] = size(Img);
d0 = 10; %截止频率为10，数值越小越平滑
center_row = fix(row/2); center_col = fix(col/2);
for i = 1:row %d0为10的高斯低/高通滤波器
for j = 1:col
d = sqrt((i - center_row)^2 + (j - center_col)^2);
hl(i,j) = exp(-d^2 / (2 * d0^2)); %高斯低通滤波器
hh(i,j) = 1 - hl(i,j); %高斯高通滤波器
gl(i,j) = hl(i,j) * Img(i,j);%高斯低通滤波
gh(i,j) = hh(i,j) * Img(i,j);%高斯高通滤波
end
end
Kl = im2uint8(real(ifft2(ifftshift(gl))));
Kh = im2uint8(real(ifft2(ifftshift(gh))));
4 嘘。TM
（1）傅立叶变换得到的频谱图中，低频系数对应（ ）
答：变化平缓的部分。
（2）图像锐化处理算法包括（ ）
答：空域中利用一阶微分算子和二阶微分算子，频域中利用高通滤波器。
