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

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

误差函数的最小值在频率域中的表达比如下：
$$\begin{array}{rl}\hat{F}(u,v)& =\left[\frac{H^{\ast }(u,v)S_f(u,v)}{S_f(u,v)|H(u,v){|}^2+S_{\eta }(u,v)}\right]G(u,v)\\ & =\left[\frac{H^{\ast }(u,v)}{|H(u,v){|}^2+S_{\eta }(u,v)/S_f(u,v)}\right]G(u,v)\\ & =\left[\frac{1}{H(u,v)}\frac{|H(u,v){|}^2}{|H(u,v{)}^2+K}\right]G(u,v)\end{array}\text{\hspace{1em}(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_tanif 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)]=1return PSF / PSF.sum()  #对点扩散函数进行归一化亮度else:for i in range(degree):offset=round(i*slope_cot)PSF[int(center_position-offset),int(center_position+offset)]=1return PSF / PSF.sum()def get_motion_dsf(image_size, motion_angle, motion_dis):"""Get motion PSFparam: image_size: input image shapeparam: motion_angle: blur motion angleparam: motion_dis: blur distant, the greater value, more blurredreturn normalize PSF"""PSF = np.zeros(image_size)  # 点扩散函数x_center = (image_size[0] - 1) / 2y_center = (image_size[1] - 1) / 2sin_val = np.sin(motion_angle * np.pi / 180)cos_val = np.cos(motion_angle * np.pi / 180)# 将对应角度上motion_dis个点置成1for 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)] = 1return PSF / PSF.sum()    # 归一化# 对图片进行运动模糊
def make_blurred(input, PSF, eps):"""blurred image with PSFparam: input: input imageparam: PSF: input PSF maskparam: eps: epsilon, very small value, to make sure not divided or multiplied by zeroreturn blurred image"""input_fft = np.fft.fft2(input)              #  image FFTPSF_fft = np.fft.fft2(PSF)+ eps             # PSF FFT plus epsilonblurred = np.fft.ifft2(input_fft * PSF_fft) # image FFT multiply PSF FFTblurred = np.abs(np.fft.fftshift(blurred))return blurreddef inverse_filter(input, PSF, eps): """inverse filter using FFT to denoiseparam: input: input imageparam: PSF: known PSFparam: eps: epsilon"""input_fft = np.fft.fft2(input)PSF_fft = np.fft.fft2(PSF) + eps           #噪声功率，这是已知的，考虑epsilonresult = np.fft.ifft2(input_fft / PSF_fft) #计算F(u,v)的傅里叶反变换result = np.abs(np.fft.fftshift(result))return resultdef wiener_filter(input, PSF, eps, K=0.01):"""wiener filter for image denoiseparam: input: input imageparam: PSF: input the PSF maskparam: eps: epsilonparam: K=0.01: K value for wiener fuctionreturn image after wiener filter"""input_fft = np.fft.fft2(input)PSF_fft = np.fft.fft2(PSF) + epsPSF_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()