如果在FFT中有很多0，如何寻找卷积核?

Question

如果在FFT中有很多0，如何寻找卷积核?

matlabimage-processingfftconvolutiondeconvolution

3

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

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

结果没有意义。我发现

包含许多NaN值，fftblur和fftorig都包含许多0值。我需要做些什么吗？例如使用fftshift吗？

编辑：为了更容易理解，我现在使用来自http://matlabgeeks.com/tips-tutorials/how-to-blur-an-image-with-a-fourier-transform-in-matlab-part-i/的图像。

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

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

这是结果：

https://imgur.com/a/b7uvj

很明显，这与链接顶部提到的高斯模糊不匹配。

- user5739619

你在这里发布的代码中使用了 /，但你想要使用 ./。然而，这不应该阻止 NaN 值出现在两个操作数都为 0 的情况下。 - Cris Luengo

在频域中进行除法是一个技巧，在你的信号/系统课程中教授后很少有用，你已经看到了其中一个原因：数值不稳定性。（你还混淆了线性卷积和循环卷积，因为没有对FFT进行零填充。）当Matlab图像处理工具箱有一系列去模糊函数时，你知道这是一个复杂的主题。 - Ahmed Fasih

1

对你刚刚发布的链接做一个评论：它说“内核可以放置在图像的任何位置；这并不重要。”这是错误的。卷积核需要居中于原点，即左上角像素。如果你没有做到这一点，在频域中的卷积核将不是纯实数，并且会改变图像频率成分的相位。你可以在他们的第一个结果中看到这一点，输出图像被移位了！ - Cris Luengo

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Cris Luengo · Accepted Answer

这就是维纳滤波器发挥作用的地方。通常用于去卷积 - 从滤波图像和卷积核估计原始未经滤波的图像。然而，由于卷积的可交换性，去卷积正是您试图解决的完全相同的问题。也就是说，如果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取决于输入图像，但有技巧可以估计它。此外，像这样的简单滤波器永远无法完全消除我在此处施加的强烈模糊。需要更先进的迭代方法才能获得更好的去卷积效果。

估计核

让我们反过来，从original和in估计psf。直接进行除法（相当于Wiener滤波器，K = 0）是可能的，但输出非常嘈杂。原始图像具有非常低的频域值的地方，估计值很差。选择正确的K可以更好地估计核。如果K太大，则结果仍然是一个很差的近似。

% 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等内容，得到了一个相当不错的结果，显示出略微偏移的高斯核。

然后我对第二个结果（填充后）进行了相同的操作，得到了一个更差的结果，但仍然非常容易辨认出是高斯核。