1.運動模糊
為了模擬圖像退化的過程，在這里創(chuàng)建了一個用于模擬運動模糊的點擴散函數(shù)，具體模糊的方向取決于輸入的motion_angle。如果運動方向接近水平，則模糊效果近似水平，如果運動方向接近垂直，則模糊效果近似垂直。具體操作如下。首先創(chuàng)建二維數(shù)組PSF，并將其所有元素初始化為零，作為點擴散函數(shù)的初始值。之后計算圖像的中心位置。然后，計算運動方向的斜率的正切值以及余切值。根據(jù)運動方向的斜率的正切值以及余切值去判斷運動方向。如果斜率的正切值小于等于1，表示運動方向接近水平，所采取的操作是在點擴散函數(shù)的相應位置設置值為1，形成一條近似水平的模糊效果。這個相應位置的計算方式是中心加上偏移量。垂直方向同理。運動模糊這一塊的代碼具體如下。

def motion_process(image_size, motion_angle):# 創(chuàng)建一個大小為image_size的二維數(shù)組PSF，并將其所有元素初始化為零，作為點擴散函數(shù)的初始值。PSF = np.zeros(image_size)print(image_size)# 計算圖像的中心位置，由于數(shù)組索引是從0開始的，所以需要減1。center_position = (image_size[0] - 1) / 2print(center_position)# 計算運動方向的斜率的正切值。這里的motion_angle是以度為單位的運動方向角度。slope_tan = math.tan(motion_angle * math.pi/ 180)# 計算斜率的余切值。slope_cot = 1/slope_tan# 如果斜率的正切值小于等于1，表示運動方向接近水平。if slope_tan <=1:for i in range(15):# 計算相對于中心位置的水平偏移。round函數(shù)用于將浮點數(shù)四舍五入為整數(shù)。offset = round(i * slope_tan) #((center_position-i)*slope_tan)# 在點擴散函數(shù)的相應位置設置值為1，形成一條近似水平的模糊效果。PSF [int (center_position + offset), int(center_position - offset)] = 1# 將點擴散函數(shù)進行歸一化，確保其總和為1。return PSF/PSF.sum() #對點擴散函數(shù)進行歸一化亮度# 如果斜率的正切值大于1，表示運動方向接近垂直。else:for i in range(15):# 計算相對于中心位置的垂直偏移。offset = round(i * slope_cot)# ：在點擴散函數(shù)的相應位置設置值為1，形成一條近似垂直的模糊效果。PSF[int(center_position - offset),int(center_position + offset)] = 1#     將點擴散函數(shù)進行歸一化，確保其總和為1。return PSF / PSF.sum()

之后，對點擴散函數(shù)（PSF）進行二維傅里葉變換。這里加上e是為了避免除零錯誤。之后 將輸入圖像的傅里葉變換與點擴散函數(shù)的傅里葉變換進行逐元素相乘，在將結果進行逆傅里葉變換，得到模糊處理后的圖像。最后用fftshift將頻譜移到圖像中心，然后取絕對值，得到最終的模糊圖像。
總結一下，第一步是為了得到點擴散函數(shù)，可以理解為一個模糊的模板，這一步，是在進行時域卷積，頻域乘積，對輸入圖像和模糊模板進行卷積，得到輸出。
得到模糊圖像的代碼如下。

def make_blurred(input, PSF, eps):# 對輸入圖像進行二維傅里葉變換。input_fft = fft.fft2(input)  # 進行二維數(shù)組的傅里葉變換# 對點擴散函數(shù)（PSF）進行二維傅里葉變換，并添加一個小的常數(shù)eps以避免除零錯誤。這是在頻域進行卷積操作的準備步驟。PSF_fft = fft.fft2(PSF) + eps# input_fft * PSF_fft, 將輸入圖像的傅里葉變換與點擴散函數(shù)的傅里葉變換進行逐元素相乘，這相當于在時域中進行卷積。# 之后對相乘結果進行逆傅里葉變換，得到模糊處理后的圖像。blurred = fft.ifft2(input_fft * PSF_fft)# 用fftshift將頻譜移到圖像中心，然后取絕對值，得到最終的模糊圖像。blurred = np.abs(fft.fftshift(blurred))return blurred

2.逆濾波
對圖像和點擴散函數(shù)（PSF）分別進行傅里葉變換，同時將將PSF的頻域表示進行平移，將其中心移到圖像中心。這是為了避免逆濾波后的圖像進行旋轉。簡單而言，就是將頻域中的輸入圖像和PSF的傅里葉變換進行逐元素相除，然后進行逆傅里葉變換，就實現(xiàn)了逆濾波操作。具體代碼如下。

def inverse (input, PSF, eps):input_fft = fft.fft2(input) # 進行二維數(shù)組的傅里葉變換，將圖像轉換到頻域PSF_fft =fft.fft2(PSF)+ eps # 噪聲功率，這是已知的，考慮epsilon# 為避免逆濾波后的圖像進行旋轉，將PSF的中心移到圖像中心PSF_fft_shifted = fft.fftshift(PSF_fft)  # 將PSF的中心移到圖像中心# 對頻域中的輸入圖像和PSF的傅里葉變換進行逐元素相除，然后進行逆傅里葉變換。這一步實現(xiàn)了逆濾波操作。result = fft.fft2(input_fft / PSF_fft_shifted)# 對逆濾波得到的頻域結果進行頻譜中心平移，取其絕對值，得到逆濾波后的時域圖像。result = np.abs(fft.fftshift(result))return result

3.維納濾波
同理，對圖像和點擴散函數(shù)（PSF）分別進行傅里葉變換。之后，計算維納濾波的頻域濾波器。將輸入圖像的傅里葉變換與維納濾波器的傅里葉變換逐元素相乘，然后進行逆傅里葉變換，得到維納濾波后的時域結果。具體代碼如下。

def wiener(input,PSF,eps,K=0.01): # 維納濾波，K=0.01# 對輸入圖像進行二維傅里葉變換，將圖像轉換到頻域。input_fft = fft.fft2(input)# 點擴散函數(shù)（PSF）進行傅里葉變換，同時考慮了噪聲功率 eps。這一步是為了在頻域中進行維納濾波操作的準備。PSF_fft = fft.fft2(PSF) + eps# 計算維納濾波的頻域濾波器。np.conj()是復共軛操作，這里計算的是維納濾波器的分母，其中 K 是維納濾波的參數(shù)，用于控制噪聲增強的程度。PSF_fft_1 = np.conj(PSF_fft) /(np.abs(PSF_fft)** 2 + K)# 將輸入圖像的傅里葉變換與維納濾波器的傅里葉變換逐元素相乘，然后進行逆傅里葉變換，得到維納濾波后的時域結果。result = fft.ifft2(input_fft * PSF_fft_1)# 將維納濾波后的頻域結果進行頻譜中心平移，取其絕對值，得到維納濾波后的時域圖像。result = np.abs(fft.fftshift(result))return result

4.函數(shù)調用與繪圖
這塊比較簡單，直接附上代碼

image = cv2.imread('R.jpg')
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
img_h =image.shape[0]
img_w =image.shape[1]
plt.figure(1)
plt.rcParams['font.sans-serif'] = ['SimHei']  # 選擇一個包含中文字符的字體
plt.rcParams['axes.unicode_minus'] = False  # 用來正常顯示負號
plt.xlabel("原始圖像")
plt.gray()
plt.imshow(image) #顯示原圖像
plt.figure(2)
plt.gray()# 進行運動模糊處理
PSF = motion_process((img_h, img_w), 60)
blurred = np.abs(make_blurred(image, PSF, 1e-3))
plt.subplot(231)
plt.xlabel(" 進行模糊 ")
plt.imshow(blurred)
result = inverse(blurred, PSF, 1e-3)
# 逆濾波，對圖像進行濾波處理
result = np.flipud(result)  
result = np.fliplr(result)
plt.subplot(232)
plt.xlabel("對進行模糊的圖像進行逆濾波")
plt.imshow(result)result = wiener(blurred, PSF, 1e-3)  # 維納濾波
plt.subplot(233)
plt.xlabel("對模糊后的圖像進行維納濾波(k=0.01)")
plt.imshow(result)# 在模糊圖像上進一步添加呈正態(tài)分布的噪聲
blurred_noisy = blurred + 0.1 * blurred.std()*\np.random.standard_normal(blurred.shape)  # 添加噪聲,standard_normal# 產(chǎn)生隨機的函數(shù)
plt.subplot(234)
plt.xlabel("在模糊圖像上引入噪聲")
plt.imshow(blurred_noisy)  # 顯示添加噪聲且運動模糊的圖像result=inverse(blurred_noisyPSF,0.1+1e-3) # 對添加噪聲的圖像進行逆濾波
plt.subplot(235)
plt.xlabel("對模糊和加噪的圖片進行逆濾波")
plt.imshow(result)
result = wiener(blurred_noisy, PSF,0.1 + 1e-3)  # 對添加噪聲的圖像進行維納濾波
plt.subplot(236)
plt.xlabel("對模糊和加噪的圖片進行維納濾波(k=0.01) ")
plt.imshow(result)
plt.tight_layout()
plt.show()

5.運行結果


6.總結
本次實驗對原來的圖像進行運動模糊后，分別采用逆濾波和維納濾波進行對圖像的恢復。之后在模糊圖像的基礎上進一步添加呈標準正態(tài)分布的噪聲。再次，分別采用逆濾波和維納濾波進行對圖像的恢復。但是從實驗結果來看，還是存在一定的振鈴效應。
因為開頭要導入一些必要的庫。完整代碼如下。

import matplotlib.pyplot as plt
import numpy as np
from numpy import fft
import math
import cv2# 對退化過程進行建模
def motion_process(image_size, motion_angle):# 創(chuàng)建一個大小為image_size的二維數(shù)組PSF，并將其所有元素初始化為零，作為點擴散函數(shù)的初始值。PSF = np.zeros(image_size)print(image_size)# 計算圖像的中心位置，由于數(shù)組索引是從0開始的，所以需要減1。center_position = (image_size[0] - 1) / 2print(center_position)# 計算運動方向的斜率的正切值。這里的motion_angle是以度為單位的運動方向角度。slope_tan = math.tan(motion_angle * math.pi/ 180)# 計算斜率的余切值。slope_cot = 1/slope_tan# 如果斜率的正切值小于等于1，表示運動方向接近水平。if slope_tan <=1:for i in range(15):# 計算相對于中心位置的水平偏移。round函數(shù)用于將浮點數(shù)四舍五入為整數(shù)。offset = round(i * slope_tan) #((center_position-i)*slope_tan)# 在點擴散函數(shù)的相應位置設置值為1，形成一條近似水平的模糊效果。PSF [int (center_position + offset), int(center_position - offset)] = 1# 將點擴散函數(shù)進行歸一化，確保其總和為1。return PSF/PSF.sum() #對點擴散函數(shù)進行歸一化亮度# 如果斜率的正切值大于1，表示運動方向接近垂直。else:for i in range(15):# 計算相對于中心位置的垂直偏移。offset = round(i * slope_cot)# ：在點擴散函數(shù)的相應位置設置值為1，形成一條近似垂直的模糊效果。PSF[int(center_position - offset),int(center_position + offset)] = 1#     將點擴散函數(shù)進行歸一化，確保其總和為1。return PSF / PSF.sum()#
def make_blurred(input, PSF, eps):# 對輸入圖像進行二維傅里葉變換。input_fft = fft.fft2(input)  # 進行二維數(shù)組的傅里葉變換# 對點擴散函數(shù)（PSF）進行二維傅里葉變換，并添加一個小的常數(shù)eps以避免除零錯誤。這是在頻域進行卷積操作的準備步驟。PSF_fft = fft.fft2(PSF) + eps# input_fft * PSF_fft, 將輸入圖像的傅里葉變換與點擴散函數(shù)的傅里葉變換進行逐元素相乘，這相當于在時域中進行卷積。# 之后對相乘結果進行逆傅里葉變換，得到模糊處理后的圖像。blurred = fft.ifft2(input_fft * PSF_fft)# 用fftshift將頻譜移到圖像中心，然后取絕對值，得到最終的模糊圖像。blurred = np.abs(fft.fftshift(blurred))return blurred# 逆濾波的目標是盡可能地從經(jīng)過模糊和添加噪聲的圖像中恢復原始圖像
# 逆濾波
def inverse (input, PSF, eps):input_fft = fft.fft2(input) # 進行二維數(shù)組的傅里葉變換，將圖像轉換到頻域PSF_fft =fft.fft2(PSF)+ eps # 噪聲功率，這是已知的，考慮epsilon# 為避免逆濾波后的圖像進行旋轉，將PSF的中心移到圖像中心PSF_fft_shifted = fft.fftshift(PSF_fft)  # 將PSF的中心移到圖像中心# 對頻域中的輸入圖像和PSF的傅里葉變換進行逐元素相除，然后進行逆傅里葉變換。這一步實現(xiàn)了逆濾波操作。result = fft.fft2(input_fft / PSF_fft_shifted)# 對逆濾波得到的頻域結果進行頻譜中心平移，取其絕對值，得到逆濾波后的時域圖像。result = np.abs(fft.fftshift(result))return result
def wiener(input,PSF,eps,K=0.01): # 維納濾波，K=0.01# 對輸入圖像進行二維傅里葉變換，將圖像轉換到頻域。input_fft = fft.fft2(input)# 點擴散函數(shù)（PSF）進行傅里葉變換，同時考慮了噪聲功率 eps。這一步是為了在頻域中進行維納濾波操作的準備。PSF_fft = fft.fft2(PSF) + eps# 計算維納濾波的頻域濾波器。np.conj()是復共軛操作，這里計算的是維納濾波器的分母，其中 K 是維納濾波的參數(shù)，用于控制噪聲增強的程度。PSF_fft_1 = np.conj(PSF_fft) /(np.abs(PSF_fft)** 2 + K)# 將輸入圖像的傅里葉變換與維納濾波器的傅里葉變換逐元素相乘，然后進行逆傅里葉變換，得到維納濾波后的時域結果。result = fft.ifft2(input_fft * PSF_fft_1)# 將維納濾波后的頻域結果進行頻譜中心平移，取其絕對值，得到維納濾波后的時域圖像。result = np.abs(fft.fftshift(result))return resultimage = cv2.imread('R.jpg')
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
img_h =image.shape[0]
img_w =image.shape[1]
plt.figure(1)
plt.rcParams['font.sans-serif'] = ['SimHei']  # 選擇一個包含中文字符的字體
plt.rcParams['axes.unicode_minus'] = False  # 用來正常顯示負號
plt.xlabel("原始圖像")
plt.gray()
plt.imshow(image) #顯示原圖像
plt.figure(2)
plt.gray()# 進行運動模糊處理
PSF = motion_process((img_h, img_w), 60)
blurred = np.abs(make_blurred(image, PSF, 1e-3))
plt.subplot(231)
plt.xlabel(" 進行模糊 ")
plt.imshow(blurred)
result = inverse(blurred, PSF, 1e-3)
# 逆濾波，對圖像進行濾波處理
result = np.flipud(result)
result = np.fliplr(result)
plt.subplot(232)
plt.xlabel("對進行模糊的圖像進行逆濾波")
plt.imshow(result)result = wiener(blurred, PSF, 1e-3)  # 維納濾波
plt.subplot(233)
plt.xlabel("對模糊后的圖像進行維納濾波(k=0.01)")
plt.imshow(result)# 在模糊圖像上進一步添加呈正態(tài)分布的噪聲
blurred_noisy = blurred + 0.1 * blurred.std()*\np.random.standard_normal(blurred.shape)  # 添加噪聲,standard_normal# 產(chǎn)生隨機的函數(shù)
plt.subplot(234)
plt.xlabel("在模糊圖像上引入噪聲")
plt.imshow(blurred_noisy)  # 顯示添加噪聲且運動模糊的圖像result=inverse(blurred_noisyPSF,0.1+1e-3) # 對添加噪聲的圖像進行逆濾波
plt.subplot(235)
plt.xlabel("對模糊和加噪的圖片進行逆濾波")
plt.imshow(result)
result = wiener(blurred_noisy, PSF,0.1 + 1e-3)  # 對添加噪聲的圖像進行維納濾波
plt.subplot(236)
plt.xlabel("對模糊和加噪的圖片進行維納濾波(k=0.01) ")
plt.imshow(result)
plt.tight_layout()
plt.show()