59.离焦去模糊滤波器（1）-理论＋demo篇

本文作者：小嗷

微信公众号：aoxiaoji

链接：https://f600lt.github.io/archives/

在本篇中，您将学习:

如何恢复模糊的图像？

什么是Wiener过滤器？

本文你会找到以下问题的答案:

Matlab的相关代码，还有去模糊的理论知识

2.1 恢复散焦和模糊的图像

这一解释是基于[72]和[231]两本书。此外，您可以参考Matlab的Matlab教程图像去模糊在Matlab和文章SmartDeblur。

[72]拉斐尔·C·冈萨雷斯和其他人。数字图像基础，数字图像处理。1987年。(冈萨雷斯)

[231]И.С. Грузман, В.С. Киричук, В.П. Косых, Г.И. Перетягин, and А.А. Спектор. Цифровая обработка изображений в информационных системах. 2000.

从理论和实践的角度来看，失真图像的恢复是图像处理中最有趣和最重要的问题之一。

有一些特殊的情况:由于不正确的焦点而产生的模糊和由于运动而产生的模糊——这些缺陷(你们每个人都很清楚，并且很难修复)被选为本文的主题。至于其他的图像缺陷(噪音，不正确的曝光，失真)，人类已经学会了如何纠正它们，任何好的照片编辑器都有这个工具。

为什么几乎没有办法校正模糊和散焦(除了不锋利的Mask)-也许根本不可能做到这一点?事实上，这是有可能的——相关的数学理论的发展大约始于70年前，但像其他图像处理算法一样，去模糊算法12年才得到广泛应用。下面是几张照片，展示了这种令人惊叹的效果:

右边的图像实际上是由左边的图像构成的，而且，没有使用48位格式的优化(在这种情况下，几乎100%的恢复源图像)——只有一个规则的PNG，在左边有合成模糊。结果令人印象深刻……但我们在实践中，并非一切都那么好。

2.2 简介

让我们从远处开始。很多人认为模糊是一种不可逆转的操作，这种情况下的信息会永远丢失，因为每个像素都变成了一个点，所有的东西都混合在一起，在一个大的模糊半径的情况下，我们会得到一个平面的颜色在图像上。但这并不是完全正确的——所有的信息只是根据某些规则重新分配，并且在某些假设下肯定可以恢复。一个例外是图像的边界，它的宽度等于模糊半径-这里不可能完全恢复。

让我们用一个一维情况的小例子来演示它——假设我们有一行像素，具有以下值:

x1 | x2 | x3 | x4…——源图像

在模糊每个像素的值后，将其添加到左边的值:X'i = Xi + Xi-1。通常，它也需要除以2，但为了简单起见，我们将把它去掉。因此，我们得到一个模糊的图像与以下像素值:

x1 + x0 | x2 + x1 | x3 + x2 | x4 + x3…——模糊图像

现在我们将尝试恢复它,我们将必然地减去值根据以下方案——第一个像素从第二个,第二个像素的结果从第三个,第三从第四个像素的结果等等,我们将得到以下:

x1 + x0 | x2 - x0 | x3 + x0 | x4 - x0…——恢复图像

因此，我们得到的不是模糊的图像，而是在每个像素上增加一个未知的常数x0和交替符号的源图像。

这是更好的,我们可以选择这个常数视觉上,我们可以假设它约等于x1,与这样一个标准,我们可以自动选择相邻的像素值变化尽可能小,等等。（总之是有规律的）

但一切都改变当我们添加噪声(总是存在在现实图像)。

对于所描述的方案，在每个阶段都会将噪声值累加到总值中，这个事实最终会产生一个绝对不可接受的结果，但是我们已经知道，即使使用这样一个原始的方法，恢复也是很可能的。

2.3 模糊过程模型

现在让我们来更正式和科学地描述这些模糊和恢复过程。我们将只考虑灰度图像，假设对于全彩色图像的处理，对每个RGB颜色通道重复所有需要的步骤就足够了。让我们介绍以下定义:

f(x, y) -源图像(非模糊)

h(x, y) -模糊函数

n(x, y) -加性噪声

g(x, y) -模糊结果图像

我们将通过以下方式形成模糊过程模型:

g(x, y) = h(x, y) * f(x, y) + n(x, y) (1)

模糊图像恢复的任务是找到源图像的最佳逼近f'(x, y)。让我们以更详细的方式考虑每个组件。对于函数f(x, y)和g(x, y)所有的东西都很清楚。但是对于h(x, y)我需要说几个词，它是什么?在模糊的过程中，源图像的每个像素都变成了散焦时的一个点，变成了通常由于运动而模糊的线段(或路径)。或者我们可以说，模糊图像的每个像素都是由源图像附近的像素“组装”而成的。

所有这些重叠在一起，导致图像模糊。根据这个原理，一个像素的传播被称为模糊函数。其他同义词- PSF (Point spread function)、kernel等。（PSF：点散播函数，小嗷等一下会说到）

这个函数的大小小于图像本身的大小——例如，当我们考虑第一个“演示”示例时，函数的大小是2，因为每个结果像素由两个像素组成。

2.4 模糊函数

让我们看看典型的模糊函数是什么样的。下面我们将使用已经成为这方面标准的工具——Matlab，它包含了图像处理最多样化实验所需的所有东西(以及其他东西)，并允许专注于算法，将所有日常工作转移到函数库。然而，这只能以牺牲性能为代价。让我们回到PSF，下面是他们的例子:

对于高斯函数的使用PSF:产生预定义的滤波器('gaussian', 30, 8);

对于运动模糊函数使用:产生预定义的滤波器(‘motion’，40,45);

呵呵哒

将模糊函数应用到另一个函数(在他的例子中，是一个图像)的过程称为卷积，即源图像的某个区域卷积到模糊图像的一个像素。它通过运算符“*”表示，但不要将其与简单的乘法混淆!从数学上讲，对于维数为M×N的图像f和维数为M×N的模糊函数h，可以写成:

其中a = (m - 1) / 2, b = (n - 1) / 2。与卷积相反的过程被称为反卷积，这种任务的求解相当罕见。

2.5 噪声模型

我们只考虑最后一个求和项，即(1)式中引起噪声的n(x, y)，数字传感器产生噪声的原因有很多，但最基本的原因是热振动(布朗运动)和暗电流。噪音量还取决于许多因素,如ISO值,传感器类型,像素尺寸,温度,磁场值,等等。在大多数情况下有高斯噪声(这是由两个参数——平均值和色散),也是添加剂,与图像不相关,不取决于像素坐标。

最后三个假设对于进一步的工作非常重要。

小嗷有一篇专门写噪声第23篇

2.6 卷积定理

让我们回到恢复的最初任务——我们需要考虑到噪声，以某种方式逆转卷积。从公式(2)中我们可以看出，从g(x, y)中得到f(x, y)并不容易，如果我们直接计算它，我们就会得到一组庞大的方程。但是

傅里叶

变换起了作用，我们不会详细讨论它，因为关于这个话题已经讲了很多了。

因此，有一个所谓的卷积定理，根据这个定理，卷积在空间域的运算等价于在频域的常规乘法(这里乘法是一个元素一个元素，而不是矩阵1)。与卷积相反的运算相当于频域中的除法，即可以表示为:

（终于说到点子上，卷积在空间域的运算等价于在频域的常规乘法）

其中H(u, v)、 F(u, v) ---- H(x,y)和F(x,y)的傅里叶函数。因此，(1)的模糊过程可以在频域中表示为:

2.7 逆滤波器

我们想把这个方程除以H(u,v)和得到以下评价F ^(u,v)的源图像:

逆向运算

这被称为逆滤波，但在实践中它几乎从未奏效。为什么如此?为了回答这个问题，让我们看看公式(5)中的最后一个求和项——如果函数H(u, v)给出的值接近于零或等于它，那么这个求和项的输入将占主导地位。这在实际例子中几乎都能看到为了解释这一点让我们记住傅里叶变换之后的频谱是什么样子的。

源图：

用Matlab将它转换成灰度图，得到光谱:

因此，我们将有两个分量:振幅和相位谱。顺便说一下，很多人忘记了这个阶段。请注意，振幅谱是用对数尺度表示的，因为它的值变化非常大——有几个数量级，在中心的值是最大值(数百万)，当它们离中心越来越远时，很快就会衰减到零。由于这一事实，反滤波将只在零或几乎零噪声值的情况下工作。让我们在实践中演示一下，使用以下脚本:（简单来说爆表，超过255或者低过1）

noisevar = 0.0000001 noisevar = 0.000005

noise_var为方差

显而易见，即使是非常小的噪声也会造成严重的失真，这实际上大大限制了该方法的应用。

2.8 现有的去卷积方法

有一些方法，考虑到图像中存在的噪声，其中最流行也是最早的方法之一是Wiener滤波器。将图像和噪声视为随机过程，对无畸变的图像f取f'，这些值的均方差最小。该偏差的最小值在频域函数上达到:

这个结果是Wiener在1942年发现的。在本文中我们不会给出他的详细结论，有兴趣的人可以在这里找到。S函数在这里分别表示噪声的能谱和源图像的能谱——由于这些值很少为人所知，那么部分Sn / Sf被某个常数K所取代，可以近似地表示为信噪比。

下一种方法是“约束最小二乘滤波”，其他名称为“Tikhonov滤波”、“Tikhonov正则化”。他的思想是将任务以矩阵的形式形成，然后求解相应的优化任务。这个方程的结果可以写成:

在y -正则化参数,аP(u,v)——Fourier-function拉普拉斯算子矩阵(3 * 3)。

另一个有趣的方法是Richardson(1972年)和Lucy独立(1974年)提出的，所以这种方法被称为Lucy-Richardson方法。它的独特之处在于它是非线性的，这与前三种不同——这可能会得到更好的结果。第二个特征——这种方法是迭代的，因此迭代停止的准则存在困难。

主要思想在于使用最大似然法，假定图像受到泊松分布的影响。计算公式非常简单，不使用傅里叶变换，一切都是在空间域内完成的:

这里的符号“*”，如前所述，表示卷积运算。这种方法在处理天文照片的程序中被广泛使用——使用反褶积(而不是像照片编辑器那样使用粗糙的掩模)是事实上的标准。该方法的计算复杂度相当高——根据迭代次数的不同，对一张普通照片进行处理需要花费数小时甚至数天的时间。

最后考虑的方法，确切地说，目前还没有开发的所有方法都是盲目反褶积。在之前的所有方法中，我们都认为模糊函数PSF是已知的，但在实践中，这是不正确的，通常我们只知道近似的PSF通过可见扭曲的类型。盲目反褶积就是要考虑到这一点。

原则非常简单,没有深入细节——有选择的第一个近似的PSF,然后使用的一个方法,执行反褶积后的质量是识别的程度根据一些标准,基于此学位PSF函数优化和迭代重复直到达到所需的结果。

2.9 实践

让我们继续练习，我们将从综合模糊和噪声图像上列出的方法的比较开始。

结果图：

15.png

16.PNG

3.0 结论

在此之前，所有的模糊都是人为的，这对实践和学习是很好的，但看看这些如何与真实的照片一起工作是非常有趣的。这里有一个这样的例子，用佳能500D相机用手动对焦拍摄(模糊):

然后运行脚本:

得到以下结果:

正如我们所看到的，图像上出现了新的细节，清晰度变得更高了，但是在锐利的边缘上出现了“环效应”。

还有一个由于运动而产生的真实模糊的例子——为了使相机固定在三脚架上，设置了相当长的曝光值，并且在曝光时的移动是均匀的，得到如下模糊:

19.png

脚本几乎是一样的，但是PSF类型是“motion”:

Result:

20.png

再一次显著提高了质量，汽车变得容易辨认。现在不同了之处，不像上一个散焦的例子。

3.1 程序exe：

这是小嗷在github找到的，演示了模糊和散焦图像的恢复。选择了FFTW库，它是傅立叶变换最快的开源实现。

程序被称为SmartDeblur, windows二进制文件和来源，你可以从GitHub:

https://github.com/Y-Vladimir/SmartDeblur/downloads

所有源文件都在GPL v3许可证之下。

主窗口截图:

主要功能

高速度。2048*1500像素的图像处理在预览模式下需要300毫秒左右(调整滑块可以移动时)。但高质量的加工可能需要几分钟

实时参数变化应用(没有任何预览按钮)

全分辨率处理(没有小预览窗口)

内核参数的深度调优

简单友好的用户界面

帮助屏幕与图像示例

卷积方法：Wiener, Tikhonov, Total Variation prior

在结束语中，我再举一个实际(非合成)焦距模糊的例子: