图像的反卷积

Question

我有源图像和结果图像。我知道，在源图像上使用了某些卷积矩阵来得到结果。是否可以计算出这个卷积矩阵？或者至少得到一个非常相似的卷积矩阵。

- Martin Perry

Answer 1

原则上是可以的。只需使用FFT将两个图像转换到频率空间中，然后将结果图像的FFT除以源图像的FFT即可得到卷积核的近似值。然后应用反FFT以获取卷积核的近似值。

要理解为什么会这样，需要注意在空间域中的卷积对应于频率域中的乘法，因此反卷积同样对应于频率域中的除法。在普通的反卷积中，人们会将卷积图像的FFT除以内核的FFT以恢复原始图像。但是，由于卷积（如乘法）是可交换的操作，内核和源的角色可以任意交换：通过内核卷积源与通过源卷积内核完全相同。

然而，正如其他答案所指出的那样，这不太可能得到精确重建卷积核的结果，原因与普通的反卷积一样，通常无法精确重建原始图像：四舍五入和裁剪将引入噪声，并且卷积有可能完全擦除某些频率（将它们乘以零），此时这些频率无法重建。

话虽如此，如果原始卷积核大小有限，则重建的卷积核通常应该在原点周围具有单个尖峰，逼近原始卷积核，被低级噪声所包围。即使您不知道原始卷积核的确切大小，提取这个尖峰并丢弃重建的其余部分也不应太难。

例如：这是灰度的Lenna测试图像，缩小到256×256像素，并在GIMP中使用5×5内核进行卷积的原图像。我将使用numpy/scipy和Python进行反卷积：

from scipy import misc
from numpy import fft

orig = misc.imread('lena256.png')
blur = misc.imread('lena256blur.png')
orig_f = fft.rfft2(orig)
blur_f = fft.rfft2(blur)

kernel_f = blur_f / orig_f         # do the deconvolution
kernel = fft.irfft2(kernel_f)      # inverse Fourier transform
kernel = fft.fftshift(kernel)      # shift origin to center of image
kernel /= kernel.max()             # normalize gray levels
misc.imsave('kernel.png', kernel)  # save reconstructed kernel

下方展示了卷积核的256×256像素重建图以及其中心7×7像素区域的放大图：

将重建后的图像与原始卷积核进行比较，可以发现它们非常相似；实际上，在重建图像上应用介于0.5到0.68之间的任何截止值都会恢复原始卷积核。 重建图像中峰值周围的微弱波纹是由于舍入和边缘效应引起的噪声。

我不确定造成重建图像中出现十字形伪影的原因是什么（但我肯定有更多这方面经验的人可以告诉您），但是我的猜测是它与图像边缘有关。 当我在GIMP中卷积原始图像时，我要求它通过扩展图像来处理边缘（实际上是复制最外层的像素），而FFT去卷积则假设图像边缘环绕。 这可能会在重建图像中引入沿x和y轴的虚假相关性。

- Ilmari Karonen

---

谢谢。虽然不是精确的内核，但对我来说结果已经足够好了。 - Martin Perry

我正在尝试将代码重写为C语言。rfft2是实数输入的FFT，也称为DFT或DCT，并返回实数输出？而这一行：kernel_f = blur_f / orig_f..它是每个元素的除法还是矩阵除法（乘以逆矩阵）？ - Martin Perry

我相信rfft2的输出是复数，但它省略了从将通用复FFT应用于实际输入中获得的对称一半输出。它不是DCT。无论如何，任何FFT例程都应该可以做到。是的，numpy数组上的/只是逐元素除法（即乘以倒数）。如果需要，我认为我可以稍后提供带有GSL的C代码（但现在不行，我还没有清醒到那个程度）。 - Ilmari Karonen

C代码会很好……我正在使用fftw和实数到实数FFT，其中标记为DCT。我尝试运行实数到复数FFT。逐元素除法是通过求逆来进行乘法吗？我认为，逐元素意味着res[0][0] = m1[0][0] / m2[0][0]等（但是用逆来做并不那么简单）。 - Martin Perry

这是复数除法，就像傅里叶空间中的普通卷积需要进行复数乘法一样。附：我没有使用FFTW，但我相信你应该使用的计划是fftw_plan_dft_r2c_2d()。 - Ilmari Karonen

Answer 2

这是一个典型的反卷积问题。你所称之为卷积矩阵通常被称为“核”。卷积操作通常用星号“*”表示（不要与乘法混淆！）。使用这个符号：

结果 = 源 * 核

上面使用FFT的答案是正确的，但在存在噪声的情况下，你不能真正使用基于FFT的反卷积。正确的方法是使用Richardson-Lucy反卷积（参见https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution）。

实现起来非常简单。此答案还提供了一个示例Matlab实现：Would Richardson–Lucy deconvolution work for recovering the latent kernel?

- yvoronen

Answer 3

我已经使用fftw3将@Ilmari Karonen answer 重写成C/C++，供需要的人参考：

循环移位函数

template<class ty>
void circshift(ty *out, const ty *in, int xdim, int ydim, int xshift, int yshift)
{
  for (int i =0; i < xdim; i++) 
  {
    int ii = (i + xshift) % xdim;
    for (int j = 0; j < ydim; j++) 
    {
      int jj = (j + yshift) % ydim;
      out[ii * ydim + jj] = in[i * ydim + j];
    }
  }
}

现在的主要代码

int width = 256;
int height = 256;

int index = 0;

MyStringAnsi imageName1 = "C://ka4ag.png";    
MyStringAnsi imageName2 = "C://KyPu2.png";

double * in1 = new double[width * height];
fftw_complex * out1 = new fftw_complex[width * height]; 

double * in2 = new double[width * height];
fftw_complex * out2 = new fftw_complex[width * height]; 

MyUtils::MyImage * im1 = MyUtils::MyImage::Load(imageName1, MyUtils::MyImage::PNG);
MyUtils::MyImage * im2 = MyUtils::MyImage::Load(imageName2, MyUtils::MyImage::PNG);

for (int i = 0; i < width * height; i++)
{
    in1[i] = ((im1->Get(i).r / (255.0 * 0.5)) - 1.0);
    in2[i] = ((im2->Get(i).r / (255.0 * 0.5)) - 1.0);
}


fftw_plan dft_plan1 = fftw_plan_dft_r2c_2d(width, height, in1, out1, FFTW_ESTIMATE);    
fftw_execute(dft_plan1);
fftw_destroy_plan(dft_plan1);

fftw_plan dft_plan2 = fftw_plan_dft_r2c_2d(width, height, in2, out2, FFTW_ESTIMATE);    
fftw_execute(dft_plan2);
fftw_destroy_plan(dft_plan2);

fftw_complex * kernel = new fftw_complex[width * height];   

for (int i = 0; i < width * height; i++)
{
    std::complex<double> c1(out1[i][0], out1[i][1]);
    std::complex<double> c2(out2[i][0], out2[i][1]);

    std::complex<double> div = c2 / c1;

    kernel[i][0] = div.real();
    kernel[i][1] = div.imag();
}

double * kernelOut = new double[width * height];

fftw_plan dft_planOut = fftw_plan_dft_c2r_2d(width, height, kernel, kernelOut, FFTW_ESTIMATE);
fftw_execute(dft_planOut);
fftw_destroy_plan(dft_planOut);

double * kernelShift = new double[width * height];

circshift(kernelShift, kernelOut, width, height, (width/2), (height/2));

double maxKernel = kernelShift[0];
for (int i = 0; i < width * height; i++)
{
    if (maxKernel < kernelShift[i]) maxKernel = kernelShift[i]; 
}

for (int i = 0; i < width * height; i++)
{
    kernelShift[i] /= maxKernel; 
}

uint8 * res = new uint8[width * height];
for (int i = 0; i < width * height; i++)
{                   
   res[i] = static_cast<uint8>((kernelShift[i]+ 1.0) * (255.0 * 0.5));
}

//now in res is similar result as in @Ilmari Karonen, but shifted by +128

代码没有内存管理，因此您必须清理您的内存！

- Martin Perry

Answer 4

如果已知卷积矩阵大小的上限，则可以生成估计值。 如果它是N，则选择N * N个点，并尝试根据源和目标数据解决线性方程组，以获得卷积系数。 鉴于颜色分量的四舍五入，该系统无法解决，但是通过线性规划，您将能够通过对这些系数进行小的修改来最小化总偏差。

- Vesper

Answer 5

你可以尝试使用源图像作为卷积核来执行反卷积。但结果可能是不可预测的 - 反卷积是非常不稳定的过程，由于噪声、边缘效应、舍入误差等原因。

- MBo