如果 FFT 有很多 0,则寻找卷积核?

Finding for convolution kernel if many 0's for FFT?

我知道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 值,并且 fftblurfftorig 都包含许多 0 值。我需要对此做些什么吗?例如使用 fftshift?

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

我决定使用 link 计算 origimageblurimageunpad 的内核:

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

结果如下:

https://imgur.com/a/b7uvj

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

这就是维纳滤波器派上用场的地方。它通常用于反卷积——从过滤后的图像和卷​​积核中估计原始的、未过滤的图像。但是,由于卷积的交换性,反卷积 与您要解决的问题 完全相同。也就是说,如果g = f * h,那么可以估计f 来自 gh(反卷积),或者您可以估计 h 来自 gf,完全相同。

解卷积

Wiener filter Wikipedia page 在数学上很重,但以一种简单的方式,您寻找内核具有较小值(与图像中的噪声相比)的频率,并且您不应用在这些频率上划分。这称为正则化,您可以在其中施加一些约束以避免噪声主导结果。

这是 MATLAB 代码,对模糊输入图像使用 in,对 spatial-domain 过滤器使用 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 取决于输入图像,尽管有一些技巧可以估计它。此外,像这样的简单滤镜永远无法完美消除我在此处设置的刺眼模糊效果。需要更高级的迭代方法才能获得更好的反卷积。

估计内核

让我们反过来,从 originalin 估计 psf。可以直接进行除法(相当于 K=0 的 Wiener 滤波器),但输出噪声很大。如果原始图像的 frequency-domain 值非常低,您将得到很差的估计值。选择正确的 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 和所有内容,得到了一个很好的结果,显示了一个稍微偏移的高斯内核。

然后我对第二个结果重复了相同的操作(填充后)并得到了更差的结果,但仍然非常容易识别为高斯分布。