数字图像处理——频率域平滑锐化图像常用滤波器

频率域滤波器

本文主要是动手实现一些常见的滤波器，并展示其频率域和空间域图像。

这里介绍的常见滤波器为：（高通/低通）理想滤波、布特沃斯滤波器和高斯滤波器，根据书本 P169 页，他们的区别与联系如下：

这三种滤波器涵盖了从非常急剧（理想）的滤波到非常平滑（高斯）的滤波范围。布特沃斯滤波器有一个参数，它成为滤波器的“阶数”。当阶数较高时，布特沃斯滤波器接近于理想滤波器。对于较低的阶数值，布特沃斯滤波器更像高斯滤波器。这样，布特沃斯滤波器可视为两种“极端”滤波器的过渡。

为了方便后续实验，对一些辅助代码模块化编写，主要以可视化部分为主，见文末附录部分。

理想滤波器

所谓"理想"是指无法通过硬件实现的硬截断

理想低通滤波器 ILPF

在圆外“阻断”所有频率，而在圆内无衰减的通过所有频率，这种二维低通滤波器称为理想低通滤波器（ILPF），由下面的函数确定

\[ H_{ILPF}(u,v) = \left \{ \begin{aligned} 1, & D(u,v) \le D_0 \\ 0, & D(u,b) > D_0 \end{aligned} \right. \]

其中\(D_0\)是一个正常数，\(D(u,v)\)表示频率域中的点\((u,v)\)距离频率域中心\((\frac{P}{2},\frac{Q}{2})\)的距离。

理想高通滤波器 IHPF

与低通类似，高通是将阈值的圆内“阻断”所有频率，而在圆外无衰减的通过所有频率，描述如下

\[ H_{IHPF}(u,v) = \left \{ \begin{aligned} 0, & D(u,v) \le D_0 \\ 1, & D(u,b) > D_0 \end{aligned} \right. \]

代码实现

代码实现低通滤波器并展示其频率域透视图、频率域图像显示、空间域图像显示和径向剖面图。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
# 获得指定大小的理想滤波器
def getIdealMask(mask_shape, filter_d0,hl_type):
    assert hl_type in ("lpf","hpf")
    rows,cols = mask_shape[0],mask_shape[1]
    crow = rows/2
    ccol = cols/2
    mask = np.zeros((rows,cols))
    for i in range(rows):
        for j in range(cols):
            dis = sqrt((i-crow)**2 + (j-ccol)**2)
            if hl_type == "lpf":
                if dis <= filter_d0:
                    mask[i,j] = 1
                else:
                    mask[i,j] = 0
            elif hl_type == "hpf":
                if dis <= filter_d0:
                    mask[i,j] = 0
                else:
                    mask[i,j] = 1
    return mask

# 测试ILPF
# 参数设置
mask_shape = (100,100)
d = 20
filter_type = "lpf"
# 获得滤波器
myfilter = getIdealMask(mask_shape,d,filter_type)
# 绘图
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "ILPF Perspective Axes3D", cmap = "gray")
# 不想进行列表解析,需要调用frompyfunc构建np可以用的分段函数
ufunc1 = np.frompyfunc(lambda x: 0 if (x-d)>0 else 1, 1, 1)
drawCurv(ax2,[ufunc1],["ILPF"],d,title = "ILPF Curv")
drawPanel(ax3,myfilter,title = "ILPF Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "ILPF Spatial Panel Axes2D")
plt.show()

# 测试IHPF
d = 20
filter_type = "hpf"
myfilter = getIdealMask(mask_shape,d,filter_type)
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "IHPF Perspective Axes3D", cmap = "gray")
# 不想进行列表解析,需要调用frompyfunc构建np可以用的分段函数
ufunc1 = np.frompyfunc(lambda x: 1 if (x-d)>0 else 0, 1, 1)
drawCurv(ax2,[ufunc1],["IHPF"],d,title = "Curv")
drawPanel(ax3,myfilter,title = "IHPF Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "IHPF Spatial Panel Axes2D")
plt.show()

理想低通 ILPF 结果四图

png
png

理想高通 IHPF 结果四图

png
png

布特沃斯滤波器

可通过硬件实现，可以通过阶数进行控制，一些资料中又称之为“巴特沃斯滤波器”。

布特沃斯低通滤波器 BLPF

截止频率位于距原点\(D_0\)的\(n\)阶布特沃斯滤波器（BLPF）的传递函数定义为：

\[H_{BLPF}(u,v) = \dfrac{1}{ 1 + {[ \dfrac {D(u,v)}{D_0} ]}^{2n} }\]

布特沃斯高通滤波器 BHPF

对应的传递函数定义为：

\[H_{BHPF}(u,v) = \dfrac{1}{ 1 + {[ \dfrac {D_0}{D_(u,v)} ]}^{2n} }\]

（分母分子颠倒）

两式中\(n\)对应了即阶参数，下面的代码给出巴特沃斯滤波器的实现，其频率域透视图、频率域图像显示、空间域图像显示和径向剖面图，曲线图绘制出不同阶下的取值。

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
# 获得指定大小的布特沃斯滤波器
def getButterworthMask(mask_shape,filter_d0,hl_type,butter_n = 1):
    assert hl_type in ("lpf","hpf")
    rows,cols = mask_shape[0],mask_shape[1]
    crow = rows/2
    ccol = cols/2
    mask = np.zeros((rows,cols))
    for i in range(rows):
        for j in range(cols):
            dis = sqrt((i-crow)**2 + (j-ccol)**2)
            if hl_type == "lpf":
                mask[i,j] = 1.0/(1+(dis/filter_d0)**(2*butter_n))
            elif hl_type == "hpf":
                # 除以0情况特判一下
                if np.abs(dis)<eps:
                    mask[i,j] = 0
                else:
                    mask[i,j] = 1.0/(1+(filter_d0/dis)**(2*butter_n))

    return mask

# 测试BLPF
# 参数设置
mask_shape = (100,100)
d = 20
filter_type = "lpf"
# 获得滤波器
myfilter = getButterworthMask(mask_shape,d,filter_type, butter_n=2)
# 绘图
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "BLPF(n=2) Perspective Axes3D", cmap = "gray")
funcs = []
labels = []
for i in range(1,6):
    labels.append("BLPF "+"n="+str(i))
funcs.append(lambda x:1.0/(1+(x/d)**(2*1)))
funcs.append(lambda x:1.0/(1+(x/d)**(2*2)))
funcs.append(lambda x:1.0/(1+(x/d)**(2*3)))
funcs.append(lambda x:1.0/(1+(x/d)**(2*4)))
funcs.append(lambda x:1.0/(1+(x/d)**(2*5)))

drawCurv(ax2,funcs,labels,d,title = "BLPF(n=1,2,3,4,5) Curv")
drawPanel(ax3,myfilter,title = "BLPF(n=2) Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "BLPF(n=2) Spatial Panel Axes2D")
plt.show()

# 测试BHPF
d = 20
filter_type = "hpf"
myfilter = getButterworthMask(mask_shape,d,filter_type, butter_n=1)
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "BHPF(n=2) Perspective Axes3D", cmap = "gray")
funcs = []
labels = []
for i in range(1,6):
#     funcs.append(lambda x:1.0/(1+(d/x)**(2*i)))
#     ufunc = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*i)), 1, 1)
#     funcs.append(ufunc)
    labels.append("BHPF "+"n="+str(i))
ufunc1 = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*1)), 1, 1)
ufunc2 = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*2)), 1, 1)
ufunc3 = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*3)), 1, 1)
ufunc4 = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*4)), 1, 1)
ufunc5 = np.frompyfunc(lambda x: 0 if np.abs(x)<eps else 1.0/(1+(d/x)**(2*5)), 1, 1)
funcs.append(ufunc1)
funcs.append(ufunc2)
funcs.append(ufunc3)
funcs.append(ufunc4)
funcs.append(ufunc5)


drawCurv(ax2,funcs,labels,d,title = "BHPF(n=1,2,3,4,5) Curv")
drawPanel(ax3,myfilter,title = "BHPF(n=2) Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "BHPF(n=2) Spatial Panel Axes2D")
plt.show()

布特沃斯低通 BLPF 结果四图

png
png

布特沃斯高通 BHPF 结果四图

png
png

高斯低通滤波器 GLPF

高斯低通滤波器二维形式由下式给处：

\[H_{GLPF}(u,v) = e^{\dfrac{-D^2(u,v)}{2 \sigma^2}}\]

\(\sigma\)描述了中心的扩散速度，和其他滤波器描述式统一，通过令\(\sigma = D_0\)，可以用表示其他滤波器的方法表示高斯滤波器。

\[H_{GLPF}(u,v) = e^{\dfrac{-D^2(u,v)}{2 D_0^2}}\]

高斯高通滤波器 GHPF

如下：

\[H_{GHPF}(u,v) =1 - e^{\dfrac{-D^2(u,v)}{2 D_0^2}}\]

代码实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
# 获得指定大小的高斯滤波器
def getGaussianMask(mask_shape,filter_d0,hl_type):
    assert hl_type in ("lpf","hpf")
    rows,cols = mask_shape[0],mask_shape[1]
    crow = rows/2
    ccol = cols/2
    mask = np.zeros((rows,cols))
    for i in range(rows):
        for j in range(cols):
            dis = sqrt((i-crow)**2 + (j-ccol)**2)
            if hl_type == "hpf":
                mask[i,j] = 1-np.exp(-(dis**2) / (2*(filter_d0**2)))
            elif hl_type == "lpf":
                mask[i,j] = np.exp(-(dis**2)/(2*(filter_d0**2)))
    return mask

# 测试GLPF
# 参数设置
mask_shape = (200,200)
d = 20
filter_type = "lpf"
# 获得滤波器
myfilter = getGaussianMask(mask_shape,d,filter_type)
# 绘图
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "GLPF Perspective Axes3D", cmap = "gray")
drawCurv(ax2,[lambda x:np.exp(-(x**2)/(2*(d**2)))],["GLPF"],d,title = "GLPF Curv")
drawPanel(ax3,myfilter,title = "GLPF Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "GLPF Spatial Panel Axes2D")
plt.show()

# 测试GHPF
d = 20
filter_type = "hpf"
myfilter = getGaussianMask(mask_shape,d,filter_type)
plt.figure(figsize=(12,12))
ax1=plt.subplot(221,projection = "3d")
ax2=plt.subplot(222)
ax3=plt.subplot(223)
ax4=plt.subplot(224)
drawPerspective(ax1,myfilter,title = "GLPF Perspective Axes3D", cmap = "gray")
drawCurv(ax2,[lambda x:1-np.exp(-(x**2)/(2*(d**2)))],["GHPF"],d,title = "GLPF Curv")
drawPanel(ax3,myfilter,title = "GLPF Frequency Panel Axes2D")
spatial_myfilter = frequencyToSpatial(myfilter)
drawPanel(ax4,spatial_myfilter,title = "GLPF Spatial Panel Axes2D")
plt.show()

高斯低通 GLPF 结果四图

png
png

高斯高通 GHPF 结果四图

png
png

总结

在规定滤波器为 100x100，阈值为 20 时，可以明显观察到，理想滤波器->高阶布特沃斯滤波器->低阶布特沃斯滤波器->高斯滤波器，可以由函数 Curv 看出对应的过渡。

我们同时也发现理想滤波器确实会存在振铃特性，这个将在后面的文章中再做分析学习。

附录

辅助代码

1
2
3
4
5
6
7
8
import numpy as np
import cv2
import imageio
import matplotlib.pyplot as plt
from math import sqrt
from mpl_toolkits.mplot3d import Axes3D

eps = np.finfo(float).eps

绘制三维透视图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def drawPerspective(handleax,input_matrix,title=None,cmap = "gray"):
    handleax.set_title(title)
    handleax.set_zlabel('Z')  # 坐标轴
    handleax.set_ylabel('Y')
    handleax.set_xlabel('X')
    x,y = input_matrix.shape
    X = np.arange(0,x,1)
    Y = np.arange(0,y,1)
    # 由于图像x,y坐标和 meshigrid出来是互反的
    # 这里需要调转一下
    # 否则会出现mismatch的现象
    X,Y = np.meshgrid(Y, X)
    handleax.plot_surface(Y, X, input_matrix, cmap=cmap)
#     handleax.plot_wireframe(Y, X, input_matrix, cmap=cmap)

绘制平面图

1
2
3
4
5
def drawPanel(handleax,input_matrix,title=None,cmap = "gray"):
    handleax.set_title(title)
    handleax.set_ylabel('Y')
    handleax.set_xlabel('X')
    handleax.imshow(input_matrix,cmap = cmap)

绘制曲线图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def drawCurv(handleax,functions,labels,filter_d0,title=None,cmap = "gray"):
    # 绘制从0到 3D_0的函数剖面图
    handleax.set_title(title)
    handleax.set_xlabel("$D(u,v)$")
    handleax.set_ylabel("$H(u,v)$")
    # 标出D_0点
    # handleax.annotate(r"$D_0$", xy = (filter_d0,0) , weight='heavy')

    for func,lab in zip(functions,labels):
        # 对每一对func和label绘图 （针对需要画多条线的情况）
        X = np.arange(0,3*filter_d0+1,0.1)
        Y = func(X)
        handleax.plot(X,Y,label = lab)
    handleax.set_xticks([0,filter_d0])
    handleax.set_xticklabels(["$0$","$D_0$"])
    handleax.legend()

频率域转空间域

1
2
3
4
5
6
7
def frequencyToSpatial(input_matrix):
    # 这里不太明白为什么shift与否最后都需要添加一个fftshift来得到想要的空间域图像
    shift_input_matrix = np.fft.ifftshift(input_matrix)
#     shift_input_matrix = input_matrix
    spatial_img = np.abs(np.fft.ifft2(shift_input_matrix))
    spatial_img = np.fft.fftshift(spatial_img)
    return spatial_img

对三维透视可视化代码进行测试

1
2
3
4
5
6
7
8
9
10
11
12
13
# 测试1
test_img1 = cv2.imread('./DIP3E_Original_Images_CH02/Fig0207(a)(gray level band).tif',0)
test_img2 = cv2.imread('./DIP3E_CH07_Original_Images/Fig0726(a).tif',0)
plt.figure(figsize=(16,12))
ax1 = plt.subplot(221,projection = "3d")
ax2 = plt.subplot(222)
ax3 = plt.subplot(223,projection = "3d")
ax4 = plt.subplot(224)
drawPerspective(ax1,test_img1,"Test Fig0726(a) Axes3D","rainbow")
drawPanel(ax2,test_img1,"Test Fig0726(a) Axes2D","rainbow")
drawPerspective(ax3,test_img2,"Test Fig0726(a) Axes3D","rainbow")
drawPanel(ax4,test_img2,"Test Fig0726(a) Axes2D","rainbow")
plt.show()

测试结果

png
png

对绘制函数曲线代码进行测试

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 测试2
test_d0 = 10
test_funcs1 = []
test_funcs1.append(lambda x: np.sin(x))
test_funcs1.append(lambda x: np.cos(x))
test_labels1 = ["$sin(x)$","$cos(x)$"]

test_funcs2 = []
test_funcs2.append(lambda x: np.sqrt(x))
test_funcs2.append(lambda x: np.exp(-(x**2)/(2*(test_d0)**2)))
test_labels2 = ["$x$","$e^{\dfrac{-D^2(u,v)}{2 D_0^2}}$"]

plt.figure(figsize=(14,7))
ax1 = plt.subplot(121)
ax2 = plt.subplot(122)
drawCurv(ax1,test_funcs1,test_labels1,test_d0,title="Test sin(x) and cos(x)")
drawCurv(ax2,test_funcs2,test_labels2,test_d0,title="Test x and gaussian")
plt.show()

测试结果

png
png