如果FFT有很多0,那么查找卷积内核?

编程入门 行业动态 更新时间:2024-10-15 18:30:56
本文介绍了如果FFT有很多0,那么查找卷积内核?的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧! 问题描述

我知道 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 取决于输入图像,但有一些技巧可以估算它。此外,像这样的简单滤镜永远无法完美地消除像我在这里所说的刺耳的模糊。需要更高级的迭代方法才能获得更好的反卷积。

估算内核

让我们来看看,并从原始和中估算 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: matlabgeeks/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:

imgur/a/b7uvj

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,那么查找卷积内核?

本文发布于:2023-11-30 03:49:57,感谢您对本站的认可!
本文链接:https://www.elefans.com/category/jswz/34/1648620.html
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系,我们将在24小时内删除。
本文标签:卷积   有很多   内核   FFT

发布评论

评论列表 (有 0 条评论)
草根站长

>www.elefans.com

编程频道|电子爱好者 - 技术资讯及电子产品介绍!