数字图像处理第五章

别来这个网址

已于 2023-05-05 23:26:43 修改

阅读量497

点赞数

分类专栏：数字图像处理文章标签： opencv 计算机视觉 python

于 2023-05-05 23:03:53 首次发布

本文链接：https://blog.csdn.net/m0_68784427/article/details/130500765

版权

数字图像处理专栏收录该内容

10 篇文章 2 订阅

订阅专栏

数字图像处理第五章

图像退化/复原处理的模型

用退化函数对退化过程建模，它和附加噪声选项一起，作用于输入图像f(x,y),产生一幅退化的图像g(x,y): 在这里插入图片描述
其中 $\eta(x,y)$ h为乘性噪声，η为加性噪声

噪声模型

比较重要的噪声：高斯噪声，瑞利噪声，爱尔兰噪声，指数噪声，均匀噪声，椒盐噪声
其直方图如下：在这里插入图片描述
针对上述噪声，对于高斯噪声做了测试，代码如下

import cv2
import numpy as np


# 读取图像images/lena_gray_512.tif
img = cv2.imread('images/lena_gray_512.tif', cv2.IMREAD_GRAYSCALE)

# 添加高斯噪声
img = img/255
noise = np.random.normal(0, 0.1, img.shape)
img_noise = img + noise

# 显示图像
cv2.imshow('img', img)
cv2.imshow('img_noise', img_noise)
cv2.waitKey()
cv2.destroyAllWindows()

运行结果如图：
在这里插入图片描述

值得一提的是，由于cv2.imshow只允许 (0-1 float or 0-255 int) 的输入，因此必须进行转换才能得到正确的显示结果，更详细的原因参见：https://blog.csdn.net/howlclat/article/details/107216722

空间滤波

前文已经介绍了很多种空间滤波器，下面我们将介绍一种自适应中值滤波器，并介绍其与传统中值滤波器的区别。

中值滤波器是一种常用的去除噪声的滤波器，它的原理是在图像的每一个像素点上取一个邻域，将邻域中的像素按照亮度值排序，取中间值作为该像素点的新值。这种方法可以有效地去除椒盐噪声等离散噪声，但是对于连续性噪声（如高斯噪声）的效果不是很好。

自适应中值滤波器是一种改进的中值滤波器，它能够适应不同的噪声强度，并选择合适的邻域大小进行滤波。在自适应中值滤波器中，不同于传统的中值滤波器使用固定的邻域大小，它会动态地调整邻域大小，以使滤波器更好地适应噪声的不同强度和尺度。同时，它还引入了两个阈值，用于检测和处理脉冲噪声。

具体来说，在自适应中值滤波器中，算法会根据像素点周围像素的亮度值差异以及像素点本身的亮度值来动态调整邻域大小，并且只有在邻域中存在脉冲噪声时才进行处理。因此，自适应中值滤波器可以更好地保留图像细节和边缘信息，同时对于不同类型的噪声也具有更好的去噪效果。

代码如下：

def adaptive_median_filter(img, ksize_max):
    # 自适应中值滤波
    img_new = np.zeros(img.shape)
    img_pad = np.pad(img, ((ksize_max // 2, ksize_max // 2), (ksize_max // 2, ksize_max // 2)), 'constant', constant_values=0)
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            ksize = 3
            while True:
                img_pad_tmp = img_pad[i:i + ksize, j:j + ksize]
                z_min = np.min(img_pad_tmp)
                z_max = np.max(img_pad_tmp)
                z_med = np.median(img_pad_tmp)
                A1 = z_med - z_min
                A2 = z_med - z_max
                if A1 > 0 and A2 < 0:
                    B1 = img_pad[i, j] - z_min
                    B2 = img_pad[i, j] - z_max
                    if B1 > 0 and B2 < 0:
                        img_new[i, j] = img_pad[i, j]
                    else:
                        img_new[i, j] = z_med
                    break
                else:
                    ksize += 2
                    if ksize > ksize_max:
                        img_new[i, j] = z_med
                        break
    return img_new

运行结果：
在这里插入图片描述
可以看到得到了不错的结果。

频域滤波削减周期噪声

用频率域技术可以有效地分析并滤除周期噪声，为消除这些周期噪声，将使用三种类型的选择性滤波器：带阻、带通和陷波滤波器。

带阻滤波器：它的主要应用之一是在频率域噪声成分的一般位置近似一直的应用中消除噪声。
带通滤波器：它执行与带阻滤波器相反的操作。
陷波滤波器：阻止（或通过）事先定义的中心频率邻域内的频率。

线性、位置不变的退化

在退化复原模型中，输入输出关系可以表示为： $H[f(x,y)+\eta(x,y)]$

如果系统H是一个线性系统，那么H应该满足可加性和均匀性：
$H[f_1(x,y)+f_2(x,y)]=H[f_1(x,y)]+H[f_2(x,y)]$
$H[af_1(x,y)] = aH[f_1(x,y)]$
其中a为比例常数，f1与f2为任意两个输入图像。
总之，具有加性噪声的线性空间不变退化系统，可在空间域建模为退化（点扩散）函数与一幅图像的卷积，然后再加上噪声。基于卷积定理，在频率域中，同样的过程可表示为图像和退化函数的变换的乘积，然后再加上噪声的变换。

许多类型的退化可近似为线性、位置不变的过程。这种方法的优点是：可以使用许多线性系统理论的工具来解决图像复原问题。

估计退化函数

在阁像复原中，主要有三种估计退化函数的方法：观察法，试验法，数学建模法。

观察法

假设提供了一幅退化图像，而没有退化函数H的知识，那么估计该函数的一个方法就是收集图像自身的信息。例如，如果图像是模糊的，可以观察包含简单节后的一小部分图像，像某一物体和背景的一部分。为了减少观察时的噪声影响，可以寻找强信号的内容区。
例如，如果我们发现不够清晰，可以去降噪，锐化等等一系列操作

试验法

试验估计是一种通过实验来估计图像退化函数的方法。它通过对已知的测试图像施加已知的退化操作，比较输入图像和输出图像的误差或距离，从而估计退化函数。尽管该方法简单易行，但其结果受测试图像和退化操作的影响，精度可能不如其他更复杂的方法。
例如，我们要对a图片复原，已知a是通过路由器传输的，那么我们可以让图像B也用同一个路由器传输得到图像b。然后根据B和b的区别来还原a。

数学建模法

看起来就是对特定场景进行建模，比如高速运动以及大气湍流场景下拍摄的图片，可以通过对其进行建模来还原拍摄到的图像

其他滤波

逆滤波

即根据滤波函数求出逆滤波函数并据此变换。但是显然易得的是，我们无法完全还原原本的图像。

维也纳滤波

维纳滤波也叫最小均方误差滤波，是一种建立在最小化统计准则的基础上的复原方法，在平均意义上，它可以看成是最优的。维纳滤波综合了退化函数和噪声统计特征两个方面进行复原处理，在认为图像和噪声是随机过程的基础上，以恢复图像和原图像的均方误差最小为准则。

from scipy.signal import wiener
import cv2
import numpy as np
import matplotlib.pyplot as plt


def gasuss_noise(image, mean=0, var=0.001):
    image = np.array(image / 255, dtype=float)
    noise = np.random.normal(mean, var ** 0.5, image.shape)
    out = image + noise
    if out.min() < 0:
        low_clip = -1.
    else:
        low_clip = 0.
    out = np.clip(out, low_clip, 1.0)
    out = np.uint8(out * 255)
    return out


if __name__ == '__main__':
    lena = cv2.imread('images/lena_gray_512.tif', cv2.IMREAD_GRAYSCALE)

    if lena.shape[-1] == 3:
        lenaGray = cv2.cvtColor(lena, cv2.COLOR_BGR2GRAY)
    else:
        lenaGray = lena.copy()

    # 添加高斯噪声
    lenaNoise = gasuss_noise(lenaGray)

    # 维纳滤波
    lenaNoise = lenaNoise.astype('float64')
    lenaWiener = wiener(lenaNoise, [3, 3])
    lenaWiener = np.uint8(lenaWiener / lenaWiener.max() * 255)

    fig, axs = plt.subplots(1, 3, figsize=(12, 4))

    # 显示原图
    axs[0].imshow(lenaGray, cmap='gray')
    axs[0].set_title('original')

    # 显示添加高斯噪声后的图像
    axs[1].imshow(lenaNoise, cmap='gray')
    axs[1].set_title('aftergasuss')

    # 显示经过维纳滤波后的图像
    axs[2].imshow(lenaWiener, cmap='gray')
    axs[2].set_title('afterwiener')

    plt.show()

对比
在这里插入图片描述

由投影重建图像

当收集到一副图像各个角度的投影后，并希望通过这些投影的图像重建原图像。通过下图也就是直接重建的图像可见，直接重建会有非常明显的“晕环”现象。在这里插入图片描述

雷登变换/傅里叶变换

在这里插入图片描述
雷登变换示例

图像在各个角度下的雷登变换的集合称为正弦图。
以第一行的图像为例。明显的，其像素点的值在90°的投影下是最多的，在0°或180°的投影下是最少的。对应于右侧的正弦图也能够体现出来。

我们可以把各个角度下的投影经过一次反投影然后求和复原图像，得到下图
在这里插入图片描述
我们可以注意到上图有明显的晕环，虽然增大采样次数可以改善，但是会影响计算速度，傅里叶切片定理可以解决这个问题

傅里叶切片定理

在这里插入图片描述
傅里叶切片定理用一句话表示就是：一个投影的一维傅里叶变换就是得到该投影原图的二维傅里叶变换的一个切片，其切片角度就是投影的角度θ。
由此可以消除晕环，这种重建方法称为滤波反投影法。

总结步骤如下：
1.计算每个投影的一维傅里叶变换
2.用一个滤波函数|w|乘以每个傅里叶变换，就是加窗。
3.得到每个滤波后的一维反傅里叶变换。
4.将步骤3得到的结果求和

代码如下：

from scipy import ndimage
import numpy as np
import cv2
import matplotlib.pyplot as plt


def discreteRadonTransform(image, steps):  # 离散雷登变换
    channels = image.shape[0]
    resRadon = np.zeros((channels, channels), dtype=np.float32)
    for s in range(steps):
        rotation = ndimage.rotate(image, -s * 180 / steps, reshape=False).astype(np.float32)
        resRadon[:, s] = sum(rotation)
    return resRadon


def inverseRadonTransform(image, steps):  # 雷登变换反投影重建图像
    channels = image.shape[0]
    backproject = np.zeros((steps, channels, channels))
    for s in range(steps):
        expandDims = np.expand_dims(image[:, s], axis=0)
        repeat = expandDims.repeat(channels, axis=0)
        backproject[s] = ndimage.rotate(repeat, s * 180 / steps, reshape=False).astype(np.float32)
    invRadon = np.sum(backproject, axis=0)
    return invRadon


def SLFilter(N, d):  # SL 滤波器, Sinc 函数对斜坡滤波器进行截断
    filterSL = np.zeros((N,))
    for i in range(N):
        filterSL[i] = - 2 / (np.pi ** 2 * d ** 2 * (4 * (i - N / 2) ** 2 - 1))
    return filterSL


def filterInvRadonTransform(image, steps):  # 滤波反投影重建图像
    channels = len(image[0])
    backproject = np.zeros((steps, channels, channels))  # 反投影
    filter = SLFilter(channels, 1)  # SL 滤波器
    for s in range(steps):
        value = image[:, s]  # 投影值
        valueFiltered = np.convolve(filter, value, "same")  # 投影值和 SL 滤波器进行卷积
        filterExpandDims = np.expand_dims(valueFiltered, axis=0)
        filterRepeat = filterExpandDims.repeat(channels, axis=0)
        backproject[s] = ndimage.rotate(filterRepeat, s * 180 / steps, reshape=False).astype(np.float32)
    filterInvRadon = np.sum(backproject, axis=0)
    return filterInvRadon


# 读取原始图像
img1 = cv2.imread("ch5/1.tif", 0)  # flags=0 读取为灰度图像
img2 = cv2.imread("ch5/2.tif", 0)

# 雷登变换
imgRadon1 = discreteRadonTransform(img1, img1.shape[0])  # Radon 变换
imgRadon2 = discreteRadonTransform(img2, img2.shape[0])

# 雷登变换反投影重建图像
imgInvRadon1 = inverseRadonTransform(imgRadon1, imgRadon1.shape[0])
imgInvRadon2 = inverseRadonTransform(imgRadon2, imgRadon2.shape[0])

# 滤波反投影重建图像
imgFilterInvRadon1 = filterInvRadonTransform(imgRadon1, imgRadon1.shape[0])
imgFilterInvRadon2 = filterInvRadonTransform(imgRadon2, imgRadon2.shape[0])

plt.figure(figsize=(9, 7))
plt.subplot(231), plt.axis('off'), plt.title("origin image"), plt.imshow(img1, 'gray')  # 绘制原始图像
plt.subplot(232), plt.axis('off'), plt.title("inverse Radon trans"), plt.imshow(imgInvRadon1, 'gray')
plt.subplot(233), plt.axis('off'), plt.title("filter inv-Radon trans"), plt.imshow(imgFilterInvRadon1, 'gray')
plt.subplot(234), plt.axis('off'), plt.title("origin image"), plt.imshow(img2, 'gray')
plt.subplot(235), plt.axis('off'), plt.title("inverse Radon trans"), plt.imshow(imgInvRadon2, 'gray')
plt.subplot(236), plt.axis('off'), plt.title("filter inv-Radon trans"), plt.imshow(imgFilterInvRadon2, 'gray')
plt.tight_layout()
plt.show()

运行结果如下
在这里插入图片描述