15

我有源图像和结果图像。我知道，已经在源上使用了一些卷积矩阵来获得结果。可以计算这个卷积矩阵吗？或者至少不是确切的一个，但非常相似。

4

5 回答 5

24

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

要了解为什么会这样，请注意空间域中的卷积对应于频域中的乘法，因此反卷积同样对应于频域中的除法。在普通的反卷积中，将卷积图像的 FFT 除以内核的 FFT 以恢复原始图像。但是，由于卷积（与乘法一样）是一种交换运算，因此内核和源的角色可以任意互换：内核对源进行卷积与源对内核进行卷积完全相同。

然而，正如其他答案所指出的那样，这不太可能产生内核的精确重建，原因与普通反卷积通常不会完全重建原始图像相同：舍入和裁剪会在过程中引入噪声，这是可能的让卷积完全消除一些频率（通过将它们与零相乘），在这种情况下，这些频率无法重建。

也就是说，如果您的原始内核大小有限（支持），则重建后的内核通常应该在原点周围有一个尖峰，近似于原始内核，并被低级噪声包围。即使您不知道原始内核的确切大小，提取这个峰值并丢弃其余的重建应该也不难。

例子：

这是灰度的Lenna 测试图像，缩小到 256×256 像素，并在 GIMP 中与 5×5 内核进行卷积：

原来的* 核心→结果

我将使用 Python 和 numpy/scipy 进行反卷积：

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 轴引入虚假相关性。

于 2013-05-14T12:05:50.553 回答
2

这是一个经典的反卷积问题。您所说的卷积矩阵通常称为“内核”。卷积运算通常用星号“*”表示（不要与乘法混淆！）。使用这个符号

Result = Source * Kernel

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

实施起来非常简单。这个答案还提供了一个 Matlab 实现示例：Richardson-Lucy 反卷积是否可以用于恢复潜在内核？

于 2015-12-07T22:03:23.170 回答
1

我已经使用 fftw3 为某人重写了@Ilmari Karonen对 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

代码没有内存管理，所以你必须清理你的内存！

于 2013-05-16T12:15:48.937 回答
1

好吧，如果卷积矩阵大小的上限已知，则可以生成估计值。如果是 N，则选择 N*N 个点并尝试根据源和目标的数据求解一个针对卷积系数的线性方程组。考虑到颜色分量的舍入，系统将无法解决，但通过线性规划，您将能够通过对这些系数进行小的更改来最小化与预期值的总偏移量。

于 2013-05-14T11:20:02.750 回答
1

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

于 2013-05-14T11:48:26.330 回答