【youcans 的 OpenCV 例程 200 篇】106. 退化图像的逆滤波
6. 退化图像复原
图像复原是对图像退化的过程进行估计，并补偿退化过程造成的失真，以便获得未经退化的原始图像或原始图像的最优估值，从而改善图像质量的一种方法。
典型的图像复原方法是根据图像退化的先验知识建立退化模型，以退化模型为基础采用滤波等手段进行处理，使复原后的图像符合一定的准则，达到改善图像质量的目的。
因此，图像复原是沿着质量降低的逆过程来重现真实的原始图像，通过去模糊函数而去除图像模糊。
6.1 退化图像的逆滤波（Inverse filter）
图像退化表示为退化算子
H
\mathcal{H}
H， 退化函数可以用观察法、试验法或建模法估计，则通过逆滤波就可以直接实现图像复原。用退化图像的傅里叶变换除以退化函数的傅里叶变换，得到原始图像的傅里叶变换估计：
F
^
(
u
,
v
)
=
G
(
u
,
v
)
H
(
u
,
v
)
\hat{F}(u,v) = \frac{G(u,v)}{H(u,v)}
F^(u,v)=H(u,v)G(u,v)
但是，由于实际上退化图像是退化算子与加性噪声项共同作用的结果，由此得到：
F
^
(
u
,
v
)
=
F
(
u
,
v
)
+
N
(
u
,
v
)
H
(
u
,
v
)
\hat{F}(u,v) = F(u,v) + \frac{N(u,v)}{H(u,v)}
F^(u,v)=F(u,v)+H(u,v)N(u,v)
这表明即使获得退化函数 H ( u , v ) H(u,v) H(u,v) 的估计，由于噪声项是未知的，因此也不能准确地复原原始图像。
进一步地，如果退化函数为 0 或很小，则噪声项的影响将非常严重（信噪比低）。这时，需要将频率限制到原点附近进行分析，可以减少遇到零值的可能性。
例程 9.20：湍流模糊退化图像的逆滤波
如前所述，通过湍流退化模型可以得到退化图像。使用该退化模型进行逆滤波，退化函数与生成退化图像所用的退化函数相反：
H
(
u
,
v
)
=
e
−
k
[
(
u
−
M
/
2
)
2
+
(
v
−
N
/
2
)
2
]
5
/
6
H(u,v) = e^{-k [(u-M/2)^2+(v-N/2)^2]^{5/6}}
H(u,v)=e−k[(u−M/2)2+(v−N/2)2]5/6
但是，直接使用退化模型 H(u,v) 逆滤波的结果（D0=full）很差，用理想低通滤波器对退化模型 H(u,v) 在半径 D0 之外截止后，则视觉效果较好。
# 9.20: 湍流模糊退化图像的逆滤波
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
def getDegradedImg(image, Huv, eps): # 根据退化模型生成退化图像
# (1) 傅里叶变换, 中心化
fft = np.fft.fft2(image.astype(np.float32)) # 傅里叶变换
fftShift = np.fft.fftshift(fft) # 将低频分量移动到频域图像中心
# (2) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数
fftShiftFilter = fftShift * Huv # Guv = Fuv * Huv
# (3) 对修正傅里叶变换 进行傅里叶逆变换，逆中心化
invShift = np.fft.ifftshift(fftShiftFilter) # 将低频分量逆转换回图像四角
imgIfft = np.fft.ifft2(invShift) # 逆傅里叶变换，返回值是复数数组
imgDegraded = np.uint8(cv2.normalize(np.abs(imgIfft), None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]
return imgDegraded
def ideaLPFilter(img, radius=10): # 理想低通滤波器
M, N = img.shape[1], img.shape[0]
u, v = np.meshgrid(np.arange(M), np.arange(N))
D = np.sqrt((u - M//2)**2 + (v - N//2)**2)
kernel = np.zeros(img.shape[:2], np.float32)
kernel[D <= radius] = 1
return kernel
def inverseFilter(image, Huv, D0): # 根据退化模型逆滤波
# (1) 傅里叶变换, 中心化
fft = np.fft.fft2(image.astype(np.float32)) # 傅里叶变换
fftShift = np.fft.fftshift(fft) # 将低频分量移动到频域图像中心
# (2) 在频率域修改傅里叶变换: 傅里叶变换 点乘 滤波器传递函数
if D0==0:
fftShiftFilter = fftShift / Huv # Guv = Fuv / Huv
else:
lpFilter = ideaLPFilter(image, radius=D0)
fftShiftFilter = fftShift / Huv * lpFilter # Guv = Fuv / Huv
# (3) 对修正傅里叶变换 进行傅里叶逆变换，逆中心化
invShift = np.fft.ifftshift(fftShiftFilter) # 将低频分量逆转换回图像四角
imgIfft = np.fft.ifft2(invShift) # 逆傅里叶变换，返回值是复数数组
imgRebuild = np.uint8(cv2.normalize(np.abs(imgIfft), None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]
return imgRebuild
# 读取原始图像
img = cv2.imread("../images/Fig0525a.tif", 0) # flags=0 读取为灰度图像
# 生成湍流模糊图像
HTurb = turbulenceBlur(img, k=0.0025)
imgBlur = np.abs(getDegradedImg(img, HTurb, 0.0))
print(imgBlur.max(), imgBlur.min())
# # 逆滤波
imgRebuild = inverseFilter(imgBlur, HTurb, 480) # Huv 全滤波器
imgRebuild1 = inverseFilter(imgBlur, HTurb, D0=40) # 在半径 D0 之外 Huv 截止
imgRebuild2 = inverseFilter(imgBlur, HTurb, D0=70)
imgRebuild3 = inverseFilter(imgBlur, HTurb, D0=100)
plt.figure(figsize=(9, 7))
plt.subplot(231), plt.title("origin"), plt.axis('off'), plt.imshow(img, 'gray')
plt.subplot(232), plt.title("turbulence blur"), plt.axis('off'), plt.imshow(imgBlur, 'gray')
plt.subplot(233), plt.title("inverse filter(D0=full)"), plt.axis('off'), plt.imshow(imgRebuild, 'gray')
plt.subplot(234), plt.title("inverse filter(D0=40)"), plt.axis('off'), plt.imshow(imgRebuild1, 'gray')
plt.subplot(235), plt.title("inverse filter(D0=70)"), plt.axis('off'), plt.imshow(imgRebuild2, 'gray')
plt.subplot(236), plt.title("inverse filter(D0=100)"), plt.axis('off'), plt.imshow(imgRebuild3, 'gray')
plt.tight_layout()
plt.show()
（本节完）
版权声明：
youcans@xupt 原创作品，转载必须标注原文链接：(https://blog.csdn.net/youcans/article/details/123027356)
Copyright 2022 youcans, XUPT
Crated：2022-2-20
