第5章 Python 数字图像处理(DIP) - 图像复原与重建15 - 最小均方误差（维纳）滤波

最新推荐文章于 2022-11-10 07:37:18 发布

jasneik

最新推荐文章于 2022-11-10 07:37:18 发布

阅读量1.5k

点赞数 7

分类专栏： # 第5章图像复原与重建数字图像处理(DIP)--Python 文章标签：图像处理 python opencv numpy 图像识别

本文链接：https://blog.csdn.net/jasneik/article/details/115797426

版权

数字图像处理(DIP)--Python 同时被 2 个专栏收录

57 篇文章 129 订阅

订阅专栏

第5章图像复原与重建

17 篇文章 22 订阅

订阅专栏

标题

- 最小均方误差（维纳）滤波

最小均方误差（维纳）滤波

目标是求未污染图像 $f$ 的一个估计 $\hat{f}$ ，使它们之间的均方误差最小。
$e^2 = E \big\{(f - \hat{f})^2 \big\} \tag{5.80}$

误差函数的最小值在频率域中的表达比如下：
$\begin{aligned} \hat{F}(u, v) & = \Bigg[\frac{H^*(u, v) S_{f}(u, v)}{S_{f}(u, v)|H(u, v)|^2 + S_{\eta}(u, v)} \Bigg] G(u, v) \\ & = \Bigg[\frac{H^*(u, v) }{|H(u, v)|^2 + S_{\eta}(u, v) / S_{f}(u, v)} \Bigg] G(u, v) \\ & = \Bigg[\frac{1}{H(u,v)} \frac{|H(u,v)|^2}{|H(u,v)^2 + K} \Bigg]G(u,v) \end{aligned} \tag{5.81}$

注：逆滤波与维纳滤波都要求未退化图像和噪声的功率谱是已知的。

# 运动模糊PSF与谱
fft_shift = np.fft.fftshift(PSF)
fft = np.fft.fft2(PSF)
spect = spectrum_fft(fft)
plt.figure(figsize=(8, 8))
plt.subplot(1,2,1), plt.imshow(PSF, 'gray')
plt.subplot(1,2,2), plt.imshow(spect, 'gray')
plt.show()

# 仿真运动模糊
def motion_process(image_size, motion_angle, degree=15):
    """
    This function has some problem
    """
    PSF = np.zeros(image_size)
#     print(image_size)
    center_position=(image_size[0]-1)/2
#     print(center_position)
 
    slope_tan=math.tan(motion_angle*math.pi/180)
    slope_cot=1/slope_tan
    if slope_tan<=1:
        for i in range(degree):
            offset=round(i*slope_tan)    #((center_position-i)*slope_tan)
            PSF[int(center_position+offset),int(center_position-offset)]=1
        return PSF / PSF.sum()  #对点扩散函数进行归一化亮度
    else:
        for i in range(degree):
            offset=round(i*slope_cot)
            PSF[int(center_position-offset),int(center_position+offset)]=1
        return PSF / PSF.sum()

def get_motion_dsf(image_size, motion_angle, motion_dis):
    """
    Get motion PSF
    param: image_size: input image shape
    param: motion_angle: blur motion angle
    param: motion_dis: blur distant, the greater value, more blurred
    return normalize PSF
    """
    
    PSF = np.zeros(image_size)  # 点扩散函数
    x_center = (image_size[0] - 1) / 2
    y_center = (image_size[1] - 1) / 2
 
    sin_val = np.sin(motion_angle * np.pi / 180)
    cos_val = np.cos(motion_angle * np.pi / 180)
 
    # 将对应角度上motion_dis个点置成1
    for i in range(motion_dis):
        x_offset = round(sin_val * i)
        y_offset = round(cos_val * i)
        PSF[int(x_center - x_offset), int(y_center + y_offset)] = 1
 
    return PSF / PSF.sum()    # 归一化
    
# 对图片进行运动模糊
def make_blurred(input, PSF, eps):
    """
    blurred image with PSF
    param: input: input image
    param: PSF: input PSF mask
    param: eps: epsilon, very small value, to make sure not divided or multiplied by zero
    return blurred image
    """
    input_fft = np.fft.fft2(input)              #  image FFT
    PSF_fft = np.fft.fft2(PSF)+ eps             # PSF FFT plus epsilon
    blurred = np.fft.ifft2(input_fft * PSF_fft) # image FFT multiply PSF FFT
    blurred = np.abs(np.fft.fftshift(blurred))
    return blurred
 
def inverse_filter(input, PSF, eps): 
    """
    inverse filter using FFT to denoise
    param: input: input image
    param: PSF: known PSF
    param: eps: epsilon
    """
    input_fft = np.fft.fft2(input)
    PSF_fft = np.fft.fft2(PSF) + eps           #噪声功率，这是已知的，考虑epsilon
    result = np.fft.ifft2(input_fft / PSF_fft) #计算F(u,v)的傅里叶反变换
    result = np.abs(np.fft.fftshift(result))
    return result
 
def wiener_filter(input, PSF, eps, K=0.01):
    """
    wiener filter for image denoise
    param: input: input image
    param: PSF: input the PSF mask
    param: eps: epsilon
    param: K=0.01: K value for wiener fuction
    return image after wiener filter
    """
    input_fft = np.fft.fft2(input)
    PSF_fft = np.fft.fft2(PSF) + eps
    PSF_fft_1 = np.conj(PSF_fft) / (np.abs(PSF_fft)**2 + K)
    # 按公式，居然得不到正确的值
#     PSF_abs = PSF * np.conj(PSF)
#     PSF_fft_1 = (1 / (PSF + eps)) * (PSF_abs / (PSF_abs + K))
    result = np.fft.ifft2(input_fft * PSF_fft_1)
    result = np.abs(np.fft.fftshift(result))
    return result

# 要实现的功能都在这里调用
if __name__ == "__main__":
    
    image = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0526(a)(original_DIP).tif', 0)

    # 显示原图像
    plt.figure(1, figsize=(6, 6))
    plt.title("Original Image"), plt.imshow(image, 'gray')
    plt.xticks([]), plt.yticks([])

    plt.figure(2, figsize=(18, 12))
    
    # 进行运动模糊处理
    PSF = get_motion_dsf(image.shape[:2], -50, 100)
    blurred = make_blurred(image, PSF, 1e-3)
    plt.subplot(231), plt.imshow(blurred, 'gray'), plt.title("Motion blurred")
    plt.xticks([]), plt.yticks([])

    # 逆滤波
    result = inverse_filter(blurred, PSF, 1e-3)   
    plt.subplot(232), plt.imshow(result, 'gray'), plt.title("inverse deblurred")
    plt.xticks([]), plt.yticks([])

    # 维纳滤波
    result = wiener_filter(blurred, PSF, 1e-3)     
    plt.subplot(233), plt.imshow(result, 'gray'), plt.title("wiener deblurred(k=0.01)")
    plt.xticks([]), plt.yticks([])

    # 添加噪声,standard_normal产生随机的函数
    blurred_noisy = blurred + 0.1 * blurred.std() * np.random.standard_normal(blurred.shape)   

    # 显示添加噪声且运动模糊的图像
    plt.subplot(234), plt.imshow(blurred_noisy, 'gray'), plt.title("motion & noisy blurred")
    plt.xticks([]), plt.yticks([])

    # 对添加噪声的图像进行逆滤波
    result = inverse_filter(blurred_noisy, PSF, 0.1 + 1e-3)    
    plt.subplot(235), plt.imshow(result, 'gray'), plt.title("inverse deblurred")
    plt.xticks([]), plt.yticks([])

    # 对添加噪声的图像进行维纳滤波
    result = wiener_filter(blurred_noisy, PSF, 0.1 + 1e-3)         
    plt.subplot(236), plt.imshow(result, 'gray'), plt.title("wiener deblurred(k=0.01)")
    plt.xticks([]), plt.yticks([])

    plt.tight_layout()
    plt.show()

在这里插入图片描述

# 要实现的功能都在这里调用
if __name__ == "__main__":
    
    image = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0526(a)(original_DIP).tif', 0)

    # 显示原图像
    plt.figure(1, figsize=(6, 6))
    plt.title("Original Image"), plt.imshow(image, 'gray')
    plt.xticks([]), plt.yticks([])
    
    plt.figure(2, figsize=(18, 12))
    
    # 进行运动模糊处理
    PSF = get_motion_dsf(image.shape[:2], -45, 95)
    blurred = cv2.imread('DIP_Figures/DIP3E_Original_Images_CH05/Fig0529(d)(medium_noise_var_pt01).tif', 0)
    plt.subplot(231), plt.imshow(blurred, 'gray'), plt.title("Motion blurred")
    plt.xticks([]), plt.yticks([])

    # 逆滤波
    result = inverse_filter(blurred, PSF, 1e-3)   
    plt.subplot(232), plt.imshow(result, 'gray'), plt.title("inverse deblurred")
    plt.xticks([]), plt.yticks([])

    # 维纳滤波
    result = wiener_filter(blurred, PSF, 1e-3, K=0.01)     
    plt.subplot(233), plt.imshow(result, 'gray'), plt.title("wiener deblurred(k=0.01)")
    plt.xticks([]), plt.yticks([])

    # 添加噪声,standard_normal产生随机的函数
    blurred_noisy = blurred + 0.1 * blurred.std() * np.random.standard_normal(blurred.shape)   

    # 显示添加噪声且运动模糊的图像
    plt.subplot(234), plt.imshow(blurred_noisy, 'gray'), plt.title("motion & noisy blurred")
    plt.xticks([]), plt.yticks([])

    # 对添加噪声的图像进行逆滤波
    result = inverse_filter(blurred_noisy, PSF, 0.1 + 1e-3)    
    plt.subplot(235), plt.imshow(result, 'gray'), plt.title("inverse deblurred")
    plt.xticks([]), plt.yticks([])

    # 对添加噪声的图像进行维纳滤波
    result = wiener_filter(blurred_noisy, PSF, 0.1 + 1e-3)         
    plt.subplot(236), plt.imshow(result, 'gray'), plt.title("wiener deblurred(k=0.01)")
    plt.xticks([]), plt.yticks([])

    plt.tight_layout()
    plt.show()