如果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 值， fftblur 和 fftorig 都包含许多0的值。我需要对此做些什么吗？例如使用 fftshift ？
编辑：
使这更容易要理解，我现在正在使用以下图像：
如您所见，选择正确的 K 非常重要。如果它太低，则没有足够的正则化。如果它太高那么人工制品太多（反卷积不足）。此参数 K 取决于输入图像，但有一些技巧可以估算它。此外，像这样的简单滤镜永远无法完美地消除像我在这里所说的刺耳的模糊。需要更高级的迭代方法才能获得更好的反卷积。
估算内核
让我们来看看，并从原始和中<$ code>估算 psf 。可以直接进行除法（相当于Wiener滤波器， K = 0 ），但输出非常嘈杂。如果原始图像具有非常低的频域值，您将得到较差的估计值。选择正确的 K 可以更好地估计内核。如果 K 太大，结果再次很差。
％按分部直接估算
psf1 = fftshift（ifft2（fft2（in）./ fft2（original）））;
％Wiener过滤器方法
K = 1e-7;
H = fft2（原始）;
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 和所有东西，得到了一个非常好的结果，显示了一个稍微偏移的高斯内核。
然后我重复了第二个结果（填充后）并得到了更糟糕的结果，但仍然非常容易识别为高斯。
I know that original_image * filter = blur_image, where * is the convolution. Thus, filter = ifft(fft(blur)/fft(original))
I have an original image, the known filter, and the known blurred image. I tried the following code. I just want to compare the computed filter using fft and ifft and compare it with the known filter. I tried in Matlab:
orig = imread("orig.png")
blur = imread("blur.png")
fftorig = fft(orig)
fftblur = fft(blur)
div = fftblur/fftorig
conv = ifft(div)
The result doesn't make sense. I see that div contains many NaN values, and fftblur and fftorig both contain many values of 0. Do I need to do something to this? such as use fftshift?
EDIT: To make this easier to understand, I'm now using the images from: http://matlabgeeks.com/tips-tutorials/how-to-blur-an-image-with-a-fourier-transform-in-matlab-part-i/
I decided to compute the kernel of the origimage and blurimageunpad from that link using:
kernelc = real(ifft2(fft2(origimage)./fft2(blurimageunpad));
imagesc(kernelc)
colormap gray
Here's the result:
Which clearly doesn't match the gaussian blur mentioned at towards the top of that link
This is where the Wiener filter comes in handy. It is typically applied for deconvolution -- estimating the original, unfiltered image from the filtered image and the convolution kernel. However, because of the commutativity of the convolution, deconvolution is exactly the same problem you are trying to solve. That is, if g = f * h, then you can estimate f from g and h (deconvolution), or you can estimate h from g and f, in exactly the same manner.
Deconvolution
The Wiener filter Wikipedia page is heavy on the math, but in a simplistic way, you look for frequencies where your kernel has small values (compared to the noise in your image), and you don't apply the division at those frequencies. This is called regularization, where you impose some constraints to avoid noise dominating the results.
This is the MATLAB code, using in for the blurred input image, and psf for the spatial-domain filter:
% 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)));
Here is the image in and the output of the Wiener filter for three different values of K:
As you can see, selecting the right K is very important. If it is too low, there is not sufficient regularization. If it is too high then there are too many artefacts (insufficient deconvolution). This parameter K depends on the input image, though there are tricks to estimate it. Also, a simple filter like this can never perfectly undo a harsh blur like I put on here. More advanced iterative approaches are necessary to obtain a better deconvolution.
Estimating the kernel
Let's turn this around, and estimate psf from original and in. It is possible to do the division directly (equivalent to the Wiener filter with K=0), but the output is quite noisy. Where the original image has very low frequency-domain values, you'll get poor estimates. Picking the right K leads to a better estimate of the kernel. With a K that is too large, the result is again a poor approximation.
% 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))));
(zoom in to observe the noise)
Edit
I downloaded the images from the website you linked. I used the first result, without padding, and cropped the frames from the images as best as I could, to leave only the data and exactly the data:
original = imread('original.png');
original = original(33:374,120:460,1);
in = imread('blur_nopad.png');
in = in(33:374,120:460,1);
Then I used the code exactly as above, with the same K and everything, and got a pretty good result, showing a slightly shifted Gaussian kernel.
I then repeated the same with the second result (after padding) and got a worse result, but still very recognizable as a Gaussian.
这篇关于如果FFT有很多0，那么查找卷积内核？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
