文章目录

前言
一、逆滤波
- 1.1 估计退化函数 H ( u , v ) H(u,v) H(u,v)
- - 1.1.1 观察法
  - 1.1.2 试验法
  - 1.1.3 建模法 ★ \bigstar ★
- 1.2 直接逆滤波
- 1.3 半径受限逆滤波
二、最小均方误差（维纳）滤波
总结
参考文献

前言

本文主要介绍退化图像复原的两种方法：逆滤波和维纳滤波。

一、逆滤波

图像退化的表达式：
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)**2kernel = 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.piif uv == 0:degeneration_img[u][v] = fshift[u][v]else:H = T * np.sin(uv) * np.exp(np.complex64(-1j)*uv) /uvdegeneration_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()

滤波效果：

总结

主要展示了维纳滤波和逆滤波两种方法，但是滤波效果并不明显，可能是估计的退化函数与真实的退化函数不一致导致的。

参考文献

[1] 阮秋琦，阮宇智译；（美）拉斐尔·C.冈萨雷斯，理查德·E.伍兹.国外电子书与通信教材系列数字图像处理第4版[M].北京：电子工业出版社,2020

图像处理---逆滤波和维纳滤波相关推荐

图像处理：图像复原与重建之逆滤波、维纳滤波、约束最小二乘滤波——Matlab实现
参考资料: 陷波滤波器-matlab实现 http://blog.sina.com.cn/s/blog_ebd29d830102wdzw.html 图像复原之约束最小二乘方滤波 https://blo ...
数字图像处理实验（五）|图像复原{逆滤波和伪逆滤波、维纳滤波deconvwnr、大气湍流扰动模型、运动模糊处理fspecial}（附matlab实验代码和截图）
文章目录一.实验目的二.实验仪器三.实验原理四.实验内容 1.逆滤波:选择MATLAB文件夹中的foggy图像作为实验图像. (1)生成退化函数: (2)复原 (a)直接逆滤波 (b)修正函数 ...
图像处理之仿真运动模糊复原【使用逆滤波、维纳滤波】
目录一.基础知识二.问题分析三.效果图四.代码一.基础知识图像去模糊是一个经典的图像复原任务.造成图像模糊的原因有很多,可以主要分为三大类离焦模糊:场景中的物体处于成像景深范围之外而变得 ...
图象恢复——（逆滤波，维纳滤波）
目的:对获取图像在频域用高斯函数进行退化并叠加白噪声,对退化图像进行逆滤波和维纳滤波恢复,比较原始图像和恢复图像,对利用逆滤波和维纳滤波恢复方法恢复图像进行比较. 一.基本原理图像复原是一种客观的操 ...
逆滤波和维纳滤波（附Matlab完整代码）
一.实验目的利用逆滤波和维纳滤波,对Lena加噪运动模糊降质图像进行复原,比较不同参数选择对复原结果的影响. 二.实验内容 1) 输入Lena图像,对图像进行运动降质:降质模型: 2) ...
数字图像处理——图像退化（大气湍流模型与运动模糊模型）与图像复原（逆滤波与维纳滤波）
一.图像退化一般来说,图像的退化模型可以表示为其中g(x,y) 表示退化后的图像,h(x,y)表示退化模型,f(x,y)表示原图像,n(x,y)表示噪声. 在频域上面可以表示为下面介绍常见的两种 ...
数字图像处理4：逆滤波及其变形和维纳滤波
逆滤波和维纳滤波简介在图像的获取.传输以及记录保存过程中,由于各种因素,如成像设备与目标物体的相对运动,大气的湍流效应,光学系统的相差,成像系统的非线性畸变,环境的随机噪声等原因都会使图像产生一定 ...
图像处理/计算机视觉/ python环境下如何用滤波器（/逆滤波/均值滤波/低通滤波/高通滤波）处理图片【附代码】
计算机视觉滤波器实操基础知识一. 计算机视觉技术中常见的几种滤波器二.滤波器相关知识应用一:算术均值.几何均值.谐波逆谐波一.问题分析二.结果图三.代码附录应用二:维纳滤波,逆滤波一 ...
4.3 Python图像处理之图像恢复-无约束滤波器（逆滤波）、有约束滤波器（维纳滤波器）
4.3 Python图像处理之图像恢复-无约束滤波器(逆滤波).有约束滤波器(维纳滤波器) 文章目录 4.3 Python图像处理之图像恢复-无约束滤波器(逆滤波).有约束滤波器(维纳滤波器) 1 算 ...

图像处理---逆滤波和维纳滤波