OpenCV—Python 图像去模糊（维纳滤波，约束最⼩⼆乘⽅滤波）
⽂章⽬录
⼀、维纳滤波
对于运动引起的图像模糊，最简单的⽅法是直接做逆滤波，但是逆滤波对加性噪声特别敏感，使得恢复的图像⼏乎不可⽤。最⼩均⽅差（维纳）滤波⽤来去除含有噪声的模糊图像，其⽬标是找到未污染图像的⼀个估计，使它们之间的均⽅差最⼩，可以去除噪声，同时清晰化模糊图像。
其中：
是卷积符号
是在时间 刻输⼊的信号（未知）
是⼀个线性时间不变系统的脉冲响应（已知）
是加性噪声，与 不相关（未知）
是我们观察到的信号
我们的⽬标是找出这样的卷积函数 ，这样我们可以如下得到估计的
这⾥是的最⼩均⽅差估计。
基于这种误差度量， 滤波器可以在频率域如下描述
这⾥：
)和 是 和 在频率域ff的傅⾥叶变换。
是输⼊信号 的功率谱。
是噪声的 的功率谱。
差旅费报销单上标 代表复数共轭。
滤波过程可以在频率域完成： 这⾥ 是 的傅⾥叶变换，通过逆傅⾥叶变化可以得到去卷积后的结果。
上⾯的式⼦可以改写成更为清晰的形式：
这⾥ 是 在频率域 的傅⾥叶变换。是信号噪声⽐。当噪声为零时（即信噪⽐趋近于⽆穷），⽅括号
内各项也就等于1，意味着此时刻维纳滤波也就简化成逆滤波过程。但是当噪声增加时，信噪⽐降低，⽅括号⾥⾯值也跟着降低。这说明，维纳滤波的带通频率依赖于信噪⽐。
代码⽰例：如下代码参考于，⾃⼰解决了对RGB图⽚的⽀持，与诸位共勉。
import matplotlib .pyplot as plt
import numpy as np
from numpy import fft
programmey (t )=h (t )x (t )+⨂n (t )
⨂x (t )t h (t )n (t )x (t )y (t )g (t )x (t )(t )=x
^g (t )∗y (t )(t )x ^x (t )G (f )==∣H (f )∣S (f )+N (f )2H (f )S (f )∗∣H (f )∣+N (f )/S (f )
2H (f )
∗G (f )H (f )g h S (f )x (t )N (f )n (t )∗(f )=X
^G (f )∗Y (f )(f )X ^(t )x ^(t )x ^G (f )==H (f )1⎣⎡∣H (f )∣+2S (f )N (f )∣H (f )∣2⎦⎤H (f )1[∣H (f )∣+2SNR (f )1∣H (f )∣2]
H (f )h f SNR (f )=S (f )/N (f )
from numpy import fft
import math
import cv2
# 仿真运动模糊
def motion_process(image_size, motion_angle):
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_tan
if slope_tan <=1:
for i in range(15):
offt =round(i * slope_tan)# ((center_position-i)*slope_tan)
PSF[int(center_position + offt),int(center_position - offt)]=1
return PSF / PSF.sum()# 对点扩散函数进⾏归⼀化亮度
tci
el:
for i in range(15):
offt =round(i * slope_cot)
PSF[int(center_position - offt),int(center_position + offt)]=1
return PSF / PSF.sum()
# 对图⽚进⾏运动模糊
def make_blurred(input, PSF, eps):
input_fft = fft.fft2(input)# 进⾏⼆维数组的傅⾥叶变换
夏威夷误报导弹PSF_fft = fft.fft2(PSF)+ eps
blurred = fft.ifft2(input_fft * PSF_fft)
blurred = np.abs(fft.fftshift(blurred))
return blurred
def inver(input, PSF, eps):# 逆滤波
input_fft = fft.fft2(input)
PSF_fft = fft.fft2(PSF)+ eps # 噪声功率，这是已知的，考虑epsilon
result = fft.ifft2(input_fft / PSF_fft)# 计算F(u,v)的傅⾥叶反变换
result = np.abs(fft.fftshift(result))
return result
airport是什么意思def wiener(input, PSF, eps, K=0.01):# 维纳滤波，K=0.01
input_fft = fft.fft2(input)
PSF_fft = fft.fft2(PSF)+ eps
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
def normal(array):
array = np.where(array <0,0, array)
array = np.where(array >255,255, array)
array = array.astype(np.int16)
return array
def main(gray):
channel =[]
img_h, img_w = gray.shape[:2]
PSF = motion_process((img_h, img_w),60)# 进⾏运动模糊处理
blurred = np.abs(make_blurred(gray, PSF,1e-3))
result_blurred = inver(blurred, PSF,1e-3)# 逆滤波
result_wiener = wiener(blurred, PSF,1e-3)# 维纳滤波
blurred_noisy = blurred +0.1* blurred.std()* \
np.random.standard_normal(blurred.shape)# 添加噪声,standard_normal产⽣随机的函数 inver_mo2no = inver(blurred_noisy, PSF,0.1+1e-3)# 对添加噪声的图像进⾏逆滤波
inver_mo2no = inver(blurred_noisy, PSF,0.1+1e-3)# 对添加噪声的图像进⾏逆滤波 wiener_mo2no = wiener(blurred_noisy, PSF,0.1+1e-3)# 对添加噪声的图像进⾏维纳滤波 channel.append((normal(blurred),normal(result_blurred),normal(result_wiener),
normal(blurred_noisy),normal(inver_mo2no),normal(wiener_mo2no)))
return channel
if __name__ =='__main__':
image = cv2.imread('./gggg/001.png')全国大学生英语四六级官网
b_gray, g_gray, r_gray = cv2.py())
Result =[]
for gray in[b_gray, g_gray, r_gray]:
channel = main(gray)
Result.append(channel)
blurred = ([Result[0][0][0], Result[1][0][0], Result[2][0][0]])
result_blurred = ([Result[0][0][1], Result[1][0][1], Result[2][0][1]])
result_wiener = ([Result[0][0][2], Result[1][0][2], Result[2][0][2]])
blurred_noisy = ([Result[0][0][3], Result[1][0][3], Result[2][0][3]])
inver_mo2no = ([Result[0][0][4], Result[1][0][4], Result[2][0][4]])
mulberrysilk
wiener_mo2no = ([Result[0][0][5], Result[1][0][5], Result[2][0][5]])
#========= 可视化 ==========
plt.figure(1)
plt.xlabel("Original Image")
watertownplt.imshow(np.flip(image, axis=2))# 显⽰原图像
plt.figure(2)
plt.figure(figsize=(8,6.5))
imgNames ={"Motion blurred":blurred,
"inver deblurred":result_blurred,
"wiener deblurred(k=0.01)":result_wiener,
"motion & noisy blurred":blurred_noisy,
"inver_mo2no":inver_mo2no,
'wiener_mo2no':wiener_mo2no}
for i,(key,imgName)in enumerate(imgNames.items()):
plt.subplot(231+i)
plt.xlabel(key)
plt.imshow(np.flip(imgName, axis=2))
plt.show()
⼆、约束最⼩⼆乘⽅滤波
约束最⼩⼆乘⽅滤波（Constrained Least Squares Filtering，aka Tikhonov filtration,Tikhonov regularization）核⼼是H对噪声的敏感性问题。减少噪声敏感新问题的⼀种⽅法是以平滑度量的最佳复原为基础的，建⽴下列约束条件：
约束条件：这⾥，是为退化图像的估计，为加性噪声，拉普拉斯算⼦在这⾥表⽰平滑程度。推导：
将上式表⽰成矩阵形式，同时将约束项转换成拉格朗⽇乘⼦项：
最⼩化上代价函数，对求导，令等零有：最后可得到：
是函数
C =[∇f (x ,y )]0∑M −10∑N −1
22
∣∣G −H ∣∣=F ^22∣∣N ∣∣22
F
^N ∇2∣∣P ∣∣−F ^22λ∣∣G −H ∣∣−∣∣N ∣∣(F ^2222
)F
^P P =∗F ^λH (G −H )∗F ^=F
^=λH H +P P ∗∗λH G ∗G [∣∣H ∣∣+λ∣∣P ∣∣2222H ∗]P P =⎣⎡0−10−14−10−10
⎦⎤
三、psf2otf ,circShift
circshift(psf,K)
当K为⼀个数字时，只对矩阵进⾏上下平移，当K为⼀个坐标时，会对矩阵进⾏上下和左右两个⽅向进⾏平移。⽰例如下：执⾏：平移坐标(-1,-1)，对矩阵进⾏上移，左移1个单位，效果如下:
⽰例2：平移坐标(1,2)，对矩阵进⾏下移1个单位，右移2个单位，效果如下:
psf2otf()
# coding: utf-8
import numpy as np
import matplotlib .pyplot as plt
from numpy import fft有朋自远方来英语
import cv2
from temp_004 import psf2otf
def motion_blur (gray , degree =7, angle =60):
gray = np .array (gray )
M = RotationMatrix2D ((round (degree / 2), round (degree / 2)), angle , 1)
motion_blur_kernel = np .diag (np .ones (degree ))
motion_blur_kernel = cv2.warpAffine (motion_blur_kernel , M , (degree , degree ))
PSF = motion_blur_kernel / degree
blurred = cv2.filter2D (gray , -1, PSF )
blurred = alize (blurred ,None , 0, 255, cv2.NORM_MINMAX )
游戏培训机构
blurred = np .array (blurred , dtype =np .uint8)
return blurred ,PSF
def inver (blurred , PF ):
IF_fft = fft .fft2(blurred )
result = fft .ifft2(IF_fft / PF )
result = np .real (result )
return result
def wiener (blurred , PF , SNR =0.01): # 维纳滤波，K=0.01
IF_fft = fft .fft2(blurred )
G_f = np .conj (PF ) / (np .abs (PF ) ** 2 + SNR )
result = fft .ifft2(IF_fft * G_f )
result = np .real (result )
return result
def CLSF (blurred ,PF ,gamma = 0.05):
outheight , outwidth = blurred .shape [:2]
kernel = np .array ([[0, -1, 0],
[-1, 4, -1],
[0, -1, 0]])⇒
⎣⎡0−10−14−10−10⎦⎤⎣⎡4−1−1−100−100⎦⎤⇒⎣⎢⎢⎡1611162712173813184914195101520⎦⎥⎥⎤⎣⎢⎢⎡1949142051015161611172712183813
⎦⎥⎥⎤
