数字图像处理MATLAB学习笔记(二)
Filtering in the Frequency Domain频域处理
主要通过Matlab语言学习傅立叶变换下的频域处理和实践
1. The 2-D Discrete Fourier Transform(2-D DFT)二维离散傅里叶变换
1.1 DFT和IDFT
其中满足的条件:
- x = 0 , 1 , 2 , … , M − 1 {x=0,1,2,\dots,M-1} x=0,1,2,…,M−1
- y = 0 , 1 , 2 , … , N − 1 {y=0,1,2,\dots,N-1} y=0,1,2,…,N−1
- u = 0 , 1 , 2 , … , M − 1 {u=0,1,2,\dots,M-1} u=0,1,2,…,M−1
- v = 0 , 1 , 2 , … , N − 1 {v=0,1,2,\dots,N-1} v=0,1,2,…,N−1
令
f
(
x
,
y
)
{f(x,y)}
f(x,y)表示一副
M
×
N
{M\times N}
M×N大小的图像,**二维离散傅里叶变换(2-D DFT)**记为
F
(
u
,
v
)
{F(u,v)}
F(u,v),则有
F
(
u
,
v
)
=
∑
x
=
0
M
−
1
∑
y
=
0
N
−
1
f
(
x
,
y
)
e
x
p
(
−
j
2
π
(
u
x
M
+
v
y
N
)
)
F(u,v)=\sum^{M-1}_{x=0}\sum_{y=0}^{N-1}f(x,y)exp(-j2\pi(\frac{ux}{M}+\frac{vy}{N}))
F(u,v)=x=0∑M−1y=0∑N−1f(x,y)exp(−j2π(Mux+Nvy))
由式得,
F
(
0
,
0
)
=
m
e
a
n
(
s
u
m
(
f
(
x
,
y
)
)
)
{F(0,0)=mean(sum(f(x,y)))}
F(0,0)=mean(sum(f(x,y)))
**离散傅立叶反变换(Inverse DFT=IDFT)**表达式为:
f
(
x
,
y
)
=
1
M
N
∑
u
=
0
M
−
1
∑
v
=
0
N
−
1
F
(
u
,
v
)
e
x
p
(
j
2
π
(
u
x
M
+
v
y
N
)
)
f(x,y)=\frac{1}{MN}\sum^{M-1}_{u=0}\sum_{v=0}^{N-1}F(u,v)exp(j2\pi(\frac{ux}{M}+\frac{vy}{N}))
f(x,y)=MN1u=0∑M−1v=0∑N−1F(u,v)exp(j2π(Mux+Nvy))
因此,当得知
F
(
u
,
v
)
{F(u,v)}
F(u,v)时,我们可以通过计算IDFT的平均值得到
f
(
x
,
y
)
{f(x,y)}
f(x,y)。
1.2 DFT性质
尽管 f ( x , y ) {f(x,y)} f(x,y)是实函数,DFT通常是一个复数变换。
先设两个函数:
- R ( u , v ) {R(u,v)} R(u,v):the real components of F ( u , v ) {F(u,v)} F(u,v);实部
- I ( u , v ) I(u,v) I(u,v):the imaginary components of F ( u , v ) {F(u,v)} F(u,v);虚部
**傅立叶频谱(Fourier spectrum)**定义为:
∣
F
(
u
,
v
)
∣
=
[
R
2
(
u
,
v
)
+
I
2
(
u
,
v
)
]
1
2
|F(u,v)|=[R^2(u,v)+I^2(u,v)]^\frac{1}{2}
∣F(u,v)∣=[R2(u,v)+I2(u,v)]21
**变换相角(phase angle of transform)**定义为:
ϕ
(
u
,
v
)
=
a
c
r
t
a
n
[
I
(
u
,
v
)
R
(
u
,
v
)
]
\phi(u,v)=acrtan[\frac{I(u,v)}{R(u,v)}]
ϕ(u,v)=acrtan[R(u,v)I(u,v)]
极坐标(polar form)下表示复数(complex function)
F
(
u
,
v
)
{F(u,v)}
F(u,v)为
F
(
u
,
v
)
=
∣
F
(
u
,
v
)
∣
e
x
p
(
j
ϕ
(
u
,
v
)
)
F(u,v)=|F(u,v)|exp(j\phi(u,v))
F(u,v)=∣F(u,v)∣exp(jϕ(u,v))
**功率谱(Power spectrum)可定义为幅度(magnitude)的平方:
P
(
u
,
v
)
=
∣
F
(
u
,
v
)
∣
2
=
R
2
(
u
,
v
)
+
I
2
(
u
,
v
)
P(u,v)=|F(u,v)|^2=R^2(u,v)+I^2(u,v)
P(u,v)=∣F(u,v)∣2=R2(u,v)+I2(u,v)
如果
f
(
x
,
y
)
{f(x,y)}
f(x,y)是real function,傅里叶变换是关于远点共轭对称的,即
F
(
u
,
v
)
=
F
∗
(
−
u
,
−
v
)
F(u,v)=F^*(-u,-v)
F(u,v)=F∗(−u,−v)
即傅里叶谱频谱关于原点对称,即
∣
F
(
u
,
v
)
∣
=
∣
F
(
−
u
,
−
v
)
∣
|F(u,v)|=|F(-u, -v)|
∣F(u,v)∣=∣F(−u,−v)∣
同时满足下列equation(
k
1
{k_1}
k1和
k
2
{k_2}
k2均为Integer):
F
(
u
,
v
)
=
F
(
u
+
k
1
M
,
v
)
+
F
(
u
,
v
+
K
2
N
)
=
F
(
u
+
k
1
M
,
v
+
k
2
N
)
F(u,v)=F(u+k_1M,v)+F(u,v+K_2N)=F(u+k_1M,v+k_2N)
F(u,v)=F(u+k1M,v)+F(u,v+K2N)=F(u+k1M,v+k2N)
也就是说,DFT在u、v方向上是无穷周期(infinitely periodic)**的,周期由M、N决定。DFT的周期性表达为
f
(
x
,
y
)
=
f
(
x
+
k
1
M
,
y
)
+
f
(
x
,
y
+
K
2
N
)
=
f
(
x
+
k
1
M
,
y
+
k
2
N
)
f(x,y)=f(x+k_1M,y)+f(x,y+K_2N)=f(x+k_1M,y+k_2N)
f(x,y)=f(x+k1M,y)+f(x,y+K2N)=f(x+k1M,y+k2N)
即,IDFT获取到的图像也是具有无穷周期的。所以,我们计算时只计算一个周期内的,即
M
×
N
{M\times N}
M×N大小的图片。
2. Computing and Visualizing the 2-D DFT in MATLAB
DFT和IDFT实现:这里有一个fast Fourier transform(FFT) algorithm,在matlab中使用函数fft2。逆变换则使用函数ifft2。
F = fft2(f)
如果需要对图像进行0填充,则语法变为
F = fft2(f, P, Q)
输出F的大小变为PxQ.
Fourier spectrum实现:使用函数abs
S = abs(F)
abs会计算F中每一个元素的magnitude(real和imaginary components的平方和的平方根)
>> f = imread('Fig0303(a).tif');
>> imshow(f);xlabel('Image');
>> F = fft2(f);
>> subplot(1,2,1);
>> imshow(f);xlabel('Image');
>> subplot(1,2,2);
>> imshow(S, []);xlabel('Fourier spectrum');
![image-20211111232917979](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211111232917979.png)
周期性处理函数:fftshift。也就是将变换的原点移动到频域矩阵的中心,语法为
Fc = fftshift(F)
Fc就是中心变换后的结果
- 如果F是一个vector,Fc就是F中前一半与后一半向量的互换
- 如果F是一个Matrix,Fc就是F横竖分四份(奇数时前半部分向上取整),对角互换的结果
- 如果F是一个multidimensional array,Fc就是F
这里还可以进行逆操作,函数为ifftshift
根据上一节总结的, F ( u , v ) {F(u,v)} F(u,v)的real和imaginary components使用如下函数表示
- Real:function real(F)
- Imaginary:function imag(F)
phase angle of transform实现,使用atan2函数,传入real components和imaginary components
phi = atan2(I, R) = atan2(imag(F),real(F))
或者使用angle函数,直接传入F,不用再转换为real和imaginary
phi = angle(F)
此时,得到DFT有:
F = S.*exp(i*phi)
Note:这里的i指的是复数,即 i = 0 + 1 × i {i=0+1\times i} i=0+1×i
![image-20211111232827221](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211111232827221.png)
3. Filtering in the Frequency Domain频域滤波
3.1 Fundamentals基本知识
空间和频域上的线性滤波都是基于Convolution操作的。可以用公式表达为
f
(
x
,
y
)
⊙
h
(
x
,
y
)
⇔
H
(
u
,
v
)
F
(
u
,
v
)
f
(
x
,
y
)
h
(
x
,
y
)
⇔
H
(
u
,
v
)
⊙
F
(
u
,
v
)
f(x,y)\odot h(x,y)\Leftrightarrow H(u,v)F(u,v) \\ f(x,y)h(x,y)\Leftrightarrow H(u,v)\odot F(u,v)
f(x,y)⊙h(x,y)⇔H(u,v)F(u,v)f(x,y)h(x,y)⇔H(u,v)⊙F(u,v)
频域滤波的main idea就是选择一个filter transfer function用特定的规则修改F(u,v)。
函数解释:
- h(x,y):filtering mask卷积核
- H(u,v):filter transfer function滤波传递函数。
在处理discrete quantities时,F和H都是周期循环的,所以离散频域上的卷积也是周期性的,即DFT的卷积操作也称为circular convolution循环卷积。为了得到相同大小结果,我们需要进行0-padding进行填充。
- 原因:原因比较复杂,这个看了两遍才看懂。简单来讲,就是因为在进行周期性卷积操作时,如果函数非零部分的延续周期很靠近,周期性卷积就会造成相邻函数的重叠。如果进行0-padding,则可以避开这些重叠的部分,因为值是0,对结果无影响。
解决办法:构造两个padded function,每个输出大小都是 P × Q {P\times Q} P×Q。假设 f ( x , y ) {f(x,y)} f(x,y)和 h ( x , y ) {h(x,y)} h(x,y)的大小分别为 A × B {A\times B} A×B和 C × D {C\times D} C×D。
则,输出大小需满足的条件为:
P
≥
A
+
C
−
1
Q
≥
B
+
D
−
1
P \ge A+C-1\\ \\ Q \ge B+D-1
P≥A+C−1Q≥B+D−1
通过函数实现:function paddedsize
function PQ = paddedsize(AB, CD, PARAM)
%PADDEDSIZE Computes padded size useful for FFT-based filtering.
% PQ = PADDEDSIZE(AB), where AB is a two-element size vector,
% computes the two-element size vector PQ = 2*AB.
%
% PQ = PADDEDSIZE(AB, 'PWR2') computes the vector PQ such that
% PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX(AB).
%
% PQ = PADDEDSIZE(AB, CD), where AB and CD are two-element size
% vectors, computes the two-element size vector PQ. The elements
% of PQ are the smallest even integers greater than or equal to
% AB + CD - 1.
%
% PQ = PADDEDSIZE(AB, CD, 'PWR2') computes the vector PQ such that
% PQ(1) = PQ(2) = 2^nextpow2(2*m), where m is MAX([AB CD]).
if nargin == 1
PQ = 2*AB;
elseif nargin == 2 && ~ischar(CD)
PQ = AB + CD - 1;
PQ = 2 * ceil(PQ / 2);
elseif nargin == 2
m = max(AB); % Maximum dimension
% Find power-of-2 at least twice m.
P = 2^nextpow2(2*m);
PQ = [P, P];
elseif (nargin == 3) && strcmpi(PARAM, 'pwr2')
m = max([AB CD]); % Maximum dimension
P = 2^nextpow2(2*m);
PQ = [P, P];
else
error('Wrong number of inputs.');
end
这里还要自定义两个函数:lpfilter和dfuv
-
function lpfilter:得到指定类型的低通滤波器
-
function dftuv:计算网格频率矩阵U和V
我们通过调用上述函数,获得对比图为:
>> f = imread('Fig0305(a).tif');
>> imshow(f);
>> [M, N] = size(f);
>> [f, revertClass] = tofloat(f);
>> F = fft2(f);
>> sig = 10;
>> H = lpfilter('gaussian', M, N, sig);
>> G = H.*F;
>> g = ifft2(G);
>> g = revertClass(g);
>> imshow(g);
>> PQ = paddedsize(size(f)); % f is floating point.
>> Fp = fft2(f, PQ(1), PQ(2)); % computing the FFT with padding
>> Hp = lpfilter('gaussian', PQ(1), PQ(2), 2*sig);
>> Gp = Hp.*Fp;
>> gp = ifft2(Gp);
>> gpc = gp(1:size(f,1), 1:size(f,2));
>> gpc = revertClass(gpc);
>> imshow(gpc);
>> title('Image lowpass-filtered with padding');
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-c5PXUya6-1647138467506)(/Users/wanghe/Library/Application Support/typora-user-images/image-20211112133152273.png)]
3.2 Basic Steps in DFT Filtering DFT滤波基本操作
-
image数据类型转换:使用tofloat函数转换为浮点型数据。
-
获得图片padding扩充大小:使用paddedsie函数获得,传入参数为f的size。
-
傅里叶变换:需要加上padding值,得到F。
-
生成一个大小为PQ(2)xPQ(2)的H滤波函数。
-
使用滤波器H乘FFT交换得到G。
-
获得G的逆FFT交换得到g。
-
裁剪左上角位置图像,也就是原图大小。
-
滤波图像数据类型转换为原类型即可。
![image-20211112141827233](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211112141827233.png)
3.3 An M-function for Filtering in the Frequency Domain
设计一个函数,使得可以自动地进行频域滤波操作。函数名为dftfilt。
function g = dftfilt(f, H, classout)
%DFTFILT Performs frequency domain filtering
% g = DFTFILT(f, H, CLASSOUT) filters f in the frequency domain
% using the filter transfer function H. The output, g, is the
% filtered image, which has the same size as f.
%
% Valid values of CLASSOUT are
%
% 'original' The output is of the same class as the input.
% This is the default if CLASSOUT is not included
% in the call
% 'fltpoint' The output is floating point of class single, unless
% both f and H are of class double, in which case the
% output also is of class double.
%
% DFTFILE automatically pads f to be the same size as H. Both f
% and H must be real. In addition, H must be an uncentered,
% circularly-symmetric filter function.
% Convert the input to floating point.
[f, revertClass] = tofloat(f);
% Obtain the FFT of the padded input.
F = fft2(f, size(H, 1), size(H, 2));
% Perform filtering
g = ifft2(H.*F);
% Crop to original size.
g = g(1:size(f, 1), 1:size(f, 2)); % g is of class single here.
% Convert the output to the same class as the input if so specified.
if nargin == 2 || strcmp(classout, 'original')
g = revertClass(g);
elseif strcmp(classout, 'fltpoint')
return
else
error('Undefined class for the output image.');
end
4 Obtaining Frequency Domain Filters from Spatial Filters
function freqz2:计算FIR(有限脉冲响应)滤波的响应频率。是一个Linear Filter。
H = freqz2(h, R, C)
参数信息:
- h:二维空间滤波器
- H:二维频域滤波器
- R, C:通常传入padding值,即PQ(1)和PQ(2)
Example:空间与频域滤波器对比
图片大小为600x600
>> f = imread('Fig0309(a).tif');
>> f = tofloat(f);
>> F = fft2(f);
>> S = fftshift(log(1 + abs(F)));
>> imshow(S, []);
![image-20211112231344837](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211112231344837.png)
本例中,我们使用Sobel空间滤波器来增强垂直边缘。
>> h = fspecial('sobel')' % 符号'表示矩阵转置
h =
1 0 -1
2 0 -2
1 0 -1
然后,我们比较sobel滤波器下空间滤波和频域滤波的区别。
>> freqz2(h)
>> PQ = paddedsize(size(f)); % 生成H滤波器
>> H = freqz2(h, PQ(1), PQ(2));
>> H1 = ifftshift(H);
>> imshow(abs(H), []);
>> imshow(abs(H1), []);
![image-20211112233039893](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211112233039893.png)
接下来生成滤波图像
>> gs = imfilter(f, h); % spatial domain
>> gf = dftfilt(f, H1); % obtained in frequency domain
>> imshow(gs, []);
>> imshow(gf, []);
生成的图片偏灰色调,主要是gs和gf中存在负数。由于sobel模版适用于相应绝对值检测垂直边缘,所以图像的绝对值更具有相关性。
>> imshow(abs(gs), []);
>> imshow(abs(gf), []);
![image-20211112234314312](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211112234314312.png)
为了边缘显示更清晰,可以采用创建阈值化的二值图像。
>> imshow(abs(gs) > 0.2*abs(max(gs(:)))); % 选用0.2只显示强度比gs最大值的20%还要大的边缘
>> imshow(abs(gf) > 0.2*abs(max(gs(:)))); % 选用0.2只显示强度比gf最大值的20%还要大的边缘
![image-20211112234716260](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211112234716260.png)
通过计算gs与gf的差距,发现最大的差距仅为1.4063e-06,最小差别为0。差距非常小,可以忽略。
5 Generating Filters Directly in the Frequency Domain
5.1 Creating Meshgrid Arrays for Use in Implementing Filters in the Frequency Domain 建立用于实现频域滤波器的网格数组
这里需要创建的函数为:dftuv(M, N)
因为FFT计算假设变换的原点位置在左上角,距离计算也是根据该点进行。
function [U, V] = dftuv(M, N)
% DFTUV Computes meshgrid frequency matrices.
% [U, V] = DFTUV(M, N) computes meshgrid frequency matrices U and V. U and
% V are useful for computing frequency-domain filter functions that can be
% used with DFTFILT. U and V are both M-by-N.
% more details to see the textbook Page 93
%
% [U,V] = DFTUV(M,N)计算网格频率矩阵U和V。 U和V对于计算可与DFTFILT一起使用的
% 频域滤波器函数很有用。 U和V都是M-by-N。
% Set up range of variables.
% 设置变量范围
u = 0 : (M - 1);
v = 0 : (N - 1);
% Compute the indices for use in meshgrid.
% 计算网格的索引,即将网络的原点转移到左上角,因为FFT计算时变换的原点在左上角。
idx = find(u > M / 2);
u(idx) = u(idx) - M;
idy = find(v > N / 2);
v(idy) = v(idy) - N;
% Compute the meshgrid arrays.
% 计算网格矩阵
[V, U] = meshgrid(v, u);
5.2 Lowpass(Smoothing) Frequency Domain FIlters
一个理想的低通滤波器(ILPF)具有如下的传递方程
H
(
u
,
v
)
=
{
1
,
if
D
(
u
,
v
)
≤
D
0
0
,
if
D
(
u
,
v
)
>
D
0
H(u,v)=\begin{cases} 1, \text{if} \space D(u,v) \le D_0 \\ 0, \text{if} \space D(u,v) \gt D_0 \end{cases}
H(u,v)={1,if D(u,v)≤D00,if D(u,v)>D0
参数信息:
- D ( u , v ) {D(u,v)} D(u,v):点(u,v)到滤波器中心的距离
- D 0 {D_0} D0:正数
- 满足 D ( u , v ) = D 0 {D(u,v)=D_0} D(u,v)=D0的轨迹为圆
因为滤波器H要乘一幅图像的Fourier transform,我们看到理想滤波器切断(0)了园外的F(u,v)的分量,保留圆上和圆内的点(1)。ILPF通常用于解释振铃(ringing)和重叠误差(wraparound error)等现象。
Butterworth lowpass filter(BLFT):布特沃斯低通滤波器,具有从滤波中心到
D
0
{D_0}
D0的截止频率,表达式为
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
这里默认n=1。当
D
(
u
,
v
)
=
D
0
{D(u,v)=D_0}
D(u,v)=D0,滤波器降到最大值1的0.5,即
H
(
u
,
v
)
=
0.5
{H(u,v)=0.5}
H(u,v)=0.5。
Gaussian Lowpass Filter(GLPF):高斯低通滤波器,用的最多的一个
H
(
u
,
v
)
=
e
x
p
(
−
D
2
(
u
,
v
)
2
σ
2
)
H(u,v)=exp(\frac{-D^2(u,v)}{2\sigma^2})
H(u,v)=exp(2σ2−D2(u,v))
参数:
σ
=
D
0
{\sigma=D_0}
σ=D0,即标准差
根据截止参数,我们可以得到
H
(
u
,
v
)
=
e
x
p
(
−
D
2
(
u
,
v
)
2
D
0
2
)
H(u,v)=exp(\frac{-D^2(u,v)}{2D_0^2})
H(u,v)=exp(2D02−D2(u,v))
当
D
(
u
,
v
)
=
D
0
{D(u,v)=D_0}
D(u,v)=D0,滤波器降到最大值1的0.607
Example:Lowpass FIlter Test
>> f = imread('Fig0313(a).tif');
>> [f, revertClass] = tofloat(f);
>> PQ = paddedsize(size(f));
>> [U, V] = dftuv(PQ(1), PQ(2));
>> D = hypot(U, V); % Compute faster than D=sqrt(U.^2+V.^2)
>> D0 = 0.05*PQ(2);
>> F = fft2(f, PQ(1), PQ(2)); % Needed for the spectrum.
>> H = exp(-(D.^2)/(2*(D0.^2))); % GLPF
>> subplot(2,2,1);
>> imshow(f);title('image');
>> subplot(2,2,2);
>> imshow(fftshift(H));
>> title('GLPF shows as an Image');
>> subplot(2,2,3);
>> imshow(log(1+abs(fftshift(F))), []);
>> title("Spectrum of original image");
>> subplot(2,2,4);
>> imshow(g);title('Filtered image');
![image-20211113104813904](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211113104813904.png)
我们也可以将上述三种低通滤波器合并至一个M-Function文件中,函数为lpfilter
function H = lpfilter(type, M, N, D0, n)
% LPFILTER Computes frequency domain lowpass filters
% H = LPFILTER(TYPE, M, N, D0, n) creates the transfer function of a
% lowpass filter, H, of the specified TYPE and size (M-by-N). To view the
% filter as an image or mesh plot, it should be centered using H =
% fftshift(H)
% Valid value for TYPE, D0, and n are:
% 'ideal' Ideal lowpass filter with cutoff frequency D0. n need not be
% supplied. D0 must be positive.
% 'btw' Butterworth lowpass filter of order n, and cutoff D0. The
% default value for n is 1.0. D0 must be positive.
% 'gaussian'Gaussian lowpass filter with cutoff (standard deviation) D0.
% n need not be supplied. D0 must be positive.
%
% 得到指定类型的低通滤波器
% Use function dftuv to set up the meshgrid arrays needed for computing the
% required distances.
[U, V] = dftuv(M, N);
% Compute the distances D(U, V)
D = sqrt(U.^2 + V.^2);
% Begin filter computations
switch type
case 'ideal'
H = double(D <= D0);
case 'btw'
if nargin == 4
n = 1;
end
H = 1 ./ (1 + (D ./ D0) .^ (2 * n));
case 'gaussian'
H = exp(-(D .^ 2) ./ (2 * (D0 ^ 2)));
otherwise
error('Unkown filter type.')
end
5.3 Wireframe and Surface Plotting线框图和表面图的绘制
主要用到的函数:
function mesh:绘制一个MxN大小的二维函数H的线框图
mesh(H)
mesh(H(1:k:end, 1:k:end))
function colormap:上色用的,默认为黑色[0 0 0]
colormap([0 0 0])
function view:观察点,可以控制观察点视角位置
view(az, el) % az表示方位角,el表示仰角
Example Wireframe plotting:
>> H = fftshift(lpfilter('gaussian', 500, 500, 50)); % 绘制高斯低通滤波器
>> mesh(double(H(1:10:500, 1:10:500))) % 为了防止数据稠密,选取50个点做线框图
>> axis tight
>> colormap([0 0 0])
>> axis off
>> grid off
>> view(-25, 30)
>> view(-25, 0)
![image-20211113110539283](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211113110539283.png)
function surf:绘制表面图
surf(H)
function colormap:颜色转换
colormap(gray) % 转换为灰度
Example Surface Plotting:
>> H = fftshift(lpfilter('gaussian', 500, 500, 50)); % 绘制高斯低通滤波器
>> surf(double(H(1:10:500, 1:10:500)))
>> subplot(1,2,1)
>> surf(double(H(1:10:500, 1:10:500)))
>> axis tight
>> subplot(1,2,2)
>> surf(double(H(1:10:500, 1:10:500)))
>> axis tight
>> colormap(gray)
>> axis off
>> shading interp % 平滑曲线
![image-20211113111143227](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211113111143227.png)
如果是绘制含有两个变量的解析函数,我们用function meshgrid产生坐标值,并从这些坐标值中产生我们用于传入mesh和surf的离散矩阵。
6 Highpass (Sharpening) Frequency Domain Filters高通(锐化)频域滤波器
高通(锐化)是低通(平滑)的反处理,即通过消弱Fourier Transform的低频并保持高频相对不变来sharpen image。
相对于Lowpass滤波器的transfer function,highpass滤波器的transfer function表达式为:
H
H
P
(
u
,
v
)
=
1
−
H
L
P
(
u
,
v
)
H_{HP}(u,v)=1-H_{LP}(u,v)
HHP(u,v)=1−HLP(u,v)
则,对于上述提到的三个Lowpass filter transfer function,我们可以得到Highpass filter transfer function:
- Ideal Highpass filter:
H ( u , v ) = { 0 , if D ( u , v ) ≤ D 0 1 , if D ( u , v ) > D 0 H(u,v)=\begin{cases} 0, \text{if} \space D(u,v) \le D_0 \\ 1, \text{if} \space D(u,v) \gt D_0 \end{cases} H(u,v)={0,if D(u,v)≤D01,if D(u,v)>D0
- Butterworth Filter:
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
- Gaussian Filter:
H ( u , v ) = 1 − e x p ( − D 2 ( u , v ) 2 D 0 2 ) H(u,v)=1-exp(\frac{-D^2(u,v)}{2D_0^2}) H(u,v)=1−exp(2D02−D2(u,v))
6.1 A Function for Highpass Filtering
根据上述获得的三个高通滤波转化函数,我们可以自定义一个函数实现上述要求,这里不需要太麻烦,可以只可调用function lpfilter,最后结果用1减去即可。
function H = hpfilter(type, M, N, D0, n)
%HPFILTER Computes frequency domain highpass filters
% H = HPFILTER(TYPE, M, N, D0, n) creates the transfer function of
% a highpass filter, H, of the specified TYPE and size (M-by-N).
% Valid values for TYPE, D0, and n are:
%
% Valid value for TYPE, D0, and n are:
% 'ideal' Ideal lowpass filter with cutoff frequency D0. n need not be
% supplied. D0 must be positive.
% 'btw' Butterworth lowpass filter of order n, and cutoff D0. The
% default value for n is 1.0. D0 must be positive.
% 'gaussian'Gaussian lowpass filter with cutoff (standard deviation) D0.
% n need not be supplied. D0 must be positive.
%
% H is of floating point class single. It is returned uncentered
% for consistency with filtering function dftfilt. To view H as an
% image or mesh plot, it should be centered using Hc = fftshift (H).
%
% The transfer function Hhp of a highpass fílter is 1 - Hlp,
% where Hlp is the transfer function of the corresponding lowpass
% filter. Thus ,we can use function lpfilter to generate highpass
% filters.
if nargin == 4
n = 1; % Default value of n.
end
% Generate highpass filter.
Hlp = lpfilter(type, M, N, D0, n);
H = 1 - Hlp;
Example Highpass filters:
Ideal Filter绘制:
>> subplot(2,3,1);
>> H = fftshift(hpfilter('ideal', 500, 500, 50));
>> mesh(double(H(1:10:500,1:10:500)));
>> axis tight
>> colormap([0 0 0])
>> axis off
>> title('ideal')
>> subplot(2,3,4)
>> imshow(H, []);
>> title('ideal image')
Butterworth Filter绘制:
>> subplot(2,3,2)
>> btwH = fftshift(hpfilter('btw', 500, 500, 50));
>> mesh(double(btwH(1:10:500,1:10:500)));
>> axis tight
>> axis off
>> title('butterworth')
>> subplot(2,3,5)
>> imshow(btwH, []);
>> title('butterworth image')
Gaussian Filter绘制:
>> subplot(2,3,3)
>> gH = fftshift(hpfilter('gaussian', 500, 500, 50));
>> mesh(double(gH(1:10:500,1:10:500)));
>> axis tight
>> axis off
>> title('gaussin')
>> subplot(2,3,6)
>> imshow(gH, []);
>> title('gaussin image')
![image-20211113144924519](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211113144924519.png)
Example Highpass filtering:利用Gaussian高通滤波函数再频域内处理图片
>> PQ = paddedsize(size(f));
>> D0 = 0.05*PQ(1);
>> H = hpfilter('gaussian', PQ(1), PQ(2), D0);
>> g = dftfilt(f, H);
>> subplot(1,2,1);imshow(f);
>> subplot(1,2,2);imshow(g);
![image-20211113145455805](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211113145455805.png)
6.2 High-Frequency Emphasis Filtering高频强度滤波
Highpass filters偏离了dc的0项,因此减少了image中平均值为0的值。一种补偿方法是给HPF加上偏移量。当偏移量与滤波器乘以某个大于1的常数结合起来时,可以高亮出图像的高频区域,所以称这种方法为高频强度滤波。
同时乘数的大小也增强了低频的部分,但只要偏移量与乘数相比较小,低频增强的影响就小于高频增强的影响。
传递函数为:
H
H
F
E
(
u
,
v
)
=
a
+
b
H
H
P
(
u
,
v
)
H_{HFE}(u,v)=a+bH_{HP}(u,v)
HHFE(u,v)=a+bHHP(u,v)
其中,a是偏移量offset,b是大于1的常数。
Example Combining High-frequency emphasis and histogram equalization:
>> PQ = paddedsize(size(f));
>> D0 = 0.05*PQ(1);
>> HBW = hpfilter('btw', PQ(1), PQ(2), D0, 2);
>> H = 0.5 + 2*HBW;
>> gbw = dftfilt(f, HBW, 'fltpoint');
>> gbw = gscale(gbw);
>> ghf = dftfilt(f, H, 'fltpoint');
>> ghf = gscale(ghf);
>> ghe = histeq(ghf, 256);
>> subplot(2,2,1);imshow(f);
>> subplot(2,2,2);imshow(gbw);
>> subplot(2,2,3);imshow(ghf);
>> subplot(2,2,4);imshow(ghe);
![image-20211113153345146](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211113153345146.png)
7 Selective Filtering
分为两类:第一类是带阻(bandreject)和带通(bandpass)滤波器,第二类是陷波带阻(notcherject)和陷波带通(notchpass)滤波器。
7.1 Bandreject and Bandpass Filters
这类滤波器可以很容易用lowpass和highpass滤波器构建。表达式可写为
H
B
P
(
u
,
v
)
=
1
−
H
B
R
(
u
,
v
)
H_{BP}(u,v)=1-H_{BR}(u,v)
HBP(u,v)=1−HBR(u,v)
则,对于上述提到的三个Lowpass filter transfer function,我们可以得到带通带阻滤波器:
- Ideal filter:
H ( u , v ) = { 0 , D 0 − W 2 ≤ D ( u , v ) ≤ D 0 + W 2 1 , otherwise H(u,v)=\begin{cases} 0, D_0-\frac{W}{2} \le D(u,v)\le D_0+\frac{W}{2} \\ 1, \text{otherwise} \end{cases} H(u,v)={0,D0−2W≤D(u,v)≤D0+2W1,otherwise
函数实现为
function H = idealReject(D, D0, W)
RI = D <= D0 - (W/2); % Points of region inside the inner
% boundary of the reject band are labeled 1.
% All other points are labeled 0.
RO = D >= D0 - (W/2); % Points of region outside the outer
% boundary of the reject band are labeled 1.
% All other points are labeled 0.
H = tofloat(RO | RI); % Ideal bandreject filter.
- Butterworth Filter:
H ( u , v ) = 1 1 + [ W D ( u , v ) D 2 ( u , v ) − D 0 2 ] 2 n H(u,v)=\frac{1}{1+[\frac{WD(u,v)}{D^2(u,v)-D^2_0}]^{2n}} H(u,v)=1+[D2(u,v)−D02WD(u,v)]2n1
函数实现为
function H = btwReject(D, D0, W, n)
H = 1./(1 + (((D*W)./(D.^2 - D0.^2)).^2*n));
- Gaussian Filter:
H ( u , v ) = 1 − e x p ( − [ D 2 ( u , v ) − D 0 2 W D ( u , v ) ] 2 ) H(u,v)=1-exp(-[\frac{D^2(u,v)-D^2_0}{WD(u,v)}]^2) H(u,v)=1−exp(−[WD(u,v)D2(u,v)−D02]2)
函数实现为
function H = gaussReject(D, D0, W)
H = 1 - exp(-((D.^2 - D0.^2)./(D.*W+eps)).^2);
其中,W是对Ideal filter是带宽,对gaussian filter是粗略的截止频率,W和n结合是对butterworth的通带的宽度。
function H = bandfilter(type, band, M, N, D0, W, n)
%BANDFILTER Computes frequency domain band filters.
%
% Parameters used in the filter definitions :
% M: Number of rows in the filter .
% N: Number of col umns in the filter .
% DO : Radius of the center of the band.
% W: "Width" of the band. W is the true width only for
% ideal filters. For the other two filters this pa r ameter
% acts more like a smooth cutoff .
% n: Order of the Butterworth filter if one is specified. W
% and n interplay to determine the effective broadness of
% the reject or pass band. Higher values of both these
% parameters result in broader bands .
% Valid values of BAND are :
%
% 'reject' Bandreject filter .
% 'pass' Bandpass filter .
%
% One of these two values must. be specified for BAND.
%
% H = BANDFILTER(' ideal', BAND, M, N, DO, W) computes an M-by-N :
% ideal bandpass or bandreject filter, depending on the value of
% BAND.
%
% H = BANDFILTER( 'btw', BAND, M, N, D0, W, n) computes an M-by-N
% Butterworth filter of order n. The filter is either bandpass or
% bandreject, depending on the value of BAND. The default value of
% n is 1.
%
% H = BANDFILTER( 'gaussian', BAND, M, N, D0, W) computes an M-by-N
% gaussian filter. The filter is either bandpass or bandreject,
% depending on BAND.
%
% H is of floating point class single. It is returned uncentered
% for consistency with filtering function dftfilt. To view H as an
% image or mesh plot, it should be centered using Hc = fftshift (H) .
% Use function dftuv to set up the meshgrid arrays needed for
% computing the required distances.
[U, V] = dftuv(M, N);
% Compute the distances D(U, V).
D = hypot(U, V);
% Determine if need to use default n.
if nargin < 7
n = 1; % Default BTW filter order
end
% Begin filter computations. All filters are computed as bandreject
% filters. At the end, they are converted to bandpass if so
% specified. Use lower (type) to protect against the input being
% capitalized.
switch lower(type)
case 'ideal'
H = idealReject(D, D0, W);
case 'btw'
H = btwReject(D, D0, W, n);
case 'gaussin'
H = gaussReject(D, D0, W);
otherwise
error('Unknown filter type');
end
% Generate a bandpass filter if one was specified
if strcmp(band, 'pass')
H = 1 - H;
end
7.2 Notchreject and Notchpass Filters
陷波滤波器是更有用的选择滤波器。NF事先定义关于频率矩形的中心的邻域内的频率。零相移滤波器(Zero-phase-shift)必须关于远点对称。
以频率(u0,v0)为中心开的notch必须有相应的位于(-u0,-v0)开的notch。
Notchreject Filter可以用中心被平移到陷波滤波器中心的HP的乘积来构造。包括一对开槽Q的general form为
H
N
R
(
u
,
v
)
=
∏
k
=
1
Q
H
k
(
u
,
v
)
H
−
k
(
u
,
v
)
H_{NR}(u,v)=\prod_{k=1}^QH_k(u,v)H_{-k}(u,v)
HNR(u,v)=k=1∏QHk(u,v)H−k(u,v)
其中,
H
k
(
u
,
v
)
{H_k(u,v)}
Hk(u,v)和
H
−
k
(
u
,
v
)
{H_{-k}(u,v)}
H−k(u,v)是分别以
(
u
k
,
v
k
)
{(u_k,v_k)}
(uk,vk)和
(
−
u
k
,
−
v
k
)
{(-u_k,-v_k)}
(−uk,−vk)为中心的HPF,这些中心是以频率矩形的中心(M/2, N/2)来确定的,所以,每个Filter的distance function有:
D
k
(
u
,
v
)
=
[
(
u
−
M
/
2
−
u
k
)
2
+
(
v
−
N
/
2
−
v
k
)
2
]
1
2
D
−
k
(
u
,
v
)
=
[
(
u
−
M
/
2
+
u
k
)
2
+
(
v
−
N
/
2
+
v
k
)
2
]
1
2
D_k(u,v)=[(u-M/2-u_k)^2+(v-N/2-v_k)^2]^{\frac{1}{2}} \\ D_{-k}(u,v)=[(u-M/2+u_k)^2+(v-N/2+v_k)^2]^{\frac{1}{2}}
Dk(u,v)=[(u−M/2−uk)2+(v−N/2−vk)2]21D−k(u,v)=[(u−M/2+uk)2+(v−N/2+vk)2]21
类似于带阻和带通滤波器,我们可以将陷波带阻和带通滤波器的表达式写为
H
N
P
(
u
,
v
)
=
1
−
H
N
R
(
u
,
v
)
H_{NP}(u,v) = 1 - H_{NR}(u,v)
HNP(u,v)=1−HNR(u,v)
Example 用陷波滤波器减少波纹模式:
>> f = imread('Fig0321(a).tif');
>> subplot(2,3,1);imshow(f);
>> [M N] = size(f);
>> [f, revertClass] = tofloat(f);
>> F = fft2(f);
>> S = gscale(log(1 + abs(fftshift(F)))); % Spectrum
>> title('Image with moiré patterns');
>> subplot(2,3,2);imshow(S);title('Spectrum');
>> % Use function imtool to obtain the coordinates of the
>> % spikes interactively.
>> C1 = [99 154; 128 163];
>> % Notch filter:
>> H1 = cnotch('gaussian', 'reject', M, N, C1, 5);
>> H1 = cnotch('gaussian', 'reject', M, N, C1, 5);
>> % Compute spectrum of the filtered transform and shot it as
>> % an image.
>> P1 = gscale(fftshift(H1).*(tofloat(S)));
>> subplot(2,3,3);imshow(P1);title('Spectrum low-frequency bursts with moiré patterns');
>> % Filter Image
>> g1 = dftfilt(f, H1);
>> g1 = revertClass(g1);
>> subplot(2,3,4);imshow(g1);title('Notch Filter Result');
>> % Repeat with the following C2 to reduce the higher
>> % frequency interference components.
>> C2 = [99 154; 128 163; 49 160; 133 233; 55 132; 108 225; 112 74];
>> H2 = cnotch('gaussian', 'reject', M, N, C2, 5);
>> % Compute the spectrum of the filtered transform and show
>> % it as an image.
>> P2 = gscale(fftshift(H2).*(tofloat(S)));
>> subplot(2,3,5);imshow(P2);title('Spectrum eliminae higher frequency noise with more filters');
>> % Filter image
>> g2 = dftfilt(f, H2);
>> g2 = revertClass(g2);
>> subplot(2,3,6);imshow(g2);title('Filtered result with less noise');
![image-20211113180031202](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211113180031202.png)
陷波滤波的特殊情况包括沿着DFT的轴的取值范围进行滤波,我们可以设计一个recnotch函数使用放置在轴上的矩形完成计算。
Example Using Notch filtering to reduce periodic interference caused by malfunctioning imaging equipment:
本例中的original image附加了水平方向上的干扰信号。
>> f = imread('Fig0322(a).tif');
>> subplot(2,2,1);imshow(f);title('image');
>> [M N] = size(f);
>> [f, revertClass] = tofloat(f);
>> F = fft2(f);
>> S = gscale(log(1 + abs(fftshift(F))));
>> subplot(2,2,2);imshow(S);title('Spectrum');
>> H = recnotch('reject', 'vertical', M, N, 3, 15, 15);
>> subplot(2,2,3);imshow(fftshift(H));title('Result: NRF multiples DFT');
>> g = dftfilt(f, H);
>> g = revertClass(g);
>> subplot(2,2,4);imshow(g);title('Result: Computing IDFT');
![image-20211113183635036](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211113183635036.png)
从实验中我们可以看出,干扰是一个近似关于垂直方向周期的,所以我们期望在谱的垂直轴上找到存在的能量脉冲。我们采用窄的、矩形的陷波滤波器进行处理,结果如图三所示。最后,我们通过IDFT的计算获得改进后的结果。
为了得到空间干扰自身,我们可以使用NPF代替NRF,在vertical上隔离出干扰的频率。然后,利用IDFT得到空间干扰模式本身:
>> Hrecpass = recnotch('pass', 'vertical', M, N, 3, 15, 15);
>> interference = dftfilt(f, Hrecpass);
>> subplot(1,2,1);imshow(fftshift(Hrecpass));
>> interference = gscale(interference);
>> subplot(1,2,2);imshow(interference);
![image-20211113184509716](/Users/wanghe/Library/Application%20Support/typora-user-images/image-20211113184509716.png)