前言
本文主要介绍退化图像复原的两种方法：逆滤波和维纳滤波。
一、逆滤波
图像退化的表达式：
g
(
x
,
y
)
=
h
(
x
,
y
)
⊙
f
(
x
,
y
)
+
η
(
x
,
y
)
\begin{aligned} g(x,y)=h(x,y)\odot f(x,y)+\eta(x,y) \end{aligned}
g(x,y)=h(x,y)⊙f(x,y)+η(x,y)
f
(
x
,
y
)
:
f(x,y):
f(x,y):输入图像
h
(
x
,
y
)
:
h(x,y):
h(x,y):退化函数
η
(
x
,
y
)
:
\eta(x,y):
η(x,y):噪声项
g
:
g:
g: 退化图像
由卷积定理转换到频域中：
G
(
u
,
v
)
=
H
(
u
,
v
)
∗
F
(
u
,
v
)
+
N
(
u
,
v
)
\begin{aligned} G(u,v)=H(u,v)*F(u,v)+N(u,v) \end{aligned}
G(u,v)=H(u,v)∗F(u,v)+N(u,v)
逆滤波运算就是将退化函数去除：
F
^
=
G
(
u
,
v
)
H
(
u
,
v
)
=
F
(
u
,
v
)
+
N
(
u
,
v
)
H
(
u
,
v
)
\begin{aligned} \hat{F}=\displaystyle\frac{G(u,v)}{H(u,v)}=F(u,v)+\frac{N(u,v)}{H(u,v)} \end{aligned}
F^=H(u,v)G(u,v)=F(u,v)+H(u,v)N(u,v)
根据上式，逆滤波存在两个关键问题：
- 退化函数 H ( u , v ) H(u,v) H(u,v)的估计的确准确，复原的越清晰。
- 当退化函数 H ( u , v ) H(u,v) H(u,v)趋向于0时，后面一项的值变得过大，最终的式子受噪声项 N ( u , v ) N(u,v) N(u,v)的影响就加大。
1.1 估计退化函数 H ( u , v ) H(u,v) H(u,v)
主要有3种方法：
(
1
)
(1)
(1):观察法；
(
2
)
(2)
(2):试验法；
(
3
)
(3)
(3):数学建模法
这些方法都是在频域上进行处理
1.1.1 观察法
寻找一个信号内容很强的区域（高对比度区域），从而忽略噪声项的影响，即：
H
s
(
u
,
v
)
=
G
s
(
u
,
v
)
F
s
^
(
u
,
v
)
\begin{aligned} {H_{s}(u,v)}=\displaystyle\frac{G_{s}(u,v)}{\hat{F_{s}}(u,v)} \end{aligned}
Hs(u,v)=Fs^(u,v)Gs(u,v)
G
s
(
u
,
v
)
:
{G_{s}}(u,v):
Gs(u,v):当前需要复原的图像的傅里叶变换。
F
s
^
(
u
,
v
)
:
{\hat{F_{s}}(u,v)}:
Fs^(u,v):已经经过处理（锐化、均值滤波等）复原后的图像。
1.1.2 试验法
简而言之就是做实验模拟真实退化的过程，再根据下面的式子得到退化函数。
H
(
u
,
v
)
=
G
(
u
,
v
)
A
\begin{aligned} {H(u,v)}=\displaystyle\frac{G(u,v)}{A} \end{aligned}
H(u,v)=AG(u,v)
假设现在有某个设备可以得到和退化图像非常近似的图像，然后对一个冲激（小光点，小光点应亮到能降低噪声的影响）成像，得到冲激退化后的FFT（对应上式的
G
(
u
,
v
)
G(u,v)
G(u,v)），而一个冲激的FFT是一个常量
A
A
A。
1.1.3 建模法 ★ \bigstar ★
(
1
)
(1)
(1)大气湍流模型
H
(
u
,
v
)
=
e
−
k
(
u
2
+
v
2
)
5
/
6
\begin{aligned} {H(u,v)}=e^{-k(u^{2}+v^{2})^{5/6}} \end{aligned}
H(u,v)=e−k(u2+v2)5/6
k 是湍流常数，反映湍流的强烈程度。
k
=
{
0.00025
轻微湍流
0.001
中等湍流
0.0025
剧烈湍流
k= \left\{ \begin{array}{ll} 0.00025\quad\quad轻微湍流\\ 0.001\quad\quad\quad中等湍流\\0.0025\quad\quad\quad剧烈湍流 \end{array} \right.
k=⎩
⎨
⎧0.00025轻微湍流0.001中等湍流0.0025剧烈湍流
代码：
import cv2
import numpy as np
from matplotlib import pyplot as plt
# 读取图像
img = cv2.imread('A.png',0)
# 进行傅里叶变换
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
#-------------大气湍流模型-----------------
k = 0.001 #湍流程度
degeneration_img = np.zeros(fshift.shape,dtype=complex) #一定要指定类型
for u in range(fshift.shape[0]):
for v in range(fshift.shape[1]):
H = np.exp(-k*np.power(np.float64((u-fshift.shape[0]//2)**2+(v-fshift.shape[1]//2)**2),5/6))
degeneration_img[u][v] = fshift[u][v] * H
#-----------------------------------------
# 将频域图像在空域进行显示做的处理
img1 = np.fft.ifftshift(degeneration_img)
img1 = np.fft.ifft2(img1)
img1 = np.abs(img1)
img1 = cv2.normalize(img1, None, 0, 255, cv2.NORM_MINMAX)
# 显示原始图像和退化的图像
plt.figure(figsize=(20, 20))
plt.subplot(121), plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(img1, cmap='gray')
plt.title('degeneration img'), plt.xticks([]), plt.yticks([])
plt.show()
注意：在对degeneration_img进行初始化时，一定要指定类型dtype=complex。我在写的过程中，就忘记对类型进行指定，导致得到了错误的退化图像（大家可以尝试一下）。一般来讲不指定类型并没有多大问题，但是，由于这里是复数类型，如果不指定类型，后面语句 degeneration_img[u][v] = fshift[u][v] * H 赋值的过程中就会导致虚部自动被舍弃掉了！参考下图：
以前不清楚动态语言的坑，现在算是掉过一次坑了，以后Python写代码一定一定加类型标注！！
另外，也可以用矩阵的方法生成：（下面代码参考文章：https://blog.csdn.net/youcans/article/details/123027287）
def turbulenceBlur(img, k=0.001): # 湍流模糊传递函数
# H(u,v) = exp(-k(u^2+v^2)^5/6)
M, N = img.shape[1], img.shape[0]
u, v = np.meshgrid(np.arange(M), np.arange(N))
radius = (u - M//2)**2 + (v - N//2)**2
kernel = np.exp(-k * np.power(radius, 5/6))
return kernel
(
2
)
(2)
(2)运动模糊模型
当图像在
x
x
x 方向做速率为
x
0
(
t
)
=
a
t
/
T
x_{0}(t) = at/T
x0(t)=at/T 的匀速直线运动 ,同时在
y
y
y 方向上做速率为
y
0
(
t
)
=
b
t
/
T
y_{0}(t) = bt/T
y0(t)=bt/T 的匀速直线运动。
H
(
u
,
v
)
=
T
π
(
u
a
+
v
b
)
s
i
n
[
π
(
u
a
+
v
b
)
]
e
−
j
π
(
u
a
+
v
b
)
\begin{aligned} H(u,v)=\frac{T}{\pi(ua+vb)}sin[\pi(ua+vb)]e^{-j\pi(ua+vb)} \end{aligned}
H(u,v)=π(ua+vb)Tsin[π(ua+vb)]e−jπ(ua+vb)
代码：
# 功能：实现运动模糊的效果
import cv2
import numpy as np
from matplotlib import pyplot as plt
# 读取图像
img = cv2.imread('0526.tif',0)
# 进行傅里叶变换
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
#-------------运动模糊模型-----------------
a = b = 0.1
T = 1
degeneration_img = np.zeros(fshift.shape,dtype=np.complex64) #一定要指定类型
for u in range(fshift.shape[0]):
for v in range(fshift.shape[1]):
uv = ((u-fshift.shape[0]//2)*a + (v-fshift.shape[1]//2)*b) * np.pi
if uv == 0:
degeneration_img[u][v] = fshift[u][v]
else:
H = T * np.sin(uv) * np.exp(np.complex64(-1j)*uv) /uv
degeneration_img[u][v] = fshift[u][v] * H
#-----------------------------------------
# 将频域图像在空域进行显示做的处理
img1 = np.fft.ifftshift(degeneration_img)
img1 = np.fft.ifft2(img1)
img1 = np.abs(img1)
img1 = cv2.normalize(img1, None, 0, 1, cv2.NORM_MINMAX)
# 显示原始图像和退化的图像
plt.figure(figsize=(20, 20))
plt.subplot(121), plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(img1, cmap='gray')
plt.title('degeneration img'), plt.xticks([]), plt.yticks([])
plt.show()
(
3
)
(3)
(3)高斯模糊模型
H
(
u
,
v
)
=
e
−
D
2
(
u
,
v
)
/
2
D
0
2
D
(
u
,
v
)
=
[
(
u
−
P
/
2
)
2
+
(
v
−
Q
/
2
)
2
]
1
/
2
\begin{aligned} {H(u,v)}=e^{-D^{2}(u,v)/2D_{0}^{2}}\quad\quad\\{D(u,v) = {[(u-P/2)^{2}+(v-Q/2)^{2}]^{1/2}}} \end{aligned}
H(u,v)=e−D2(u,v)/2D02D(u,v)=[(u−P/2)2+(v−Q/2)2]1/2
代码：
import cv2
import numpy as np
# 读取输入图像
img = cv2.imread('A.png', 0)
# 进行傅里叶变换
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
# 构造高斯滤波器
rows, cols = img.shape
crow, ccol = rows//2, cols//2
d = 60 # 滤波器大小
gauss = np.zeros((rows, cols))
for i in range(rows):
for j in range(cols):
gauss[i, j] = np.exp(-((i-crow)**2+(j-ccol)**2)/(2*d**2))
# 将高斯滤波器应用于频域图像
filtered = np.multiply(fshift, gauss)
# 进行反傅里叶变换
f_ishift = np.fft.ifftshift(filtered)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)
img_back = cv2.normalize(img_back, None, 0, 1, cv2.NORM_MINMAX)
# 显示结果
cv2.imshow('Input', img)
cv2.imshow('Output', img_back)
cv2.waitKey(0)
cv2.destroyAllWindows()
1.2 直接逆滤波
不考虑噪声：
F
^
=
G
(
u
,
v
)
H
(
u
,
v
)
\begin{aligned} \hat{F}=\displaystyle\frac{G(u,v)}{H(u,v)} \end{aligned}
F^=H(u,v)G(u,v)
【noise1.png】:经过高斯模糊并加上高斯噪声的图片。
在没有任何先验知识的情况下，使用大气湍流模型估计图片的退化函数，然后进行逆滤波处理。
代码如下：
import cv2
import numpy as np
from matplotlib import pyplot as plt
# 读取图像
img = cv2.imread('noise1.png',0)
# 进行傅里叶变换
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
#-------------利用大气湍流模型估计退化传递函数-----------------
k = 0.001 #湍流程度
degeneration = np.zeros(fshift.shape,dtype=complex)
for u in range(fshift.shape[0]):
for v in range(fshift.shape[1]):
H = np.exp(-k*np.power(np.float64((u+fshift.shape[0]//2)**2+(v-fshift.shape[1]//2)**2),5/6))
degeneration[u][v] = H
#-----------------------------------------
# 逆滤波
img1 = np.divide(fshift,degeneration)
# 将频域图像在空域进行显示做的处理
img1 = np.fft.ifftshift(img1)
img1 = np.fft.ifft2(img1)
img1 = np.abs(img1)
img1 = cv2.normalize(img1, None, 0, 255, cv2.NORM_MINMAX)
# 显示原始图像和恢复后的图像
plt.figure(figsize=(20, 20))
plt.subplot(121), plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(img1, cmap='gray')
plt.title('rec img'), plt.xticks([]), plt.yticks([])
plt.show()
滤波效果：
但是，只需要将上面语句degeneration[u][v] = H
→
\rightarrow
→degeneration[u][v] = H + 0.001，将0.001作为噪声的估计，就得到一定的效果：
1.3 半径受限逆滤波
考虑噪声：
F
^
=
F
(
u
,
v
)
+
N
(
u
,
v
)
H
(
u
,
v
)
\begin{aligned} \hat{F}=F(u,v)+\frac{N(u,v)}{H(u,v)} \end{aligned}
F^=F(u,v)+H(u,v)N(u,v)
使用低通滤波器过滤掉高频的噪声之后再除以H(u, v)
代码：
import cv2
import numpy as np
from matplotlib import pyplot as plt
# 读取图像
img = cv2.imread('noise3.png',0)
# 进行傅里叶变换
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
#-------------大气湍流模型-----------------
k = 0.001 #湍流程度
degeneration = np.zeros(fshift.shape,dtype=complex) #要指定类型
for u in range(fshift.shape[0]):
for v in range(fshift.shape[1]):
H = np.exp(-k*np.power(np.float64((u+fshift.shape[0]//2)**2+(v-fshift.shape[1]//2)**2),5/6))
degeneration[u][v] = H
#-----------------------------------------
# 构建低通滤波器---理想低通
# rows,cols = fshift.shape
# radius = 1
# mask = np.zeros((rows, cols),dtype=np.float64)
# for i in range(rows):
# for j in range(cols):
# if np.sqrt((i-rows//2)**2 + (j-cols//2)**2) <= radius:
# mask[i, j] = 0.1
# fshift = np.multiply(fshift, mask)
# # 构建高斯低通滤波器
# rows, cols = img.shape
# crow, ccol = rows//2, cols//2
# d = 50 # 滤波器大小
# gauss = np.zeros((rows, cols))
# for i in range(rows):
# for j in range(cols):
# gauss[i, j] = np.exp(-((i-crow)**2+(j-ccol)**2)/(2*d**2))
# 构建巴特沃斯低通滤波器
rows, cols = img.shape
crow, ccol = rows//2, cols//2
d = 70 # 滤波器大小
butter = np.zeros((rows, cols))
for i in range(rows):
for j in range(cols):
butter[i, j] = 1/(1 + np.power((np.power((i-crow)**2+(j-ccol)**2,0.5)/d),20))
fshift = np.multiply(fshift, butter)
# 逆滤波
img1 = np.divide(fshift,butter)
# 将频域图像在空域进行显示做的处理
img1 = np.fft.ifftshift(img1)
img1 = np.fft.ifft2(img1)
img1 = np.abs(img1)
img1 = cv2.normalize(img1, None, 0, 1, cv2.NORM_MINMAX)
# 显示原始图像和恢复后的图像
plt.figure(figsize=(20, 20))
plt.subplot(121), plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(img1, cmap='gray')
plt.title('rec img'), plt.xticks([]), plt.yticks([])
plt.show()
cv2.imshow('org',img)
cv2.imshow('rev',img1)
cv2.waitKey(0)
滤波效果：
二、最小均方误差（维纳）滤波
退化图像的估计的傅里叶变换为：
F
(
u
,
v
)
^
=
[
1
H
(
u
,
v
)
∣
H
(
u
,
v
)
∣
2
∣
H
(
u
,
v
)
∣
2
+
K
]
G
(
u
,
v
)
\begin{aligned} \hat{F(u,v)}=[\frac{1}{H(u,v)}\frac{|H(u,v)|^{2}}{|H(u,v)|^{2}+K}]{G(u,v)} \end{aligned}
F(u,v)^=[H(u,v)1∣H(u,v)∣2+K∣H(u,v)∣2]G(u,v)
import cv2
import numpy as np
from matplotlib import pyplot as plt
# 读取图像
img = cv2.imread('noise2.png',0)
# 进行傅里叶变换
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
#-------------大气湍流模型-----------------
k = 0.001 #湍流程度
degeneration = np.zeros(fshift.shape,dtype=complex) #要指定类型
for u in range(fshift.shape[0]):
for v in range(fshift.shape[1]):
H = np.exp(-k*np.power(np.float64((u+fshift.shape[0]//2)**2+(v-fshift.shape[1]//2)**2),5/6))
degeneration[u][v] = H
#-----------------------------------------
# 维纳滤波
K = 0.001
degeneration += 0.1 #加上估计的噪声
img1 = (np.conj(degeneration)/(np.conj(degeneration)*degeneration + K)) * fshift
# 将频域图像在空域进行显示做的处理
img1 = np.fft.ifftshift(img1)
img1 = np.fft.ifft2(img1)
img1 = np.abs(img1)
img1 = cv2.normalize(img1, None, 0, 1, cv2.NORM_MINMAX)
# 显示原始图像和恢复后的图像
plt.figure(figsize=(20, 20))
plt.subplot(121), plt.imshow(img, cmap='gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122), plt.imshow(img1, cmap='gray')
plt.title('rec img'), plt.xticks([]), plt.yticks([])
plt.show()
滤波效果：
总结
主要展示了维纳滤波和逆滤波两种方法，但是滤波效果并不明显，可能是估计的退化函数与真实的退化函数不一致导致的。文章来源：https://www.toymoban.com/news/detail-763095.html
参考文献
[1] 阮秋琦，阮宇智译；（美）拉斐尔·C.冈萨雷斯，理查德·E.伍兹.国外电子书与通信教材系列 数字图像处理 第4版[M].北京：电子工业出版社,2020文章来源地址https://www.toymoban.com/news/detail-763095.html
到了这里，关于图像处理---逆滤波和维纳滤波的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
![[笔记]Python计算机视觉编程《一》 基本的图像操作和处理](https://imgs.yssmx.com/Uploads/2024/02/433920-1.png)
