如果 FFT 中有很多 0，则寻找卷积核？

Question

我知道original_image * filter = blur_image，*卷积在哪里。因此，filter = ifft(fft(blur)/fft(original))

我有一个原始图像、已知的滤镜和已知的模糊图像。我尝试了以下代码。我只想比较使用 fft 和 ifft 计算的滤波器，并将其与已知的滤波器进行比较。我在Matlab中尝试过：

orig = imread("orig.png")
blur = imread("blur.png")
fftorig = fft(orig)
fftblur = fft(blur)
div = fftblur/fftorig
conv = ifft(div)

结果没有意义。我看到div包含很多NaN值，并且fftblur和fftorig都包含很多 0 值。我需要对此做些什么吗？比如使用fftshift？

编辑：为了使这一点更容易理解，我现在使用以下图像： http: //matlabgeeks.com/tips-tutorials/how-to-blur-an-image-with-a-fourier-transform-in- matlab-第一部分/

我决定使用以下链接计算origimage和的内核：blurimageunpad

kernelc = real(ifft2(fft2(origimage)./fft2(blurimageunpad));
imagesc(kernelc)
colormap gray

结果如下：

https://i.stack.imgur.com/kJavg.jpg

这显然与该链接顶部提到的高斯模糊不匹配

Answer 1

这就是维纳滤波器派上用场的地方。它通常用于反卷积——根据滤波后的图像和卷​​积核估计原始的、未滤波的图像。然而，由于卷积的交换性，反卷积与您要解决的问题完全相同。也就是说，如果g = f * h ，那么您可以根据g和h估计f （反卷积），或者您可以以完全相同的方式根据g和f估计h 。

反卷积

维纳滤波器维基百科页面的数学内容较多，但以一种简单化的方式，您会寻找内核具有较小值的频率（与图像中的噪声相比），并且不会在这些频率处应用除法。这称为正则化，您可以施加一些约束以避免噪声主导结果。

这是 MATLAB 代码，用于in模糊输入图像和psf空间域滤波器：

% Create a PSF and a blurred image:
original = imread('cameraman.tif');
sz = size(original);
psf = zeros(sz);
psf(fix(sz(1)/2)+(-5:5),fix(sz(1)/2)+(-10:10)) = 1;
psf = psf/sum(psf(:));
% in = conv2(original,psf,'same'); % gives boundary issues, less of a problem for smaller psf
in = uint8(ifft2(fft2(original).*fft2(ifftshift(psf))));

% Input parameter:
K = 1e-2; 

% Compute frequency-domain PSF:
H = fft2(ifftshift(psf));

% Compute the Wiener filter:
cH = conj(H);
HcH = H .* cH;
K = K * max(max(abs(HcH)));
w = cH ./ (HcH + K);

% Apply the Wiener filter in the Frequency domain:
out = real(ifft2(w .* fft2(in)));

以下是 的三个不同值的维纳滤波器的图像in和输出K：

正如您所看到的，选择正确的K非常重要。如果太低，则正则化不充分。如果它太高，则存在太多伪影（反卷积不足）。该参数K取决于输入图像，尽管有一些技巧可以估计它。另外，像这样的简单滤镜永远无法完美消除像我在这里设置的严酷模糊。为了获得更好的反卷积，需要更先进的迭代方法。

估计内核

让我们反过来，psf从original和进行估计in。可以直接进行除法（相当于 的维纳滤波器K=0），但输出的噪声很大。如果原始图像的频域值非常低，您将得到较差的估计。选择正确的K可以更好地估计内核。如果 aK太大，结果又是一个很差的近似值。

% Direct estimation by division
psf1 = fftshift(ifft2(fft2(in)./fft2(original)));

% Wiener filter approach
K = 1e-7;
H = fft2(original);
cH = conj(H);
HcH = H .* cH;
K = K * max(max(abs(HcH)));
w = cH ./ (HcH + K);
psf2 = fftshift(real(ifft2(w .* fft2(in))));

（放大观察噪音）

编辑

我从您链接的网站下载了图像。我使用了第一个结果，没有填充，并尽可能地从图像中裁剪帧，只留下数据和准确的数据：

original = imread('original.png');
original = original(33:374,120:460,1);

in = imread('blur_nopad.png');
in = in(33:374,120:460,1);

然后我使用了与上面完全相同的代码，具有相同K的所有内容，并得到了相当好的结果，显示了稍微偏移的高斯内核。

然后，我对第二个结果（填充后）重复相同的操作，得到了更糟糕的结果，但仍然很容易被识别为高斯分布。