【OpenCV 例程 300 篇】106. 退化图像的逆滤波

专栏地址：『youcans 的 OpenCV 例程 300篇 - 总目录』

【第 7 章：图像复原与重建】
106. 退化图像的逆滤波
107. 退化图像的维纳滤波
108. 约束最小二乘方滤波
109. 几何均值滤波

【youcans 的 OpenCV 例程 300 篇】106. 退化图像的逆滤波

6. 退化图像复原

图像复原是对图像退化的过程进行估计，并补偿退化过程造成的失真，以便获得未经退化的原始图像或原始图像的最优估值，从而改善图像质量的一种方法。

典型的图像复原方法是根据图像退化的先验知识建立退化模型，以退化模型为基础采用滤波等手段进行处理，使复原后的图像符合一定的准则，达到改善图像质量的目的。

因此，图像复原是沿着质量降低的逆过程来重现真实的原始图像，通过去模糊函数而去除图像模糊。

6.1 退化图像的逆滤波（Inverse filter）

图像退化表示为退化算子 $\mathcal{H}$， 退化函数可以用观察法、试验法或建模法估计，则通过逆滤波就可以直接实现图像复原。用退化图像的傅里叶变换除以退化函数的傅里叶变换，得到原始图像的傅里叶变换估计：
$$\hat{F}(u,v)=\frac{G(u,v)}{H(u,v)}$$

但是，由于实际上退化图像是退化算子与加性噪声项共同作用的结果，由此得到：
$$\hat{F}(u,v)=F(u,v)+\frac{N(u,v)}{H(u,v)}$$

这表明即使获得退化函数 $H(u,v)$ 的估计，由于噪声项是未知的，因此也不能准确地复原原始图像。

进一步地，如果退化函数为 0 或很小，则噪声项的影响将非常严重（信噪比低）。这时，需要将频率限制到原点附近进行分析，可以减少遇到零值的可能性。

例程 9.20：湍流模糊退化图像的逆滤波

如前所述，通过湍流退化模型可以得到退化图像。使用该退化模型进行逆滤波，退化函数与生成退化图像所用的退化函数相反：
$$H(u,v)=e^{-k[(u-M/2{)}^2+(v-N/2{)}^2{]}^{5/6}}$$

但是，直接使用退化模型 H(u,v) 逆滤波的结果（D0=full）很差，用理想低通滤波器对退化模型 H(u,v) 在半径 D0 之外截止后，则视觉效果较好。

    # 9.20: 湍流模糊退化图像的逆滤波
    def turbulenceBlur(img, k=0.001):  # 湍流模糊传递函数: H(u,v) = exp(-k(u^2+v^2)^5/6)
        M, N = img.shape[1], img.shape[0]
        u, v = np.meshgrid(np.arange(M), np.arange(N))
        radius = (u - M//2)**2 + (v - N//2)**2
        kernel = np.exp(-k * np.power(radius, 5/6))
        return kernel

    def getDegradedImg(image, Huv, eps):  # 根据退化模型生成退化图像
        # (1) 傅里叶变换, 中心化
        fft = np.fft.fft2(image.astype(np.float32))  # 傅里叶变换
        fftShift = np.fft.fftshift(fft)  # 将低频分量移动到频域图像中心
        # (2) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数
        fftShiftFilter = fftShift * Huv  # Guv = Fuv * Huv
        # (3) 对修正傅里叶变换 进行傅里叶逆变换，逆中心化
        invShift = np.fft.ifftshift(fftShiftFilter)  # 将低频分量逆转换回图像四角
        imgIfft = np.fft.ifft2(invShift)  # 逆傅里叶变换，返回值是复数数组
        imgDegraded = np.uint8(cv2.normalize(np.abs(imgIfft), None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]
        return imgDegraded

    def ideaLPFilter(img, radius=10):  # 理想低通滤波器
        M, N = img.shape[1], img.shape[0]
        u, v = np.meshgrid(np.arange(M), np.arange(N))
        D = np.sqrt((u - M//2)**2 + (v - N//2)**2)
        kernel = np.zeros(img.shape[:2], np.float32)
        kernel[D <= radius] = 1
        return kernel

    def inverseFilter(image, Huv, D0):  # 根据退化模型逆滤波
        # (1) 傅里叶变换, 中心化
        fft = np.fft.fft2(image.astype(np.float32))  # 傅里叶变换
        fftShift = np.fft.fftshift(fft)  # 将低频分量移动到频域图像中心
        # (2) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数
        if D0==0:
            fftShiftFilter = fftShift / Huv  # Guv = Fuv / Huv
        else:
            lpFilter = ideaLPFilter(image, radius=D0)
            fftShiftFilter = fftShift / Huv * lpFilter  # Guv = Fuv / Huv
        # (3) 对修正傅里叶变换 进行傅里叶逆变换，逆中心化
        invShift = np.fft.ifftshift(fftShiftFilter)  # 将低频分量逆转换回图像四角
        imgIfft = np.fft.ifft2(invShift)  # 逆傅里叶变换，返回值是复数数组
        imgRebuild = np.uint8(cv2.normalize(np.abs(imgIfft), None, 0, 255, cv2.NORM_MINMAX))  # 归一化为 [0,255]
        return imgRebuild


    # 读取原始图像
    img = cv2.imread("../images/Fig0525a.tif", 0)  # flags=0 读取为灰度图像

    # 生成湍流模糊图像
    HTurb = turbulenceBlur(img, k=0.0025)
    imgBlur = np.abs(getDegradedImg(img, HTurb, 0.0))
    print(imgBlur.max(), imgBlur.min())

    # # 逆滤波
    imgRebuild = inverseFilter(imgBlur, HTurb, 480)  # Huv 全滤波器
    imgRebuild1 = inverseFilter(imgBlur, HTurb, D0=40)  # 在半径 D0 之外 Huv 截止
    imgRebuild2 = inverseFilter(imgBlur, HTurb, D0=70)
    imgRebuild3 = inverseFilter(imgBlur, HTurb, D0=100)

    plt.figure(figsize=(9, 7))
    plt.subplot(231), plt.title("origin"), plt.axis('off'), plt.imshow(img, 'gray')
    plt.subplot(232), plt.title("turbulence blur"), plt.axis('off'), plt.imshow(imgBlur, 'gray')
    plt.subplot(233), plt.title("inverse filter(D0=full)"), plt.axis('off'), plt.imshow(imgRebuild, 'gray')
    plt.subplot(234), plt.title("inverse filter(D0=40)"), plt.axis('off'), plt.imshow(imgRebuild1, 'gray')
    plt.subplot(235), plt.title("inverse filter(D0=70)"), plt.axis('off'), plt.imshow(imgRebuild2, 'gray')
    plt.subplot(236), plt.title("inverse filter(D0=100)"), plt.axis('off'), plt.imshow(imgRebuild3, 'gray')
    plt.tight_layout()
    plt.show()

（本节完）

---

版权声明：
youcans@xupt 原创作品，转载必须标注原文链接：(https://blog.csdn.net/youcans/article/details/123027356)
Copyright 2022 youcans, XUPT
Crated：2022-2-20