模糊后恢复原始图像

Question

我有以下代码将高斯滤波器应用于任意图像 25 次。每次应用过滤器时，都会对生成的图像进行归一化。

kernel = np.array([[1.0,2.0,1.0],
                  [2.0,4.0,2.0],
                  [1.0,2.0,1.0]])

for i in range(25):
    # handle each component separately
    img[:,:,0] = convolve(img[:,:,0], kernel, mode='same')
    img[:,:,1] = convolve(img[:,:,1], kernel, mode='same')
    img[:,:,2] = convolve(img[:,:,2], kernel, mode='same')
    img = img / 16 # normalize

扭转这一过程的最佳方法是什么？ IE。如果我有一个模糊的图像（执行上面的代码的结果）并且想要获得原始图像。

编辑 1：

例子

原文：

模糊：

编辑 2：

尝试复制克里斯的答案

我安装了dipimage_2.9。我正在使用macOS 10.14.2 和Matlab R2016a。

我花了一段时间才知道如何为卷积指定边界条件，因为 DIPimage 的 convolve.m 只接受 image_in 和 kernel 参数。我最终为此使用了dip_setboundary（DIPimage User Manual 第 9.2 节）。

这是代码（我只是相应地添加了dip_setboundary 和cut 的裁剪区域的起源）：

% Get data
a = readim('https://i.stack.imgur.com/OfSx2.png'); % using local path in real code
a = a{1}; % Keep only red channel

%% Create kernel
kernel = [1.0,2.0,1.0
          2.0,4.0,2.0
          1.0,2.0,1.0] / 16;
tmp = deltaim((size(kernel)-1)*25+1);
dip_setboundary('add_zeros');
for ii=1:25
    tmp = convolve(tmp,kernel);
end
kernel = tmp;

%% Apply convolution
dip_setboundary('periodic');
b = convolve(a,kernel);
dip_setboundary('symmetric'); % change back to default
% Find inverse operation
%   1- pad stuff so image and kernel have the same size
%      we maintain the periodic boundary condition for image b
b = repmat(b,ceil(imsize(kernel)./imsize(b)));
kernel = extend(kernel,imsize(b));
%   2- apply something similar to Wiener deconvolution
c = real(ift(ft(b)/(ft(kernel)+1e-6))); % Not exactly Wiener, but there's no noise!
%   3- undo padding
c = cut(c,imsize(a), [0, 0]); % upper left corner

这是生成的图像c：

Answer 1

我们看问题中单通道的代码，假设img是灰度图像——这里的所有东西都可以按通道应用，所以我们不需要重复3次：

for i in range(25):
    img = ndimage.convolve(img, kernel)
    img = img / 16 # normalize

我们将在一分钟内撤消卷积。首先让我们简化应用的操作。

简化处理

以上与以下内容相同（在数值精度范围内）：

kernel = kernel / 16 # Normalize
for i in range(25):
    img = ndimage.convolve(img, kernel)

只要img 不是发生剪裁和/或舍入的整数类型，这是正确的。一般来说，* 是卷积，C 是一些常数，

g = C (f * h) = f * (C h)

接下来，我们知道应用卷积 25 次与使用复合内核应用卷积一次是一样的，

g = (((f * h) * h) * h) * h = f * (h * h * h * h)

我们如何获得复合内核？将卷积应用于全为零且中间像素为 1 的图像再次生成内核，因此

delta = np.zeros(kernel.shape)
delta[delta.shape[0]//2, delta.shape[1]//2] = 1
kernel2 = ndimage.convolve(delta, kernel)
kernel2 == kernel # is true everywhere, up to numerical precision

因此，以下代码在问题中找到用于平滑图像的内核：

kernel = np.array([[1.0,2.0,1.0],
                  [2.0,4.0,2.0],
                  [1.0,2.0,1.0]]) / 16
delta = np.zeros(((kernel.shape[0]-1)*25+1, (kernel.shape[1]-1)*25+1))
delta[delta.shape[0]//2, delta.shape[1]//2] = 1
for i in range(25):
    delta = ndimage.convolve(delta, kernel)
kernel = delta

由于中心极限定理，这个核与高斯核非常相似。

现在我们可以通过单个卷积获得与问题中相同的输出：

output = ndimage.convolve(img, kernel)

反转卷积

逆滤波的过程称为反卷积。理论上这是一个非常微不足道的过程，但实际上由于噪声、对内核的不精确知识等原因，这非常困难。

我们知道我们可以通过傅里叶域计算卷积：

output = np.convolve(img, kernel, mode='wrap')

与

相同

output = np.real(np.fft.ifft2( np.fft.fft2(img) * np.fft.fft2(np.fft.ifftshift(kernel)) ))

（假设kernel 与img 大小相同，我们通常必须先用零填充它）。空间域和频域运算结果之间的任何差异都是由使用convolve 时图像超出其边界的方式造成的。傅里叶方法假设一个周期性边界条件，这就是为什么我在这里使用'wrap' 模式进行卷积。

逆运算只是傅里叶域中的一个除法：

img = np.real(np.fft.ifft2( np.fft.fft2(output) / np.fft.fft2(np.fft.ifftshift(kernel)) ))

为此，我们需要知道kernel 的确切值，并且在此过程中不应添加任何噪音。对于如上计算的output，理论上应该给出准确的结果

但是，对于某些频率分量，某些内核可能正好为零（即 np.fft.fft2(np.fft.ifftshift(kernel)) 包含零）。这些频率无法恢复，除以 0 会导致 NaN 值在逆变换中传遍整个图像，逆变换的图像将全是 NaN。

对于高斯核没有零，所以这不应该发生。然而，会有许多频率非常接近于零。因此，output 的傅里叶变换对于这些元素也将具有非常小的值。然后，逆过程是将一个非常小的值除以另一个非常小的值，从而导致数值精度问题。

而且你可以看到这个过程如何，如果只有一点点噪音，就会大大增强这种噪音，使得输出几乎完全由这种噪音提供。

Wiener 反卷积包括正则化，以防止这些带有噪声和数值不精确的问题。基本上，您可以通过向kernel 的傅里叶变换添加一个正值来防止除以非常小的数。 Wikipedia 很好地描述了维纳反卷积。

演示

我在这里使用带有 DIPimage 3 的 MATLAB 来做一个快速演示（对我来说比启动 Python 并弄清楚如何在那里完成所有这些工作要少得多）。这是代码：

% Get data
a = readim('https://i.stack.imgur.com/OfSx2.png');
a = a{1}; % Keep only red channel

% Create kernel
kernel = [1.0,2.0,1.0
          2.0,4.0,2.0
          1.0,2.0,1.0] / 16;
tmp = deltaim((size(kernel)-1)*25+1);
for ii=1:25
    tmp = convolve(tmp,kernel,'add zeros');
end
kernel = tmp;

% Apply convolution
b = convolve(a,kernel,'periodic');

% Find inverse operation
%   1- pad stuff so image and kernel have the same size
%      we maintain the periodic boundary condition for image b
b = repmat(b,ceil(imsize(kernel)./imsize(b)));
kernel = extend(kernel,imsize(b));
%   2- apply something similar to Wiener deconvolution
c = ift(ft(b)/(ft(kernel)+1e-6),'real'); % Not exactly Wiener, but there's no noise!
%   3- undo padding
c = cut(c,imsize(a),'top left');

这是输出，上三分之一是输入图像，中间三分之一是模糊图像，下三分之一是输出图像：

这里需要注意的重要一点是，我对初始卷积使用了周期性边界条件，这与傅里叶变换中发生的情况相匹配。其他边界条件将在边缘附近的逆变换中引起伪影。因为内核大小比图像大，所以整个图像将是一个大人工制品，您将无法恢复任何东西。另请注意，为了将内核用零填充到图像的大小，我必须复制图像，因为内核大于图像。复制图像再次匹配傅里叶变换强加的周期性边界条件。如果输入图像比卷积核大得多，这两个技巧都可以忽略，正如您在正常情况下所期望的那样。

还要注意，如果没有反卷积中的正则化，输出都是 NaN，因为我们将非常小的值除以非常小的值。内核的傅立叶变换中有很多接近零的地方，因为模糊非常严重。

最后，请注意，即使在模糊图像中添加少量噪点，也无法以可以阅读文本的方式对图像进行反卷积。逆变换看起来会很漂亮，但文字笔画会变形到足以使字母不再容易识别：

---

上面的代码使用了 DIPimage 3，它还没有官方的二进制文件可以安装，它需要从source 构建。要使用 DIPimage 2.x 运行代码，需要进行一些更改：

1. 必须使用dip_setboundary 设置边界条件，而不是直接将其传递给convolve 函数。字符串'add zeros' 和'periodic' 是边界条件。

2. ft 和 ift 函数使用对称归一化，每个函数将它们的输出乘以 1/sqrt(prod(imsize(image)))，而在 DIPimage 3 中，归一化是更常见的乘以 1/prod(imsize(image)) ift 和 1为ft。这意味着kernel的傅里叶变换必须乘以sqrt(prod(imsize(kernel)))才能匹配DIPimage 3的结果：

   c = real(ift(ft(b)/((ft(kernel)*sqrt(prod(imsize(kernel))))+1e-6)));
   

Answer 2

你不能 - 模糊通过平均失去信息。

考虑一个 1-dim 示例：

[1  2  1]   on    [1,2,3,4,5,6,7]  assuming 0 for missing "pixel" on convolution

结果为@​​987654323@。 8 可能来自[1,2,3]，也可能来自[2,1,4] - 所以你已经有两种不同的解决方案。无论您采用哪种方式，都会影响任何可能成为 12 来源的值。

这是一个过于简单的例子 - 你可以解决这个问题 - 但在图像处理中，你可能会处理 3000*2000 像素和 3x3、5x5、7x7 的二维卷积，...矩阵使得反转不切实际。

使这个二维你可能能够在数学上解决它 - 但更常见的是，如果你将它应用于二维，你不会得到大量的解决方案和非常复杂的约束来解决它卷积和 3000*2000 像素。