7.2 图像复原之空间滤波
图像复原（只存在噪声的复原）之空间滤波
文章目录
- 图像复原（只存在噪声的复原）之空间滤波
- 前言
- 1. 均值滤波器
- 1.1 算术平均滤波器
- 1.2 几何均值滤波器
- 1.3 谐波平均滤波器
- 1.4 反谐波平均滤波器
- 总结
前言
当一幅图像仅被加性噪声退化时，其退化模型可简化为：
g
(
x
,
y
)
=
f
(
x
,
y
)
+
η
(
x
,
y
)
G
(
u
,
v
)
=
F
(
u
,
v
)
+
N
(
u
,
v
)
\begin{aligned} &g(x,y)=f(x,y)+\eta(x,y) \\ &G(u,v)=F(u,v)+N(u,v) \end{aligned}
g(x,y)=f(x,y)+η(x,y)G(u,v)=F(u,v)+N(u,v)其中：
g
(
x
,
y
)
[
或
G
(
u
,
v
)
]
表示退化后的图像
f
(
x
,
y
)
[
或
F
(
u
,
v
)
]
表示原图像
η
(
x
,
y
)
[
或
N
(
u
,
v
)
]
表示图像的噪声
\begin{aligned} &g(x,y)[或G(u,v)]\quad表示退化后的图像\\ &f(x,y)[或F(u,v)]\quad表示原图像 \\ &\eta(x,y)[或N(u,v)]\quad表示图像的噪声\\ \end{aligned}
g(x,y)[或G(u,v)]表示退化后的图像f(x,y)[或F(u,v)]表示原图像η(x,y)[或N(u,v)]表示图像的噪声由于噪声项通常是未知的，因此我们 不能 想当然的以为得到原图像的方式为：
f
(
x
,
y
)
=
g
(
x
,
y
)
−
η
(
x
,
y
)
F
(
u
,
v
)
=
G
(
u
,
v
)
−
N
(
u
,
v
)
\begin{aligned} &f(x,y)=g(x,y)-\eta(x,y) \\ &F(u,v)=G(u,v)-N(u,v) \end{aligned}
f(x,y)=g(x,y)−η(x,y)F(u,v)=G(u,v)−N(u,v)对于只存在加性噪声的图像，可以使用空间滤波的方法来估计
f
(
x
,
y
)
f(x,y)
f(x,y)
1. 均值滤波器
1.1 算术平均滤波器
算术平均滤波器是最简单的均值滤波器，与空间域滤波中的盒式滤波器完全相同。令：
S
x
y
表示中心为
(
x
,
y
)
、大小为
m
×
n
的邻域
S_{xy}表示中心为(x,y)、大小为m\times n的邻域
Sxy表示中心为(x,y)、大小为m×n的邻域则复原的图像 ：
f
^
(
x
,
y
)
=
1
m
n
∑
(
r
,
c
)
∈
S
x
y
g
(
r
,
c
)
\hat{f}(x,y)=\frac {1} {mn}\sum_{(r,c)\in S_{xy}}g(r,c)
f^(x,y)=mn1(r,c)∈Sxy∑g(r,c)用白话来说就是使用
(
x
,
y
)
(x,y)
(x,y) 的
m
×
n
m\times n
m×n 邻域的均值作为复原图像的像素值**均值滤波会降低图像的噪声，但是也会模糊图像**。
算术平均滤波器复原图像的代码实现如下：
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
def imageRestoration(src, kernel):
"""
图像复原-算术平均滤波器
:param src: 需要复原的原图像
:param kernel: 卷积核
:return: 返回复原图像
"""
img = cv.boxFilter(src, -1, kernel) # cv2.boxFilter 方法 (默认normalize=True)
return img
if __name__ == '__main__':
img = cv.imread('Image/Fig0702.tif', 0)
img1 = imageRestoration(img, (3, 3))
img2 = imageRestoration(img, (9, 9))
img3 = imageRestoration(img, (15, 15))
# 图像显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['font.size'] = 20
plt.figure(figsize=(10, 10))
titleList = ["Original Image", "3*3", "9*9", "15*15"]
imageList = [img, img1, img2, img3]
for i in range(4):
plt.subplot(2, 2, i + 1), plt.title(titleList[i]), plt.axis('off')
plt.imshow(imageList[i], vmin=0, vmax=255, cmap='gray')
plt.tight_layout()
plt.savefig("Image/tmp.png")
plt.show()
结果分析：
如上图所示：原图像是一张被加性高斯噪声污染的图像，其分别使用3*3、9*9、15*15的算术平均滤波器进行复原处理，可以看出虽然算术平均滤波器可以进行图像复原，但同时也会造成图像模糊。
1.2 几何均值滤波器
几何均值滤波器图像复原的公式如下，令：
S
x
y
表示中心为
(
x
,
y
)
、大小为
m
×
n
的邻域
S_{xy}表示中心为(x,y)、大小为m\times n的邻域
Sxy表示中心为(x,y)、大小为m×n的邻域则复原的图像 ：
f
^
(
x
,
y
)
=
[
∏
(
r
,
c
)
∈
S
x
y
g
(
r
,
c
)
]
1
m
n
\hat{f}(x,y)={\left[\prod_{(r,c)\in S_{xy}}g(r,c)\right]}^{\frac{1}{mn}}
f^(x,y)=
(r,c)∈Sxy∏g(r,c)
mn1简单来说几何均值就是先求取中心点
(
x
,
y
)
(x,y)
(x,y) 的
m
×
n
m\times n
m×n 领域的像素之积，然后对像素之积求
1
m
n
\frac{1}{mn}
mn1 次幂。
算术平均滤波器复原图像的代码实现如下：
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
def ArithmenticMeanFilter(src, kernel):
"""
图像复原-算数平均滤波器
:param src: 需要复原的原图像
:param kernel: 卷积核
:return: 返回复原图像
"""
img = cv.boxFilter(src, -1, kernel) # cv2.boxFilter 方法 (默认normalize=True)
return img
def GeometricMeanFilter(src, kernel):
"""
图像复原-几何平均滤波器
:param src: 原图像
:param kernel: 卷积核的大小
:return:
"""
# 将图像转换为浮点数，以避免计算过程中出现溢出
src = src.astype(np.float32) + 1e-10 # 防止取对数时出现零
rows, cols = src.shape # 图像的大小
# 进行几何平均值计算时，边缘像素会无法处理，因此需要对边缘像素进行填充
pad_size = kernel // 2 # 边缘填充的大小
padded_image = np.pad(src, pad_size, mode='reflect') # 对边缘进行填充
# 遍历图像每个像素，并计算几何均值
output = np.zeros_like(src)
for i in range(rows):
for j in range(cols):
# 获取邻域
region = padded_image[i:i + kernel, j:j + kernel]
# 计算邻域内像素的几何均值
mean = np.exp(np.mean(np.log(region)))
output[i, j] = mean
# 计算出来的结果是浮点数，需要将浮点数转换为8位整数
output = np.clip(output, 0, 255).astype(np.uint8)
return output
if __name__ == '__main__':
img = cv.imread('Image/Fig0702.tif', 0)
img1 = ArithmenticMeanFilter(img, (3, 3))
img2 = GeometricMeanFilter(img, 3)
img3 = ArithmenticMeanFilter(img, (9, 9))
img4 = GeometricMeanFilter(img, 9)
# 图像显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['font.size'] = 20
plt.figure(figsize=(30, 20))
titleList = ["Original Image", "Arithmentic Mean Filter 3*3", "Geometric Mean Filter 3*3", "Original Image", "Arithmentic Mean Filter 9*9", "Geometric Mean Filter 9*9"]
imageList = [img, img1, img2, img, img3, img4]
for i in range(6):
plt.subplot(2, 3, i + 1), plt.title(titleList[i]), plt.axis('off')
plt.imshow(imageList[i], vmin=0, vmax=255, cmap='gray')
plt.tight_layout()
plt.savefig("Image/tmp.png")
plt.show()
结果分析：
如上图所示：原图像依旧是那张被加性高斯噪声污染的图像，其分别使用算术平均滤波器和几何平均滤波器进行复原处理，可以看出这两种滤波器都可以进行降噪处理，但是几何平均滤波器对图像的模糊效果要少一些，特别是当卷积核大的时候更加明显。
1.3 谐波平均滤波器
谐波平均滤波器图像复原的公式如下，令：
S
x
y
表示中心为
(
x
,
y
)
、大小为
m
×
n
的邻域
S_{xy}表示中心为(x,y)、大小为m\times n的邻域
Sxy表示中心为(x,y)、大小为m×n的邻域则复原的图像 ：
f
^
(
x
,
y
)
=
m
n
∑
(
r
,
c
)
∈
S
x
y
1
g
(
r
,
c
)
\hat{f}(x,y)=\frac{mn}{\sum_{(r,c)\in S_{xy}}\frac{1}{g(r,c)}}
f^(x,y)=∑(r,c)∈Sxyg(r,c)1mn简单解释一下。
- 第一步：用 1 1 1 除以领域内每个像素的值，即 1 g ( r , c ) \frac{1}{g(r,c)} g(r,c)1 。
- 第二步：对第一步的结果求和，即 ∑ ( r , c ) ∈ S x y 1 g ( r , c ) \sum_{(r,c)\in S_{xy}} \frac{1}{g(r,c)} ∑(r,c)∈Sxyg(r,c)1。
- 第三步：用 m n mn mn 除以第二步的结果。
谐波平均滤波器复原图像的代码实现如下：
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
def ArithmenticMeanFilter(src, kernel):
"""
图像复原-算数平均滤波器
:param src: 需要复原的原图像
:param kernel: 卷积核
:return: 返回复原图像
"""
img = cv.boxFilter(src, -1, kernel) # cv2.boxFilter 方法 (默认normalize=True)
return img
def GeometricMeanFilter(src, kernel):
"""
图像复原-几何平均滤波器
:param src: 原图像
:param kernel: 卷积核的大小
:return:
"""
# 将图像转换为浮点数，以避免计算过程中出现溢出
src = src.astype(np.float32) + 1e-10 # 防止取对数时出现零
rows, cols = src.shape # 图像的大小
# 进行几何平均值计算时，边缘像素会无法处理，因此需要对边缘像素进行填充
pad_size = kernel // 2 # 边缘填充的大小
padded_image = np.pad(src, pad_size, mode='reflect') # 对边缘进行填充
# 遍历图像每个像素，并计算几何均值
output = np.zeros_like(src)
for i in range(rows):
for j in range(cols):
# 获取邻域
region = padded_image[i:i + kernel, j:j + kernel]
# 计算邻域内像素的几何均值
mean = np.exp(np.mean(np.log(region)))
output[i, j] = mean
# 计算出来的结果是浮点数，需要将浮点数转换为8位整数
output = np.clip(output, 0, 255).astype(np.uint8)
return output
def HarmonicMeanFilter(src, kernel):
"""
图像复原-谐波平均滤波器
:param src: 原图像
:param kernel: 卷积核的大小
:return:
"""
# 将图像转换为浮点数，以避免计算过程中出现溢出
src = src.astype(np.float32) + 1e-10 # 防止取对数时出现零
rows, cols = src.shape # 图像的大小
# 进行几何平均值计算时，边缘像素会无法处理，因此需要对边缘像素进行填充
pad_size = kernel // 2 # 边缘填充的大小
padded_image = np.pad(src, pad_size, mode='reflect') # 对边缘进行填充
# 遍历图像每个像素，并计算几何均值
output = np.zeros_like(src)
for i in range(rows):
for j in range(cols):
# 获取邻域
region = padded_image[i:i + kernel, j:j + kernel]
# 计算邻域内像素的几何均值
mean = (kernel * kernel) / np.sum(1.0 / (region + 1e-10)) # 加上小值避免除零
output[i, j] = mean
# 计算出来的结果是浮点数，需要将浮点数转换为8位整数
output = np.clip(output, 0, 255).astype(np.uint8)
return output
if __name__ == '__main__':
img = cv.imread('Image/Fig0702.tif', 0)
img1 = ArithmenticMeanFilter(img, (3, 3))
img2 = GeometricMeanFilter(img, 3)
img3 = HarmonicMeanFilter(img, 3)
# 图像显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['font.size'] = 20
plt.figure(figsize=(20, 20))
titleList = ["Original Image", "Arithmentic Mean Filter 3*3", "Geometric Mean Filter 3*3",
"Harmonic Mean Filter 3*3"]
imageList = [img, img1, img2, img3]
for i in range(4):
plt.subplot(2, 2, i + 1), plt.title(titleList[i]), plt.axis('off')
plt.imshow(imageList[i], vmin=0, vmax=255, cmap='gray')
plt.tight_layout()
plt.savefig("Image/tmp.png")
plt.show()
结果分析：
如上图所示：原图像依旧是那张被加性高斯噪声污染的图像，其分别使用算术平均滤波器、几何平均滤波器和谐波平均滤波器进行复原处理，可以看出这三种滤波器都可以对高斯噪声进行降噪处理。
1.4 反谐波平均滤波器
谐波平均滤波器图像复原的公式如下，令：
S
x
y
表示中心为
(
x
,
y
)
、大小为
m
×
n
的邻域
S_{xy}表示中心为(x,y)、大小为m\times n的邻域
Sxy表示中心为(x,y)、大小为m×n的邻域则复原的图像 ：
f
^
(
x
,
y
)
=
∑
(
r
,
c
)
∈
S
x
y
g
(
r
,
c
)
Q
+
1
∑
(
r
,
c
)
∈
S
x
y
g
(
r
,
c
)
Q
\hat{f}(x,y)=\frac{\sum_{(r,c)\in S_{xy}}g(r,c)^{Q+1}}{\sum_{(r,c)\in S_{xy}}g(r,c)^{Q}}
f^(x,y)=∑(r,c)∈Sxyg(r,c)Q∑(r,c)∈Sxyg(r,c)Q+1简单解释一下，式中：
- Q Q Q 称为滤波器的阶数，这种滤波器适用于降低或消除噪声 。
- Q Q Q 为正值时，该滤波器消除胡椒噪声。
- Q Q Q 为负值时，该滤波器消除盐粒噪声。
- 然而，该滤波器不能同时消除这两种噪声。
- Q = 0 Q=0 Q=0 时，反谐波平均滤波器简化为算数平均滤波器。
- Q = − 1 Q=-1 Q=−1时，反谐波平均滤波器简化为谐波平均滤波器。
反谐波平均滤波器复原图像的代码实现如下：
import cv2 as cv
import numpy as np
import matplotlib.pyplot as plt
def Inv_HarmonicMeanFilter(src, kernel, Q):
"""
图像复原-反谐波平均滤波器
:param src: 原图像
:param kernel: 卷积核的大小
:param Q: 滤波器的阶数
:return:
"""
# 将图像转换为浮点数，以避免计算过程中出现溢出
src = src.astype(np.float32) + 1e-10 # 防止取对数时出现零
rows, cols = src.shape # 图像的大小
# 进行几何平均值计算时，边缘像素会无法处理，因此需要对边缘像素进行填充
pad_size = kernel // 2 # 边缘填充的大小
padded_image = np.pad(src, pad_size, mode='reflect') # 对边缘进行填充
# 遍历图像每个像素，并计算几何均值
output = np.zeros_like(src)
for i in range(rows):
for j in range(cols):
# 获取邻域
region = padded_image[i:i + kernel, j:j + kernel]
# 计算分子和分母
numerator = np.sum(np.power(region, Q + 1))
denominator = np.sum(np.power(region, Q) + 1e-10) # 加小值避免除零
# 计算反谐波均值
output[i, j] = numerator / denominator
# 计算出来的结果是浮点数，需要将浮点数转换为8位整数
output = np.clip(output, 0, 255).astype(np.uint8)
return output
if __name__ == '__main__':
pepperImg = cv.imread('Image/Fig0703.tif', 0) # 胡椒噪声图像
saltImage = cv.imread('Image/Fig0704.tif', 0) # 盐粒噪声图像
img1 = Inv_HarmonicMeanFilter(pepperImg, 3, 1.5)
img2 = Inv_HarmonicMeanFilter(saltImage, 3, -1.5)
# 当阶数取反，则会出现灾难性后果
img3 = Inv_HarmonicMeanFilter(pepperImg, 3, -1.5)
img4 = Inv_HarmonicMeanFilter(saltImage, 3, 1.5)
# 图像显示
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['font.size'] = 20
plt.figure(figsize=(20, 30))
titleList = ["Original Pepper Image", "Original Salt Image", "Inv_Harmonic Mean Filter Q = 1.5",
"Inv_Harmonic Mean Filter Q = -1.5", "Q = -1.5", "Q = 1.5"]
imageList = [pepperImg, saltImage, img1, img2, img3, img4]
for i in range(6):
plt.subplot(3, 2, i + 1), plt.title(titleList[i]), plt.axis('off')
plt.imshow(imageList[i], vmin=0, vmax=255, cmap='gray')
plt.tight_layout()
plt.savefig("Image/tmp.png")
plt.show()
结果分析：
上图中第一和第二张图像分别是胡椒噪声图像和盐粒噪声图像，接下来使用了反谐波平均滤波器对图像进行了复原处理，可以看出效果相当不错，但是反谐波平均滤波器有一个缺点，就是必须要提前知道是胡椒噪声还是盐粒噪声，否则会出现上面第五、第六张图像的问题。
总结
| 滤波器类型 | 定义 | 适用噪声类型 | 对噪声的处理效果 | 保留图像细节的能力 |
|---|---|---|---|---|
| 算术平均滤波器 | 对邻域内像素值求算术平均 | 高斯噪声 | 对高斯噪声效果较好 | 易模糊边缘细节 |
| 几何平均滤波器 | 对邻域内像素值求几何平均 | 乘性噪声、对数噪声 | 对高斯噪声效果一般 | 能保留更多细节 |
| 谐波平均滤波器 | 对邻域内像素值求倒数平均 | 盐噪声（高亮噪声） | 抑制盐噪声效果好 | 边缘保留较好 |
| 反谐波平均滤波器 | 根据参数 Q Q Q 有选择性地去噪 | 椒盐噪声（混合噪声） | 盐噪声或椒噪声（依赖 Q Q Q 值） | 灵活去噪能力强 |