第5章 Python 数字图像处理(DIP) - 图像复原与重建9 - 空间滤波 - 均值滤波器 - 算术平均、几何平均、谐波平均、反谐波平均滤波器

最新推荐文章于 2023-09-21 10:20:59 发布

jasneik

最新推荐文章于 2023-09-21 10:20:59 发布

阅读量1.8k

点赞数 8

分类专栏： # 第5章图像复原与重建数字图像处理(DIP)--Python 文章标签：图像处理 opencv numpy python 图像识别

本文链接：https://blog.csdn.net/jasneik/article/details/115444357

版权

数字图像处理(DIP)--Python 同时被 2 个专栏收录

57 篇文章 126 订阅

订阅专栏

第5章图像复原与重建

17 篇文章 22 订阅

订阅专栏

标题

- 只存在噪声的复原 - 空间滤波
- - 均值滤波器

只存在噪声的复原 - 空间滤波

仅被加性噪声退化
$\eta(x, y) \tag{5.21}$
和
$\tag{5.22}$

噪声项通常是未知的，因此不能直接从退化图像中减去噪声项来得到复原的图像。

对于周期噪声，只是个例外而不是规律，我们可以用谱 $G (u, v)$ 来估计 $N (u, v)$ ，从 $G (u, v)$ 中减去 $N (u, v)$ 能够得到原图像的一个估计。

均值滤波器

算术平均滤波器

算术平均滤波器是最简单的均值滤波器。
$\hat{f}(x, y) = \frac{1}{mn} \sum_{(r,c)\in S_{xy}} g(r,c) \tag{5.23}$

$S_{xy}$ 表示中心为 $(x, y)$ 、大小为 $m\times{n}$ 的矩形子图像窗口（邻域）的一组坐标。算术平均滤波器在由 $S_{xy}$ 定义的区域中，计算被污染图像 $g (x, y)$ 的平均值。复原的图像 $\hat{f}$ 在 $(x, y)$ 处的值，是使用 $S_{xy}$ 定义的邻域中的像素算出的算术平均值。

均值滤波平滑图像中的局部变化，它会降低图像中的噪声，但会模糊图像

def arithmentic_mean(image, kernel):
    """
    define arithmentic mean filter, math: $$\hat{f}(x, y) = \frac{1}{mn} \sum_{(r,c)\in S_{xy}} g(r,c)$$
    param image: input image
    param kerne: input kernel, actually use kernel shape
    return: image after arithmentic mean filter, 
    """
    img_h = image.shape[0]
    img_w = image.shape[1]

    m, n = kernel.shape[:2]

    padding_h = int((m -1)/2)
    padding_w = int((n -1)/2)
    
    # 这样的填充方式，可以奇数核或者偶数核都能正确填充
    image_pad = np.pad(image, ((padding_h, m - 1 - padding_h), \
                                  (padding_w, n - 1 - padding_w)), mode="edge")

    image_mean = image.copy()
    for i in range(padding_h, img_h + padding_h):
        for j in range(padding_w, img_w + padding_w):
            temp = np.sum(image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1])
            image_mean[i - padding_h][j - padding_w] = 1/(m * n) * temp
    return image_mean

# 算术平均滤波器
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接读为灰度图像

mean_kernal = np.ones([3, 3])
mean_kernal = mean_kernal / mean_kernal.size

img_arithmentic = arithmentic_mean(img_ori, kernel=mean_kernal)

img_cv2_mean = cv2.filter2D(img_ori, ddepth= -1, kernel=mean_kernal)

plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_arithmentic, 'gray'), plt.title('Self Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_cv2_mean, 'gray'), plt.title('CV2 mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

几何均值滤波器

$\hat{f}(x, y) = \Bigg[\prod_{(r,c)\in S_{xy}} g(r,c) \Bigg]^{\frac{1}{mn}} \tag{5.24}$

几何均值滤波器实现的平滑可与算术平均滤波器的相比，但损失的图像细节更少

def geometric_mean(image, kernel):
    """
    define geometric mean filter, math: $$\hat{f}(x, y) = \Bigg[\prod_{(r,c)\in S_{xy}} g(r,c) \Bigg]^{\frac{1}{mn}}$$
    param image:  input image
    param kerne:  input kernel, actually use kernel shape
    return: image after geometric mean filter, 
    """
    img_h = image.shape[0]
    img_w = image.shape[1]

    m, n = kernel.shape[:2]
    order = 1 / (kernel.size)

    padding_h = int((m -1)/2)
    padding_w = int((n -1)/2)
    
    # 这样的填充方式，可以奇数核或者偶数核都能正确填充
    image_pad = np.pad(image.copy(), ((padding_h, m - 1 - padding_h), \
                                  (padding_w, n - 1 - padding_w)), mode="edge")

    image_mean = image.copy()
    # 这里要指定数据类型，但指定是uint64或者float64，但结果不正确，反而乘以1.0，也是float64，但却让结果正确
    for i in range(padding_h, img_h + padding_h):
        for j in range(padding_w, img_w + padding_w):
            prod = np.prod(image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1]*1.0)
            image_mean[i - padding_h][j - padding_w] = np.power(prod, order)

    return image_mean

def geometric_mean(image, kernel):
    """
    :param image: input image
    :param kernel: input kernel
    :return: image after convolution
    """
    
    img_h = image.shape[0]
    img_w = image.shape[1]

    kernel_h = kernel.shape[0]
    kernel_w = kernel.shape[1]

    # padding
    padding_h = int((kernel_h -1)/2)
    padding_w = int((kernel_w -1)/2)

    image_pad = np.pad(image.copy(), (padding_h, padding_w), mode="constant", constant_values=1)

    image_convol = image.copy()
    for i in range(padding_h, img_h + padding_h):
        for j in range(padding_w, img_w + padding_w):
            temp = np.prod(image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1] * kernel)
            image_convol[i - padding_h][j - padding_w] = np.power(temp, 1/kernel.size)

    return image_convol

# 几何均值滤波器
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接读为灰度图像

mean_kernal = np.ones([3, 3])
arithmetic_kernel = mean_kernal / mean_kernal.size

img_geometric = geometric_mean(img_ori, kernel=mean_kernal)

img_arithmentic = arithmentic_mean(img_ori, kernel=arithmetic_kernel)

plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_arithmentic, 'gray'), plt.title('Arithmetic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
# plt.show()

在这里插入图片描述

谐波平均滤波器

$\hat{f}(x, y) = \cfrac{mn}{\sum_{(r,c)\in S_{xy}} \cfrac{1}{g(r,c)}} \tag{5.25}$

可以处理盐粒噪声，又能处理类似于高斯噪声的其他噪声，但不能处理胡椒噪声

def harmonic_mean(image, kernel):
    """
    define harmonic mean filter, math: $$\hat{f}(x, y) = \Bigg[\prod_{(r,c)\in S_{xy}} g(r,c) \Bigg]^{\frac{1}{mn}}$$
    param image:  input image
    param kerne:  input kernel, actually use kernel shape
    return: image after harmonic mean filter, 
    """
    epsilon = 1e-8
    img_h = image.shape[0]
    img_w = image.shape[1]

    m, n = kernel.shape[:2]
    order = kernel.size

    padding_h = int((m -1)/2)
    padding_w = int((n -1)/2)
    
    # 这样的填充方式，可以奇数核或者偶数核都能正确填充
    image_pad = np.pad(image.copy(), ((padding_h, m - 1 - padding_h), \
                                  (padding_w, n - 1 - padding_w)), mode="edge")

    image_mean = image.copy()
    # 这里要指定数据类型，但指定是uint64或者float64，但结果不正确，反而乘以1.0，也是float64，但却让结果正确
    # 要加上epsilon，防止除0
    for i in range(padding_h, img_h + padding_h):
        for j in range(padding_w, img_w + padding_w):
            temp = np.sum(1 / (image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1]*1.0 + epsilon))
            image_mean[i - padding_h][j - padding_w] = order / temp

    return image_mean

# 谐波滤波处理高斯噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接读为灰度图像

mean_kernal = np.ones([3, 3])
arithmetic_kernel = mean_kernal / mean_kernal.size

img_geometric = geometric_mean(img_ori, kernel=mean_kernal)

img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernal)

plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

# 谐波滤波处理盐粒噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(b)(circuit-board-salt-prob-pt1).tif', 0) #直接读为灰度图像

mean_kernal = np.ones([3, 3])
arithmetic_kernel = mean_kernal / mean_kernal.size

img_geometric = geometric_mean(img_ori, kernel=mean_kernal)

img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernal)

plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

# 谐波滤波处理胡椒噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) #直接读为灰度图像

mean_kernal = np.ones([3, 3])
arithmetic_kernel = mean_kernal / mean_kernal.size

img_geometric = geometric_mean(img_ori, kernel=mean_kernal)

img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernal)

plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

反（逆）谐波平均滤波器

$\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} \tag{5.26}$

Q称为滤波器的阶数。这种滤波器适用于降低或消除椒盐噪声。Q值为正时，该滤波器消除胡椒噪声；Q值为负时，该滤波器消除盐粒噪声。然而，该滤波器不能同时消除这两种噪声。注意当 $Q = 0$ 时，简化为算术平均滤波器；当 $Q = - 1$ 时，简化为谐波平均滤波器。

def inverse_harmonic_mean(image, kernel, Q=0):
    """
    define inverse harmonic mean filter, 
    math: $$\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}$$
    param image : input image
    param kernel: input kernel, actually use kernel shape
    param Q     : input order of the filter, default is 0, which equal to arithmentic mean filter, while is -1 is harmonic mean filter
    return: image after inverse harmonic mean filter, 
    """
    epsilon = 1e-8
    img_h = image.shape[0]
    img_w = image.shape[1]

    m, n = kernel.shape[:2]

    padding_h = int((m -1)/2)
    padding_w = int((n -1)/2)
    
    # 这样的填充方式，可以奇数核或者偶数核都能正确填充
    image_pad = np.pad(image.copy(), ((padding_h, m - 1 - padding_h), \
                                  (padding_w, n - 1 - padding_w)), mode="edge")

    image_mean = image.copy()
    # 这里要指定数据类型，但指定是uint64或者float64，但结果不正确，反而乘以1.0，也是float64，但却让结果正确
    # 要加上epsilon，防止除0
    for i in range(padding_h, img_h + padding_h):
        for j in range(padding_w, img_w + padding_w):
            temp = image_pad[i-padding_h:i+padding_h+1, j-padding_w:j+padding_w+1] * 1.0 + epsilon
            # image_mean[i - padding_h][j - padding_w] = np.sum(temp**(Q+1)) / np.sum(temp**Q + epsilon)
            image_mean[i - padding_h][j - padding_w] = np.sum(np.power(temp, (Q+1))) / np.sum(np.power(temp, Q) + epsilon)

    return image_mean

# 反谐波滤波处理胡椒噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) #直接读为灰度图像

mean_kernel = np.ones([3, 3])
arithmetic_kernel = mean_kernel / mean_kernel.size

img_inverse_harmonic = inverse_harmonic_mean(img_ori, kernel=mean_kernel, Q=1.5)

img_harmonic_mean = harmonic_mean(img_ori, kernel=mean_kernel)

plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_inverse_harmonic, 'gray'), plt.title('Inverse Harmonic Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_harmonic_mean, 'gray'), plt.title('Harmonic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

# 反谐波滤波处理椒盐噪声
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0510(a)(ckt-board-saltpep-prob.pt05).tif', 0) #直接读为灰度图像

mean_kernel = np.ones([3, 3])
arithmetic_kernel = mean_kernel / mean_kernel.size

img_inverse_harmonic = inverse_harmonic_mean(img_ori, kernel=mean_kernel, Q=1.5)
img_arithmentic_mean = inverse_harmonic_mean(img_ori, kernel=mean_kernel, Q=0)
# img_arithmentic_mean = arithmentic_mean(img_ori, kernel=mean_kernel)

plt.figure(figsize=(18, 6))
plt.subplot(131), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(132), plt.imshow(img_inverse_harmonic, 'gray'), plt.title('Inverse Harmonic Mean'), plt.xticks([]), plt.yticks([])
plt.subplot(133), plt.imshow(img_arithmentic_mean, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

下面是各种滤波器的对比与总结

# 算术平均滤波器和几何均值滤波器对比
img_ori = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(a)(ckt-board-orig).tif', 0) #直接读为灰度图像
img_noise = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0507(b)(ckt-board-gauss-var-400).tif', 0) #直接读为灰度图像

mean_kernel = np.ones([3, 3])
arithmetic_kernel = mean_kernel / mean_kernel.size

img_arithmentic = arithmentic_mean(img_noise, kernel=mean_kernel)
img_geometric   = geometric_mean(img_noise, kernel=mean_kernel)

plt.figure(figsize=(12, 12))
plt.subplot(221), plt.imshow(img_ori, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(222), plt.imshow(img_noise, 'gray'), plt.title('Noisy'), plt.xticks([]), plt.yticks([])
plt.subplot(223), plt.imshow(img_arithmentic, 'gray'), plt.title('Arithmentic mean'), plt.xticks([]), plt.yticks([])
plt.subplot(224), plt.imshow(img_geometric, 'gray'), plt.title('Geomentric mean'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

反谐波滤波器可以处理胡椒与盐粒噪声。在处理胡椒噪声的时候，Q值一般为正值；在处理盐粒噪声时，Q值一般为负值。错误的选择Q值，不仅不能得到好的滤波效果，反而带来灾难性的结果。

# 反谐波滤波器分别使用不同的Q值对胡椒与盐粒噪声的滤波器
img_pepper = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) 
img_salt = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(b)(circuit-board-salt-prob-pt1).tif', 0)

mean_kernel = np.ones([3, 3])
arithmetic_kernel = mean_kernel / mean_kernel.size

inverse_harmonic_15 = inverse_harmonic_mean(img_pepper, kernel=mean_kernel, Q=1.5)
inverse_harmonic_15_ = inverse_harmonic_mean(img_salt, kernel=mean_kernel, Q=-1.5)

plt.figure(figsize=(12, 12))
plt.subplot(221), plt.imshow(img_pepper, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(222), plt.imshow(img_salt, 'gray'), plt.title('Noisy'), plt.xticks([]), plt.yticks([])
plt.subplot(223), plt.imshow(inverse_harmonic_15, 'gray'), plt.title('Pepper Noise Q = 1.5'), plt.xticks([]), plt.yticks([])
plt.subplot(224), plt.imshow(inverse_harmonic_15_, 'gray'), plt.title('Salt Noise Q = -1.5'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

# 错误的使用Q值，反谐波滤波器得到的是严重的后果
img_pepper = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(a)(circuit-board-pepper-prob-pt1).tif', 0) 
img_salt = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0508(b)(circuit-board-salt-prob-pt1).tif', 0)

mean_kernel = np.ones([3, 3])
arithmetic_kernel = mean_kernel / mean_kernel.size

inverse_harmonic_15 = inverse_harmonic_mean(img_pepper, kernel=mean_kernel, Q=-1.5)
inverse_harmonic_15_ = inverse_harmonic_mean(img_salt, kernel=mean_kernel, Q=1.5)

plt.figure(figsize=(12, 12))
plt.subplot(221), plt.imshow(img_pepper, 'gray'), plt.title('Original'), plt.xticks([]), plt.yticks([])
plt.subplot(222), plt.imshow(img_salt, 'gray'), plt.title('Noisy'), plt.xticks([]), plt.yticks([])
plt.subplot(223), plt.imshow(inverse_harmonic_15, 'gray'), plt.title('Pepper Noise Q = -1.5'), plt.xticks([]), plt.yticks([])
plt.subplot(224), plt.imshow(inverse_harmonic_15_, 'gray'), plt.title('Salt Noise Q = 1.5'), plt.xticks([]), plt.yticks([])
plt.tight_layout()
plt.show()

在这里插入图片描述

jasneik

关注

8
点赞
踩
39

收藏

觉得还不错? 一键收藏
打赏
1
评论
第5章 Python 数字图像处理(DIP) - 图像复原与重建9 - 空间滤波 - 均值滤波器 - 算术平均、几何平均、谐波平均、反谐波平均滤波器

标题只存在噪声的复原 - 空间滤波均值滤波器算术平均滤波器几何均值滤波器谐波平均滤波器反（逆）谐波平均滤波器只存在噪声的复原 - 空间滤波仅被加性噪声退化g(x,y)=f(x,y)+η(x,y)(5.21)g(x, y) = f(x, y) + \eta(x, y) \tag{5.21}g(x,y)=f(x,y)+η(x,y)(5.21)和G(u,v)=F(u,v)+N(u,v)(5.22)G(u, v) = F(u, v) + N(u, v) \tag{5.22}G(u,v)=F(u,v)+N(
复制链接

扫一扫