DIP——添加运动模糊与滤波

1.运动模糊
为了模拟图像退化的过程，在这里创建了一个用于模拟运动模糊的点扩散函数，具体模糊的方向取决于输入的motion_angle。如果运动方向接近水平，则模糊效果近似水平，如果运动方向接近垂直，则模糊效果近似垂直。具体操作如下。首先创建二维数组PSF，并将其所有元素初始化为零，作为点扩散函数的初始值。之后计算图像的中心位置。然后，计算运动方向的斜率的正切值以及余切值。根据运动方向的斜率的正切值以及余切值去判断运动方向。如果斜率的正切值小于等于1，表示运动方向接近水平，所采取的操作是在点扩散函数的相应位置设置值为1，形成一条近似水平的模糊效果。这个相应位置的计算方式是中心加上偏移量。垂直方向同理。运动模糊这一块的代码具体如下。

def motion_process(image_size, motion_angle):# 创建一个大小为image_size的二维数组PSF，并将其所有元素初始化为零，作为点扩散函数的初始值。PSF = np.zeros(image_size)print(image_size)# 计算图像的中心位置，由于数组索引是从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函数用于将浮点数四舍五入为整数。offset = round(i * slope_tan) #((center_position-i)*slope_tan)# 在点扩散函数的相应位置设置值为1，形成一条近似水平的模糊效果。PSF [int (center_position + offset), int(center_position - offset)] = 1# 将点扩散函数进行归一化，确保其总和为1。return PSF/PSF.sum() #对点扩散函数进行归一化亮度# 如果斜率的正切值大于1，表示运动方向接近垂直。else:for i in range(15):# 计算相对于中心位置的垂直偏移。offset = round(i * slope_cot)# ：在点扩散函数的相应位置设置值为1，形成一条近似垂直的模糊效果。PSF[int(center_position - offset),int(center_position + offset)] = 1#     将点扩散函数进行归一化，确保其总和为1。return PSF / PSF.sum()

之后，对点扩散函数（PSF）进行二维傅里叶变换。这里加上e是为了避免除零错误。之后 将输入图像的傅里叶变换与点扩散函数的傅里叶变换进行逐元素相乘，在将结果进行逆傅里叶变换，得到模糊处理后的图像。最后用fftshift将频谱移到图像中心，然后取绝对值，得到最终的模糊图像。
总结一下，第一步是为了得到点扩散函数，可以理解为一个模糊的模板，这一步，是在进行时域卷积，频域乘积，对输入图像和模糊模板进行卷积，得到输出。
得到模糊图像的代码如下。

def make_blurred(input, PSF, eps):# 对输入图像进行二维傅里叶变换。input_fft = fft.fft2(input)  # 进行二维数组的傅里叶变换# 对点扩散函数（PSF）进行二维傅里叶变换，并添加一个小的常数eps以避免除零错误。这是在频域进行卷积操作的准备步骤。PSF_fft = fft.fft2(PSF) + eps# input_fft * PSF_fft, 将输入图像的傅里叶变换与点扩散函数的傅里叶变换进行逐元素相乘，这相当于在时域中进行卷积。# 之后对相乘结果进行逆傅里叶变换，得到模糊处理后的图像。blurred = fft.ifft2(input_fft * PSF_fft)# 用fftshift将频谱移到图像中心，然后取绝对值，得到最终的模糊图像。blurred = np.abs(fft.fftshift(blurred))return blurred

2.逆滤波
对图像和点扩散函数（PSF）分别进行傅里叶变换，同时将将PSF的频域表示进行平移，将其中心移到图像中心。这是为了避免逆滤波后的图像进行旋转。简单而言，就是将频域中的输入图像和PSF的傅里叶变换进行逐元素相除，然后进行逆傅里叶变换，就实现了逆滤波操作。具体代码如下。

def inverse (input, PSF, eps):input_fft = fft.fft2(input) # 进行二维数组的傅里叶变换，将图像转换到频域PSF_fft =fft.fft2(PSF)+ eps # 噪声功率，这是已知的，考虑epsilon# 为避免逆滤波后的图像进行旋转，将PSF的中心移到图像中心PSF_fft_shifted = fft.fftshift(PSF_fft)  # 将PSF的中心移到图像中心# 对频域中的输入图像和PSF的傅里叶变换进行逐元素相除，然后进行逆傅里叶变换。这一步实现了逆滤波操作。result = fft.fft2(input_fft / PSF_fft_shifted)# 对逆滤波得到的频域结果进行频谱中心平移，取其绝对值，得到逆滤波后的时域图像。result = np.abs(fft.fftshift(result))return result

3.维纳滤波
同理，对图像和点扩散函数（PSF）分别进行傅里叶变换。之后，计算维纳滤波的频域滤波器。将输入图像的傅里叶变换与维纳滤波器的傅里叶变换逐元素相乘，然后进行逆傅里叶变换，得到维纳滤波后的时域结果。具体代码如下。

def wiener(input,PSF,eps,K=0.01): # 维纳滤波，K=0.01# 对输入图像进行二维傅里叶变换，将图像转换到频域。input_fft = fft.fft2(input)# 点扩散函数（PSF）进行傅里叶变换，同时考虑了噪声功率 eps。这一步是为了在频域中进行维纳滤波操作的准备。PSF_fft = fft.fft2(PSF) + eps# 计算维纳滤波的频域滤波器。np.conj()是复共轭操作，这里计算的是维纳滤波器的分母，其中 K 是维纳滤波的参数，用于控制噪声增强的程度。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.函数调用与绘图
这块比较简单，直接附上代码

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)# 在模糊图像上进一步添加呈正态分布的噪声
blurred_noisy = blurred + 0.1 * blurred.std()*\np.random.standard_normal(blurred.shape)  # 添加噪声,standard_normal# 产生随机的函数
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.总结
本次实验对原来的图像进行运动模糊后，分别采用逆滤波和维纳滤波进行对图像的恢复。之后在模糊图像的基础上进一步添加呈标准正态分布的噪声。再次，分别采用逆滤波和维纳滤波进行对图像的恢复。但是从实验结果来看，还是存在一定的振铃效应。
因为开头要导入一些必要的库。完整代码如下。

import matplotlib.pyplot as plt
import numpy as np
from numpy import fft
import math
import cv2# 对退化过程进行建模
def motion_process(image_size, motion_angle):# 创建一个大小为image_size的二维数组PSF，并将其所有元素初始化为零，作为点扩散函数的初始值。PSF = np.zeros(image_size)print(image_size)# 计算图像的中心位置，由于数组索引是从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函数用于将浮点数四舍五入为整数。offset = round(i * slope_tan) #((center_position-i)*slope_tan)# 在点扩散函数的相应位置设置值为1，形成一条近似水平的模糊效果。PSF [int (center_position + offset), int(center_position - offset)] = 1# 将点扩散函数进行归一化，确保其总和为1。return PSF/PSF.sum() #对点扩散函数进行归一化亮度# 如果斜率的正切值大于1，表示运动方向接近垂直。else:for i in range(15):# 计算相对于中心位置的垂直偏移。offset = round(i * slope_cot)# ：在点扩散函数的相应位置设置值为1，形成一条近似垂直的模糊效果。PSF[int(center_position - offset),int(center_position + offset)] = 1#     将点扩散函数进行归一化，确保其总和为1。return PSF / PSF.sum()#
def make_blurred(input, PSF, eps):# 对输入图像进行二维傅里叶变换。input_fft = fft.fft2(input)  # 进行二维数组的傅里叶变换# 对点扩散函数（PSF）进行二维傅里叶变换，并添加一个小的常数eps以避免除零错误。这是在频域进行卷积操作的准备步骤。PSF_fft = fft.fft2(PSF) + eps# input_fft * PSF_fft, 将输入图像的傅里叶变换与点扩散函数的傅里叶变换进行逐元素相乘，这相当于在时域中进行卷积。# 之后对相乘结果进行逆傅里叶变换，得到模糊处理后的图像。blurred = fft.ifft2(input_fft * PSF_fft)# 用fftshift将频谱移到图像中心，然后取绝对值，得到最终的模糊图像。blurred = np.abs(fft.fftshift(blurred))return blurred# 逆滤波的目标是尽可能地从经过模糊和添加噪声的图像中恢复原始图像
# 逆滤波
def inverse (input, PSF, eps):input_fft = fft.fft2(input) # 进行二维数组的傅里叶变换，将图像转换到频域PSF_fft =fft.fft2(PSF)+ eps # 噪声功率，这是已知的，考虑epsilon# 为避免逆滤波后的图像进行旋转，将PSF的中心移到图像中心PSF_fft_shifted = fft.fftshift(PSF_fft)  # 将PSF的中心移到图像中心# 对频域中的输入图像和PSF的傅里叶变换进行逐元素相除，然后进行逆傅里叶变换。这一步实现了逆滤波操作。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)# 点扩散函数（PSF）进行傅里叶变换，同时考虑了噪声功率 eps。这一步是为了在频域中进行维纳滤波操作的准备。PSF_fft = fft.fft2(PSF) + eps# 计算维纳滤波的频域滤波器。np.conj()是复共轭操作，这里计算的是维纳滤波器的分母，其中 K 是维纳滤波的参数，用于控制噪声增强的程度。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)# 在模糊图像上进一步添加呈正态分布的噪声
blurred_noisy = blurred + 0.1 * blurred.std()*\np.random.standard_normal(blurred.shape)  # 添加噪声,standard_normal# 产生随机的函数
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()