最小均方誤差（維納）濾波

目標是求未污染圖像$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^{?}(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^{?}(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()