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')

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144

3 效果