標(biāo)題

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

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

目標(biāo)是求未污染圖像$f$的一個(gè)估計(jì)$\hat{f}$，使它們之間的均方誤差最小。
$$e^2=E\{(f?\hat{f}{)}^2\}\text{\hspace{1em}(5.80)}$$

誤差函數(shù)的最小值在頻率域中的表達(dá)比如下：
$$\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)}$$

注：逆濾波與維納濾波都要求未退化圖像和噪聲的功率譜是已知的。

# 運(yùn)動(dòng)模糊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() # 仿真運(yùn)動(dòng)模糊 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() #對(duì)點(diǎn)擴(kuò)散函數(shù)進(jìn)行歸一化亮度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) # 點(diǎn)擴(kuò)散函數(shù)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)# 將對(duì)應(yīng)角度上motion_dis個(gè)點(diǎn)置成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() # 歸一化# 對(duì)圖片進(jìn)行運(yùn)動(dòng)模糊 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) #計(jì)算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# 要實(shí)現(xiàn)的功能都在這里調(diào)用 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))# 進(jìn)行運(yùn)動(dòng)模糊處理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產(chǎn)生隨機(jī)的函數(shù)blurred_noisy = blurred + 0.1 * blurred.std() * np.random.standard_normal(blurred.shape) # 顯示添加噪聲且運(yùn)動(dòng)模糊的圖像plt.subplot(234), plt.imshow(blurred_noisy, 'gray'), plt.title("motion & noisy blurred")plt.xticks([]), plt.yticks([])# 對(duì)添加噪聲的圖像進(jìn)行逆濾波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([])# 對(duì)添加噪聲的圖像進(jìn)行維納濾波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()

# 要實(shí)現(xiàn)的功能都在這里調(diào)用 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))# 進(jìn)行運(yùn)動(dòng)模糊處理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產(chǎn)生隨機(jī)的函數(shù)blurred_noisy = blurred + 0.1 * blurred.std() * np.random.standard_normal(blurred.shape) # 顯示添加噪聲且運(yùn)動(dòng)模糊的圖像plt.subplot(234), plt.imshow(blurred_noisy, 'gray'), plt.title("motion & noisy blurred")plt.xticks([]), plt.yticks([])# 對(duì)添加噪聲的圖像進(jìn)行逆濾波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([])# 對(duì)添加噪聲的圖像進(jìn)行維納濾波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()

總結(jié)

以上是生活随笔為你收集整理的第5章 Python 数字图像处理(DIP) - 图像复原与重建15 - 最小均方误差（维纳）滤波的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。