如果 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
?
编辑:
为了更容易理解,我现在使用的图片来自:http://matlabgeeks.com/tips-tutorials/how-to-blur-an-image-with-a-fourier-transform-in-matlab-part-i/
我决定使用 link 计算 origimage
和 blurimageunpad
的内核:
kernelc = real(ifft2(fft2(origimage)./fft2(blurimageunpad));
imagesc(kernelc)
colormap gray
结果如下:
这显然与顶部提到的高斯模糊不匹配 link
这就是维纳滤波器派上用场的地方。它通常用于反卷积——从过滤后的图像和卷积核中估计原始的、未过滤的图像。但是,由于卷积的交换性,反卷积 与您要解决的问题 完全相同。也就是说,如果g = f * h,那么可以估计f 来自 g 和 h(反卷积),或者您可以估计 h 来自 g 和 f,完全相同。
解卷积
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
取决于输入图像,尽管有一些技巧可以估计它。此外,像这样的简单滤镜永远无法完美消除我在此处设置的刺眼模糊效果。需要更高级的迭代方法才能获得更好的反卷积。
估计内核
让我们反过来,从 original
和 in
估计 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
和所有内容,得到了一个很好的结果,显示了一个稍微偏移的高斯内核。
然后我对第二个结果重复了相同的操作(填充后)并得到了更差的结果,但仍然非常容易识别为高斯分布。
我知道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-part-i/
我决定使用 link 计算 origimage
和 blurimageunpad
的内核:
kernelc = real(ifft2(fft2(origimage)./fft2(blurimageunpad));
imagesc(kernelc)
colormap gray
结果如下:
这显然与顶部提到的高斯模糊不匹配 link
这就是维纳滤波器派上用场的地方。它通常用于反卷积——从过滤后的图像和卷积核中估计原始的、未过滤的图像。但是,由于卷积的交换性,反卷积 与您要解决的问题 完全相同。也就是说,如果g = f * h,那么可以估计f 来自 g 和 h(反卷积),或者您可以估计 h 来自 g 和 f,完全相同。
解卷积
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
取决于输入图像,尽管有一些技巧可以估计它。此外,像这样的简单滤镜永远无法完美消除我在此处设置的刺眼模糊效果。需要更高级的迭代方法才能获得更好的反卷积。
估计内核
让我们反过来,从 original
和 in
估计 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
和所有内容,得到了一个很好的结果,显示了一个稍微偏移的高斯内核。
然后我对第二个结果重复了相同的操作(填充后)并得到了更差的结果,但仍然非常容易识别为高斯分布。