7.2 图像复原之空间滤波

图像复原（只存在噪声的复原）之空间滤波

---

---

前言

当一幅图像仅被加性噪声退化时，其退化模型可简化为：$$\begin{array}{rl}& g(x,y)=f(x,y)+\eta (x,y)\\ & G(u,v)=F(u,v)+N(u,v)\end{array}$$其中：$$\begin{array}{rl}& g(x,y)[或G(u,v)]表示退化后的图像\\ & f(x,y)[或F(u,v)]表示原图像\\ & \eta (x,y)[或N(u,v)]表示图像的噪声\end{array}$$由于噪声项通常是未知的，因此我们 不能 想当然的以为得到原图像的方式为：$$\begin{array}{rl}& f(x,y)=g(x,y)-\eta (x,y)\\ & F(u,v)=G(u,v)-N(u,v)\end{array}$$对于只存在加性噪声的图像，可以使用空间滤波的方法来估计 $f(x,y)$

1. 均值滤波器

1.1 算术平均滤波器

算术平均滤波器是最简单的均值滤波器，与空间域滤波中的盒式滤波器完全相同。令：$$S_{xy}表示中心为(x,y)、大小为m\times n的邻域$$则复原的图像 ：$$\hat{f}(x,y)=\frac{1}{mn}\sum _{(r,c)\in S_{xy}}g(r,c)$$用白话来说就是使用 $(x,y)$ 的 $m\times 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_{xy}表示中心为(x,y)、大小为m\times n的邻域$$则复原的图像 ：$$\hat{f}(x,y)={\left[\prod _{(r,c)\in S_{xy}}g(r,c)\right]}^{\frac{1}{mn}}$$简单来说几何均值就是先求取中心点 $(x,y)$ 的 $m\times n$ 领域的像素之积，然后对像素之积求 $\frac{1}{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

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_{xy}表示中心为(x,y)、大小为m\times n的邻域$$则复原的图像 ：$$\hat{f}(x,y)=\frac{mn}{{\sum }_{(r,c)\in S_{xy}}\frac{1}{g(r,c)}}$$简单解释一下。

• 第一步：用 $1$ 除以领域内每个像素的值，即 $\frac{1}{g(r,c)}$ 。
• 第二步：对第一步的结果求和，即 ${\sum }_{(r,c)\in S_{xy}}\frac{1}{g(r,c)}$。
• 第三步：用 $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_{xy}表示中心为(x,y)、大小为m\times n的邻域$$则复原的图像 ：$$\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}$$简单解释一下，式中：

• $Q$ 称为滤波器的阶数，这种滤波器适用于降低或消除噪声 。
• $Q$ 为正值时，该滤波器消除胡椒噪声。
• $Q$ 为负值时，该滤波器消除盐粒噪声。
• 然而，该滤波器不能同时消除这两种噪声。
• $Q=0$ 时，反谐波平均滤波器简化为算数平均滤波器。
• $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$ 值）	灵活去噪能力强