使用 SimpleITK 反卷积滤波器实现非盲去模糊图像反卷积是指一类试图从模糊图像中恢复清晰图像的算法，其使用模糊图像

开启掘金成长之旅！这是我参与「掘金日新计划 · 12 月更文挑战」的第29天，点击查看活动详情

前言

图像反卷积是指一类试图从模糊图像中恢复清晰图像的算法，其使用模糊图像和可能的卷积核 h (假设卷积核是已知的或估计的)作为输入。如果模糊核 h 已知，则去模糊过程称为非盲取模糊 (non-blind deblurring)，反之称为盲去模糊 (blind deblurring)。在本节中，我们将学习如何使用 SimpleITK 库的反卷积滤波器执行非盲去模糊。

实现非盲去模糊

(1) 首先，读取输入图像，并使用失焦模糊核模糊该图像：

import SimpleITK as sitk
from scipy import signal
from skimage.metrics import peak_signal_noise_ratio
from skimage.io import imread
from skimage.color import rgb2gray
import numpy as np
import matplotlib.pylab as plt
import cv2

def get_out_of_focus_kernel(r, sz=15):
    kern = np.zeros((sz, sz), np.uint8)
    cv2.circle(kern, (sz, sz), r, 255, -1, cv2.LINE_AA, shift=1)
    kern = np.float32(kern) / 255
    return kern

im = rgb2gray(imread('1.png'))
psf = get_out_of_focus_kernel(7, 9).astype(np.float)
im_blur = signal.convolve2d(im, psf, 'same')
im_blur = im_blur / np.max(im_blur)

(2) 使用逆滤波器恢复图像是最直接的反卷积方法。根据卷积定理，空间域中的两个图像的卷积等价于两个图像的傅立叶变换乘积，因此逆滤波器需要对乘法求逆。换句话说，逆滤波器计算公式如下所示：

\hat F(u,v)=\left\{ \begin{array}{rcl} \frac {G(u,v)} {H(u,v)} & & {|H(u,v)|≥\epsilon}\\ 0 & & {otherwise} \end{array} \right.

(3) 使用 SimpleITK.InverseDeconvolutionImageFilter() 函数实现直接线性反卷积滤波器：

tkfilter = sitk.InverseDeconvolutionImageFilter()
tkfilter.SetNormalize(True)
im_res_IN = sitk.GetArrayFromImage(tkfilter.Execute (sitk.GetImageFromArray(im_blur), sitk.GetImageFromArray(psf)))

(4) 使用 Wiener 反卷积图像滤波器恢复与模糊核卷积后得到的图像，同时将噪声降至最低。

Wiener 滤波器旨在最大程度地减少低信噪比 (signal-to-noise ratio, SNR) 频率引起的噪声。Wiener 滤波器核是在频域中定义的，如下公式所示( SNR 越高，恢复的图像质量越高)：

W(u,v)=\frac {H^*(u,v)} {|H(u,v)|^2+K(u,v)}=\frac {H^*(u,v)} {|H(u,v)|^2+\frac 1 {SNR(u,v)}}\\ K(u,v)=\frac {S_{\eta}(u,v)}{S_f(u,v)}

使用 SimpleITK.WienerDeconvolutionImageFilter() 函数实现 Wiener 滤波器。在这里，模糊后未向原始图像添加额外的噪音。因此，噪声方差设置为零：

tkfilter = sitk.WienerDeconvolutionImageFilter()
tkfilter.SetNoiseVariance(0)
tkfilter.SetNormalize(True)
im_res_WN = sitk.GetArrayFromImage(tkfilter.Execute (sitk.GetImageFromArray(im_blur), sitk.GetImageFromArray(psf)))

(5) 使用 Tikhonov 方法恢复图像。该方法通过在分母中添加正则化项来改进估计的反卷积滤波器(去模糊)，该滤波器最小化以下公式中：

\underset {\hat f} {min} ||\hat f \otimes h-g||^2_{L_2}+ \mu ||\hat f||^2 \\ T(u,v)=\frac {H^*(u,v)} {|H(u,v)|^2+\mu}

其中， $\hat f$ 表示去模糊图像 $f$ 的估计值； $h$ 表示模糊核； $g$ 表示模糊图像；项 $μ$ 在此滤波器中称为正则化，如果将μ设置为零，则该滤波器等价于逆滤波器。以下代码使用 simpleitk.tikhonovdeconvolutionImageFilter() 函数实现 Tikhonov 方法正则化的逆反卷积滤波器。

其中，将正则化常数设置为 0.008，我们也可以尝试不同的值，并观察不同值到对去模糊输出图像的影响：

tkfilter = sitk.TikhonovDeconvolutionImageFilter()
tkfilter.SetRegularizationConstant(0.008) #0.06)
tkfilter.SetNormalize(True)
im_res_TK = sitk.GetArrayFromImage(tkfilter.Execute (sitk.GetImageFromArray(im_blur), sitk.GetImageFromArray(psf)))

(6) 使用 Richardson-Lucy 算法恢复图像，该滤波器实现了 Richardson-Lucy 反卷积算法。该算法假定输入图像是由具有已知核的线性移位不变系统形成的。迭代算法使用贝叶斯的方法来改善每次迭代的估计（去模糊）图像：

f_{i+1}(x)=\{[\frac {c(x)} {f_i(x)\otimes g(x)}]\otimes g(-x)\}f_i(x)

其中 $\otimes$ 表示卷积操作， $g(x)$ 表示 PSF， $c(x)$ 表示模糊图像， $f_i(x)$ 表示去模糊图像。

以下代码使用 SimpleItk.RichardsonLucyDeconVolutionImageFilter() 函数利用 Richardson-Lucy 反卷积算法对输入图像执行反卷积。我们可以尝试使用不同的迭代次数，并观察对输出图像的影响：

tkfilter = sitk.RichardsonLucyDeconvolutionImageFilter()
tkfilter.SetNumberOfIterations(100)
tkfilter.SetNormalize(True)
im_res_RL = sitk.GetArrayFromImage(tkfilter.Execute (sitk.GetImageFromArray(im_blur), sitk.GetImageFromArray(psf)))

(7) 最后，绘制所有原始、模糊和使用不同反卷积算法得到的去模糊图像并计算 PSNR 值(以比较使用不同反卷积算法得到的去模糊图像的质量)：

plt.figure(figsize=(20, 60))
plt.subplots_adjust(0,0,1,1,0.07,0.07)
plt.gray()

plt.subplot(231), plt.imshow(im), plt.axis('off'), plt.title('Original Image', size=10)
plt.subplot(232), plt.imshow(im_blur), plt.axis('off'), plt.title('Blurred (out-of-focus) Image, PSNR={:.3f}'.format(peak_signal_noise_ratio(im, im_blur)), size=10)
plt.subplot(233), plt.imshow(im_res_IN, vmin=im_blur.min(), vmax=im_blur.max()), plt.axis('off')
plt.title('Deconvolution using SimpleITK (Inverse Deconv.), PSNR={:.3f}'.format(peak_signal_noise_ratio(im, im_res_IN)), size=10)
plt.subplot(234), plt.imshow(im_res_WN, vmin=im_blur.min(), vmax=im_blur.max()), plt.axis('off')
plt.title('Deconvolution using SimpleITK (Wiener Deconv.), PSNR={:.3f}'.format(peak_signal_noise_ratio(im, im_res_WN)), size=10)
plt.subplot(235), plt.imshow(im_res_RL, vmin=im_blur.min(), vmax=im_blur.max()), plt.axis('off')
plt.title('Deconvolution using SimpleITK (Richardson-Lucy), PSNR={:.3f}'.format(peak_signal_noise_ratio(im, im_res_RL)), size=10)
plt.subplot(236), plt.imshow(im_res_TK, vmin=im_blur.min(), vmax=im_blur.max()), plt.axis('off')
plt.title('Deconvolution using SimpleITK (Tikhonov Deconv.), PSNR={:.3f}'.format(peak_signal_noise_ratio(im, im_res_TK)), size=10)

plt.show()