组合滤波器

import os
import cv2
import numpy as np
from math import sqrt
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False
def ButterworthLowPassFilter(image,d,n,s1):
    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):
                    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):
    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.cvtColor(img,cv2.COLOR_BGR2GRAY)
    #傅里叶变换
    f=np.fft.fft2(img)
    fshift=np.fft.fftshift(f)
    #取绝对值后将复数变化为实数
    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('C:/Users/xpp/Desktop/result.png')
    plt.show()
put(r'C:/Users/xpp/Desktop/Lena.png')