【图像处理】-010 图像频域处理
图像的频域处理能够在频率域内对图像进行滤波、重建、判断平移旋转等操作。这一篇博客主要用于记录我对图像频率域处理的学习历程，因此，这篇博客会进行持续更新。
文章目录
1 傅立叶变换的理论依据
通常情况下，我们通过像素下标 ( x , y ) (x,y) (x,y)访问图像的内容，进行各种处理，这是在空间域内的访问方式，那么怎么才能访问图像的频域信息呢。为了理解图像的频率域表示，需要首先理解图像从空间域进行频率域的工具——傅立叶变换（Fourier Transform）。
1.1 傅立叶级数
傅立叶变换是由法国数学家吉恩·巴普提斯特·约瑟夫·傅立叶（Jean Baptiste Joseph Fourier）于1768年在巴黎提出的，发表于《热分析理论》一书中。傅立叶指出：任何周期函数都可以表示成不同频率的正弦或余弦之和的形式，每个正弦项或余弦项乘以不同的系数（现在称该和为傅立叶级数）。
非周期函数(但该曲线下的面积是有限的)也可以用正弦和/或余弦乘以加权函数的积分来表示。在这种情况下的公式就是傅立叶变换。
用傅立叶级数或变换表示的函数特征完全可以通过傅立叶反变换来重建，而不会丢失任何信息。 它可以使我们工作于傅立叶域(频域)，而在返回到函数原有的域中的时候不会丢失任何信息。
1.2 基础概念
1.2.1 复数
复数
C
C
C的定义如下：
(1.2-1)
C
=
R
+
j
I
C = R+jI \tag{1.2-1}
C=R+jI(1.2-1)
其中，
R
R
R和
I
I
I是实数，
j
j
j是一个等于
−
1
-1
−1的平方根的虚数。
R
R
R表示复数的实部，
I
I
I表示复数的虚部。实数是虚部
I
=
0
I=0
I=0的复数。一个复数
C
C
C的共轭表示为
C
∗
C^{*}
C∗，其定义为：
(1.2-2)
C
∗
=
R
−
j
I
C^{*}=R-jI \tag{1.2-2}
C∗=R−jI(1.2-2)
在极坐标表示下复数为：
(1.2-3)
C
=
∣
C
∣
(
c
o
s
θ
+
j
s
i
n
θ
)
C=|C|(cos\theta + jsin\theta) \tag{1.2-3}
C=∣C∣(cosθ+jsinθ)(1.2-3)
其中，
∣
C
∣
=
R
2
+
I
2
|C|=\sqrt{R^2 + I^2}
∣C∣=R2+I2是复平面的圆点到点
(
R
,
I
)
(R,I)
(R,I)的向量的长度，
θ
\theta
θ是该向量与实轴的夹角。
使用欧拉公式：
(1.2-4)
e
j
θ
=
c
o
s
θ
+
j
s
i
n
θ
e^{j\theta}=cos\theta + jsin\theta \tag{1.2-4}
ejθ=cosθ+jsinθ(1.2-4)
其中，
e
=
2.71828
…
e=2.71828 \dots
e=2.71828…，可以给出复数的形式如下：
(1.2-5)
C
=
∣
C
∣
e
j
θ
C=|C|e^{j\theta} \tag{1.2-5}
C=∣C∣ejθ(1.2-5)
1.2.2 傅立叶级数
在上面我们提到过，任何周期函数都可以表示成不同频率的正弦或余弦之和的形式，每个正弦项或余弦项乘以不同的系数（现在称该和为傅立叶级数）。具有周期
T
T
T的连续变量
t
t
t的周期函数
f
(
t
)
f(t)
f(t)可以被描述为乘以适当系数的正弦和余弦和，这个和就是傅立叶级数，形式如下：
(1.2-6)
f
(
t
)
=
∑
n
=
−
∞
∞
c
n
e
j
2
π
n
T
t
f(t)=\sum_{n=-\infty}^{\infty}c_{n}e^{j\frac{2\pi n}{T}t} \tag{1.2-6}
f(t)=n=−∞∑∞cnejT2πnt(1.2-6)
其中，
(1.2-7)
c
n
=
1
T
∫
−
T
/
2
T
/
2
f
(
t
)
e
−
j
2
π
n
T
t
d
t
,
n
=
0
,
±
1
,
±
2
,
…
c_n = \frac{1}{T}\int_{-T/2}^{T/2}f(t)e^{-j\frac{2\pi n}{T}t}dt, n=0,\pm 1, \pm 2, \dots \tag{1.2-7}
cn=T1∫−T/2T/2f(t)e−jT2πntdt,n=0,±1,±2,…(1.2-7)
是系数。
1.2.3 冲激
连续变量
t
t
t在
t
=
0
t=0
t=0处的单位冲激表示为
δ
(
t
)
\delta(t)
δ(t)，其定义为：
(1.2-8)
δ
(
t
)
=
{
∞
,
t
=
0
0
,
t
̸
=
0
\delta (t)=\begin{cases} \infty,&t=0& \\ 0,&t\not=0& \end{cases} \tag{1.2-8}
δ(t)={∞,0,t=0t̸=0(1.2-8)
其中，该函数还必须满足：
(1.2-9)
∫
−
∞
∞
δ
(
t
)
d
t
=
1
\int_{-\infty}^{\infty}\delta(t)dt = 1 \tag{1.2-9}
∫−∞∞δ(t)dt=1(1.2-9)
取样特性
(1.2-10)
∫
−
∞
∞
f
(
t
)
δ
(
t
−
t
0
)
d
t
=
f
(
t
0
)
\int_{-\infty}^{\infty}f(t)\delta(t-t_0)dt=f(t_0) \tag{1.2-10}
∫−∞∞f(t)δ(t−t0)dt=f(t0)(1.2-10)
它在冲激位置
t
0
t_0
t0位置得到一个函数值。
领
x
x
x表示一个离散变量，单位离散冲激
δ
(
x
)
\delta(x)
δ(x)在离散系统中的作用与处理连续变量的冲激
δ
(
t
)
\delta(t)
δ(t)的作用相同。其定义如下：
(1.2-11)
δ
(
x
)
=
{
1
,
x
=
0
0
,
x
̸
=
0
\delta (x)=\begin{cases} 1,&x=0& \\ 0,&x\not=0& \end{cases} \tag{1.2-11}
δ(x)={1,0,x=0x̸=0(1.2-11)
且其满足：
(1.2-12)
∑
x
=
−
∞
∞
δ
(
x
)
=
1
\sum_{x=-\infty}^{\infty}\delta(x)=1 \tag{1.2-12}
x=−∞∑∞δ(x)=1(1.2-12)
在位置
x
=
x
0
x=x_0
x=x0处的离散冲激，
(1.2-13)
∑
x
=
−
∞
∞
f
(
x
)
δ
(
x
−
x
0
)
=
f
(
x
0
)
\sum_{x=-\infty}^{\infty}f(x)\delta(x-x_0)=f(x_0) \tag{1.2-13}
x=−∞∑∞f(x)δ(x−x0)=f(x0)(1.2-13)
与一维离散冲激类似，对于离散变量
x
x
x和
y
y
y，二维离散冲激定义为：
(1.2-14)
δ
(
x
,
y
)
=
{
1
,
x
=
y
=
0
0
,
其
他
\delta(x,y)=\begin{cases} 1,&x=y=0& \\ 0,&其他& \end{cases} \tag{1.2-14}
δ(x,y)={1,0,x=y=0其他(1.2-14)
其特定满足：
(1.2-15)
∑
x
=
−
∞
∞
∑
y
=
−
∞
∞
f
(
x
,
y
)
δ
(
x
−
x
0
,
y
−
y
0
)
=
f
(
x
0
,
y
0
)
\sum_{x=-\infty}^{\infty}\sum_{y=-\infty}^{\infty}f(x,y)\delta(x-x_0 , y-y_0)=f(x_0 , y_0) \tag{1.2-15}
x=−∞∑∞y=−∞∑∞f(x,y)δ(x−x0,y−y0)=f(x0,y0)(1.2-15)
离散冲激的取样特性在该冲激所在位置产生离散函数
f
(
x
,
y
)
f(x,y)
f(x,y)的值。
1.3 二维离散傅立叶变换
二维离散傅立叶变换(DFT):
(1.3-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.3-1}
F(u,v)=x=0∑M−1y=0∑N−1f(x,y)e−j2π(ux/M+vy/N)(1.3-1)
其中，
f
(
x
,
y
)
f(x,y)
f(x,y)是大小为
M
×
N
M \times N
M×N的数字图像。给出变换
F
(
u
,
v
)
F(u,v)
F(u,v)，可以使用傅立叶反变换（IDFT）得到
f
(
x
,
y
)
f(x,y)
f(x,y):
(1.3-2)
f
(
x
,
y
)
=
1
M
N
∑
u
=
0
M
−
1
∑
v
=
0
N
−
1
F
(
u
,
v
)
e
j
2
π
(
u
x
/
M
+
v
y
/
N
)
f(x,y)=\frac{1}{MN}\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}F(u,v)e^{j2\pi (ux/M+vy/N)} \tag{1.3-2}
f(x,y)=MN1u=0∑M−1v=0∑N−1F(u,v)ej2π(ux/M+vy/N)(1.3-2)
2 傅立叶变换的性质
2.1 空间和频率间隔的关系
对连续函数
f
(
t
,
z
)
f(t,z)
f(t,z)取样生成了一幅数字图像
f
(
x
,
y
)
f(x,y)
f(x,y)，它由分别在
t
t
t,
z
z
z方向所取的
M
×
N
M \times N
M×N个样本点组成，领
Δ
T
\Delta T
ΔT和
Δ
Z
\Delta Z
ΔZ表示样本间的间隔，那么响应的离散频率域变量间的间隔分别由
(2.1-1)
Δ
u
=
1
M
Δ
T
\Delta u=\frac{1}{M\Delta T} \tag{2.1-1}
Δu=MΔT1(2.1-1)
(2.1-2)
Δ
v
=
1
N
Δ
Z
\Delta v = \frac{1}{N \Delta Z} \tag{2.1-2}
Δv=NΔZ1(2.1-2)
给出。频率域样本间的间隔与空间样本间的间距和样本数成反比。
2.2 平移和旋转
傅立叶变换的平移特性满足：
(2.2-1)
f
(
x
,
y
)
e
j
2
π
(
u
0
x
/
M
+
v
0
y
/
N
)
⇔
F
(
u
−
u
0
,
v
−
v
0
)
f(x,y)e^{j2\pi (u_0 x/M+v_0 y/N)} \Leftrightarrow F(u-u_0 , v- v_0) \tag{2.2-1}
f(x,y)ej2π(u0x/M+v0y/N)⇔F(u−u0,v−v0)(2.2-1)
(2.2-2)
f
(
x
−
x
0
,
y
−
y
0
)
⇔
F
(
u
,
v
)
e
−
j
2
π
(
u
x
0
/
M
+
v
y
0
/
N
)
f(x-x_0 , y-y_0) \Leftrightarrow F(u,v)e^{-j2\pi (ux_0 /M+vy_0 /N)} \tag{2.2-2}
f(x−x0,y−y0)⇔F(u,v)e−j2π(ux0/M+vy0/N)(2.2-2)
用指数项乘以
f
(
x
,
y
)
f(x,y)
f(x,y)将使DFT的原点移到点
(
u
0
,
v
0
)
(u_0 ,v_0)
(u0,v0)，用负指数项乘以
F
(
u
,
v
)
F(u,v)
F(u,v)将使
f
(
x
,
y
)
f(x,y)
f(x,y)的原点移动到点
(
x
0
,
y
0
)
(x_0 ,y_0)
(x0,y0)，平移不改变
F
(
u
,
v
)
F(u,v)
F(u,v)的幅度（谱）。
若 f ( x , y ) f(x,y) f(x,y)旋转角度 θ 0 \theta _0 θ0角度，则 F ( u , v ) F(u,v) F(u,v)也旋转形同的角度，反之，若 F ( u , v ) F(u,v) F(u,v)旋转了一个角度， f ( x , y ) f(x,y) f(x,y)也旋转相同的角度。
clc;
clear;
close all;
% 原图
src = zeros(256,256);
src(97:160,97:160)=1;
rsrc = imrotate(src,30);
tsrc = imtranslate(src,[50,0]);
rtsrc = imtranslate(rsrc,[50,0]);
ssrc = zeros(256,256);
ssrc(65:192,65:192)=1;
src_fft = fft2(src);
src_p = atan2(imag(src_fft),real(src_fft));
src_fft=fftshift(src_fft);
src_fft=abs(src_fft);
src_fft=log(src_fft+1);
rsrc_fft = fft2(rsrc);
rsrc_p = atan2(imag(rsrc_fft),real(rsrc_fft));
rsrc_fft=fftshift(rsrc_fft);
rsrc_fft=abs(rsrc_fft);
rsrc_fft=log(rsrc_fft+1);
tsrc_fft = fft2(tsrc);
tsrc_p = atan2(imag(tsrc_fft),real(tsrc_fft));
tsrc_fft=fftshift(tsrc_fft);
tsrc_fft=abs(tsrc_fft);
tsrc_fft=log(tsrc_fft+1);
rtsrc_fft = fft2(rtsrc);
rtsrc_p = atan2(imag(rtsrc_fft),real(rtsrc_fft));
rtsrc_fft=fftshift(rtsrc_fft);
rtsrc_fft=abs(rtsrc_fft);
rtsrc_fft=log(rtsrc_fft+1);
ssrc_fft = fft2(ssrc);
ssrc_p = atan2(imag(ssrc_fft),real(ssrc_fft));
ssrc_fft=fftshift(ssrc_fft);
ssrc_fft=abs(ssrc_fft);
ssrc_fft=log(ssrc_fft+1);
figure(1),subplot(3,5,1),imshow(src),title('原图'),subplot(3,5,6),imshow(src_fft),title('原图FFT幅值'),subplot(3,5,11),imshow(src_p),title('原图FFT相位');
subplot(3,5,2),imshow(rsrc),title('原图旋转30度'),subplot(3,5,7),imshow(rsrc_fft),title('原图旋转30度FFT幅值'),subplot(3,5,12),imshow(rsrc_p),title('原图旋转30度FFT相位');
subplot(3,5,3),imshow(tsrc),title('原图平移50像素'),subplot(3,5,8),imshow(tsrc_fft),title('原图平移50像素FFT幅值'),subplot(3,5,13),imshow(tsrc_p),title('原图平移50像素FFT相位');
subplot(3,5,4),imshow(rtsrc),title('原图旋转30度平移50像素'),subplot(3,5,9),imshow(rtsrc_fft),title('原图旋转30度平移50像素FFT幅值'),subplot(3,5,14),imshow(rtsrc_p),title('原图旋转30度平移50像素FFT相位');
subplot(3,5,5),imshow(ssrc),title('原图放大'),subplot(3,5,10),imshow(ssrc_fft),title('原图放大FFT幅值'),subplot(3,5,15),imshow(ssrc_p),title('原图放大FFT相位');
从结果可以看出，
- 平移不改变图像的幅度谱，但会改变相位；
- 旋转会是幅度谱旋转相同角度，相位会被改变；
- 缩放会同时改变幅度谱和相位；
3 学习傅立叶变换的资料
- (1) https://www.cnblogs.com/h2zZhou/p/8405717.html
这一片博文写得很好。里面用很多生动的图，描述了傅立叶变换各频率正弦余弦分量叠加的过程。原文出处： 韩昊 作 者：韩 昊 知 乎：Heinrich 微 博：@花生油工人 知乎专栏：与时间无关的故事 谨以此文献给大连海事大学的吴楠老师，柳晓鸣老师，王新年老师以及张晶泊老师。 转载的同学请保留上面这句话，谢谢。如果还能保留文章来源就更感激不尽了。
4 频域处理的应用
- (1) 频域滤波处理
主要是在频域中进行低通和高通滤波。我将这部分整理成了一篇文章。
https://blog.csdn.net/freehawkzk/article/details/85278501 - (2) 同态滤波
同态滤波在处理光照不均匀的图像时效果不错。
https://blog.csdn.net/freehawkzk/article/details/85263006
8957
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
