4.2 Python图像的图像恢复-组合滤波器

4.2 Python图像的图像恢复-组合滤波器

1 算法原理

1.1 混合滤波器

本文以中值滤波+巴特沃斯低通滤波+同态滤波混合滤波器为例，读者可以自由组合其他滤波器。

本混合滤波器算法设计思路：先对图片添加椒盐噪声，先采用中值滤波器进行去噪处理，然后采用巴特沃斯低通滤波和同态滤波混合滤波器进行图像处理。巴特沃斯低通滤波器和同态滤波器在上一节已经介绍过其原理，这里不再重复给出。

中值滤波–cv2.medianBlur(src, ksize)。这里的核大小 ksize 必须是奇数，将该点周围的像素点包括本身，按次序排列，取中位数作为点的像素值。中值滤波器的特点为：随着核大小逐渐变大，会让图像变得更加模糊；核必须是大于 1 的奇数，如 3、5、7 等；在代码 dst = cv2.medianBlur(src, ksize) 中 填写核大小时，只需填写一个数即可，如 3、5、7 ，9等，对比均值滤波函数用法。

1.2 选择性滤波器

选择性滤波主要有两类：

带阻滤波器或带通滤波器； 陷波滤波器。陷波滤波器是更有用的选择性滤波器可以有效的去除周期性噪声。陷波滤波器可以用高通滤波器的乘积来构造。一般形式为：

其中，Hk(u,v)和 H−k(u,v)是高通滤波器。距离的计算：

和：

本次实验采用陷波滤波器作为选择性滤波器的算法实现，它包含两个陷波对：

2 代码

运行代码说明

1.要改变代码中的图片地址（地址不能有中文）

更改put(path)函数中的路径put(r'../image/image1.jpg')

2.注意最后的plt.savefig('1.new.jpg')是保存plt图像，如果不使用可以注释掉

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# author：mgboy time:2021/6/26

import os
import numpy as np
import cv2
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 巴特沃斯低通滤波器
def ButterworthLowPassFilter(image, d, n, s1):
    """
    Butterworth低通滤波器
    """
    f = np.fft.fft2(image)
    fshift = np.fft.fftshift(f)

    def make_transform_matrix(d):
        transform_matrix = np.zeros(image.shape)
        center_point = tuple(map(lambda x: (x - 1) / 2, s1.shape))
        for i in range(transform_matrix.shape[0]):
            for j in range(transform_matrix.shape[1]):

                def cal_distance(pa, pb):
                    from math import sqrt

                    dis = sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
                    return dis
                dis = cal_distance(center_point, (i, j))
                transform_matrix[i, j] = 1 / (1 + (dis / d) ** (2 * n))
        return transform_matrix


    d_matrix = make_transform_matrix(d)
    new_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift * d_matrix)))
    return new_img

# 同态滤波器
def homomorphic_filter(src, d0=10, r1=0.5, rh=2, c=4, h=2.0, l=0.5):
    gray = src.copy()
    if len(src.shape) > 2:
        gray = cv2.cvtColor(src, cv2.COLOR_BGR2GRAY)
    gray = np.float64(gray)
    rows, cols = gray.shape

    gray_fft = np.fft.fft2(gray)
    gray_fftshift = np.fft.fftshift(gray_fft)
    dst_fftshift = np.zeros_like(gray_fftshift)
    M, N = np.meshgrid(np.arange(-cols // 2, cols // 2), np.arange(-rows // 2, rows // 2))
    D = np.sqrt(M ** 2 + N ** 2)
    Z = (rh - r1) * (1 - np.exp(-c * (D ** 2 / d0 ** 2))) + r1
    dst_fftshift = Z * gray_fftshift
    dst_fftshift = (h - l) * dst_fftshift + l
    dst_ifftshift = np.fft.ifftshift(dst_fftshift)
    dst_ifft = np.fft.ifft2(dst_ifftshift)
    dst = np.real(dst_ifft)
    dst = np.uint8(np.clip(dst, 0, 255))
    return dst

def salt_pepper(image, salt, pepper):
    """
    添加椒盐噪声的图像
    :param image: 输入图像
    :param salt: 盐比例
    :param pepper: 椒比例
    :return: 添加了椒盐噪声的图像
    """
    height = image.shape[0]
    width = image.shape[1]
    pertotal = salt + pepper    #总噪声占比
    noise_image = image.copy()
    noise_num = int(pertotal * height * width)
    for i in range(noise_num):
        rows = np.random.randint(0, height-1)
        cols = np.random.randint(0,width-1)
        if(np.random.randint(0,100)<salt*100):
            noise_image[rows][cols] = 255
        else:
            noise_image[rows][cols] = 0
    return noise_image

def put(path):
    img = cv2.imread(path, 1)
    # img = cv2.imread(os.path.join(base, path), 1)
    img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    # 傅里叶变换
    f = np.fft.fft2(img)
    fshift = np.fft.fftshift(f)
    # 取绝对值后将复数变化为实数，取对数的目的是将数据变换到0~255
    s1 = np.log(np.abs(fshift))
    # 加入椒盐噪声
    noise = salt_pepper(img, 0.05, 0.05)
    # 经过中值滤波器消除噪声
    res = cv2.medianBlur(noise, 5)
    # 图像经过巴特沃斯低通滤波器处理
    # 设置D=100,N=1
    res = ButterworthLowPassFilter(res, 100, 1,s1)
    # 处理过的图像再经过同态滤波处理
    res = homomorphic_filter(res, 100, 1)

    res = ButterworthLowPassFilter(res, 100, 1,s1)

    f = np.fft.fft2(noise)
    fshift = np.fft.fftshift(f)
    w, h = img.shape
    """设计理想陷波滤波器"""
    flt = np.zeros(noise.shape)
    rx1 = w / 4
    ry1 = h / 2
    rx2 = w * 3 / 4
    ry2 = h / 2;
    r = min(w, h) / 6  # 半径
    for i in range(1, w):
        for j in range(1, h):
            if ((i - rx1) ** 2 + (j - ry1) ** 2 >= r ** 2) and ((i - rx2) ** 2 + (j - ry2) ** 2 >= r ** 2):
                flt[i, j] = 1
    # 选择滤波
    new_img = np.abs(np.fft.ifft2(np.fft.ifftshift(fshift * flt)))

    plt.subplot(231)
    plt.axis('off')
    plt.title('原始图像')
    plt.imshow(img, cmap='gray')
    plt.subplot(232)
    plt.axis('off')
    plt.title('噪声图像')
    plt.imshow(noise, cmap='gray')
    plt.subplot(233)
    plt.axis('off')
    plt.title('混合滤波器')
    plt.imshow(res, cmap='gray')
    plt.subplot(234),plt.axis('off'),plt.imshow(s1,'gray'),plt.title('中心频率域')
    plt.subplot(235), plt.axis('off'), plt.imshow(flt, 'gray'), plt.title('陷波滤波器')
    plt.subplot(236),plt.axis('off'),plt.imshow(new_img,'gray'),plt.title('选择性滤波器')

    # plt.savefig('2.new.jpg')
    plt.show()


# 图像处理函数，要传入路径
put(r'../image/image3.jpg')

3 效果