找到合适的陷波滤波器以从图像中去除图案
Find proper notch filter to remove pattern from image
我想在图像上应用陷波滤波器,它可以抑制图像中的图案,但尽可能保持图像的其余部分完好无损。
我执行以下步骤:
I = imread('...img....');
ft = fftshift(fft2(I));
[m,n] = size(ft);
filt = ones(m,n);
%filt(......) = 0; % my problem is here
ft = ft .* filt;
ifft_ = ifft2(ifftshift( ft));
所以我不知道要将什么设置为零才能获得正确的结果。
如果您查看图像的 fft,您可以清楚地看到导致图像图案的强频率。
您需要创建一个陷波滤波器,将这些高峰值周围的区域清零。我尝试使用高斯陷波滤波器进行此操作,结果频谱看起来像这样。
ifft图像(对比度增强)结果是
这是一些用于构建和应用过滤器的 MATLAB 代码
I = imread('YmW3f.png');
ft = fftshift(fft2(I));
[m,n] = size(ft);
% define some functions
norm_img = @(img) (img - min(img(:))) / (max(img(:)) - min(img(:)));
show_spec = @(img) imshow(norm_img(log(abs(img)-min(abs(img(:)))+1.0001)));
gNotch = @(v,mu,cov) 1-exp(-0.5*sum((bsxfun(@minus,v,mu).*(cov\bsxfun(@minus,v,mu)))));
% show spectrum before
figure();
show_spec(ft);
% by inspection
cx = 129;
cy = 129;
% distance of noise from center
wx1 = 149.5-129;
wx2 = 165.5-129;
wy = 157.5-129;
% create notch filter
filt = ones(m,n);
% use gaussian notch with standard deviation of 5
sigma = 5;
[y,x] = meshgrid(1:n, 1:m);
X = [y(:) x(:)].';
filt = filt .* reshape(gNotch(X,[cx+wx1;cy+wy],eye(2)*sigma^2),[m,n]);
filt = filt .* reshape(gNotch(X,[cx+wx2;cy+wy],eye(2)*sigma^2),[m,n]);
filt = filt .* reshape(gNotch(X,[cx-wx1;cy-wy],eye(2)*sigma^2),[m,n]);
filt = filt .* reshape(gNotch(X,[cx-wx2;cy-wy],eye(2)*sigma^2),[m,n]);
% apply filter
ft = ft .* filt;
% show spectrum after
figure();
show_spec(ft);
% compute inverse
ifft_ = ifft2(ifftshift( ft));
img_res = histeq(norm_img(ifft_));
figure();
imshow(img_res);
编辑:已交换 meshgrid
的参数,原因由 Todd Gillette 指出。
jodag 的回答对我来说效果很好,除了在 meshgrid 命令中交换了 m 和 n。应该是
[y,x] = meshgrid(1:n, 1:m);
它在示例中有效,因为图像是方形的,但对于矩形图像,它无法正常工作。
[我本来想发表评论,但我还没有名气。]
我想在图像上应用陷波滤波器,它可以抑制图像中的图案,但尽可能保持图像的其余部分完好无损。 我执行以下步骤:
I = imread('...img....');
ft = fftshift(fft2(I));
[m,n] = size(ft);
filt = ones(m,n);
%filt(......) = 0; % my problem is here
ft = ft .* filt;
ifft_ = ifft2(ifftshift( ft));
所以我不知道要将什么设置为零才能获得正确的结果。
如果您查看图像的 fft,您可以清楚地看到导致图像图案的强频率。
您需要创建一个陷波滤波器,将这些高峰值周围的区域清零。我尝试使用高斯陷波滤波器进行此操作,结果频谱看起来像这样。
ifft图像(对比度增强)结果是
这是一些用于构建和应用过滤器的 MATLAB 代码
I = imread('YmW3f.png');
ft = fftshift(fft2(I));
[m,n] = size(ft);
% define some functions
norm_img = @(img) (img - min(img(:))) / (max(img(:)) - min(img(:)));
show_spec = @(img) imshow(norm_img(log(abs(img)-min(abs(img(:)))+1.0001)));
gNotch = @(v,mu,cov) 1-exp(-0.5*sum((bsxfun(@minus,v,mu).*(cov\bsxfun(@minus,v,mu)))));
% show spectrum before
figure();
show_spec(ft);
% by inspection
cx = 129;
cy = 129;
% distance of noise from center
wx1 = 149.5-129;
wx2 = 165.5-129;
wy = 157.5-129;
% create notch filter
filt = ones(m,n);
% use gaussian notch with standard deviation of 5
sigma = 5;
[y,x] = meshgrid(1:n, 1:m);
X = [y(:) x(:)].';
filt = filt .* reshape(gNotch(X,[cx+wx1;cy+wy],eye(2)*sigma^2),[m,n]);
filt = filt .* reshape(gNotch(X,[cx+wx2;cy+wy],eye(2)*sigma^2),[m,n]);
filt = filt .* reshape(gNotch(X,[cx-wx1;cy-wy],eye(2)*sigma^2),[m,n]);
filt = filt .* reshape(gNotch(X,[cx-wx2;cy-wy],eye(2)*sigma^2),[m,n]);
% apply filter
ft = ft .* filt;
% show spectrum after
figure();
show_spec(ft);
% compute inverse
ifft_ = ifft2(ifftshift( ft));
img_res = histeq(norm_img(ifft_));
figure();
imshow(img_res);
编辑:已交换 meshgrid
的参数,原因由 Todd Gillette 指出。
jodag 的回答对我来说效果很好,除了在 meshgrid 命令中交换了 m 和 n。应该是
[y,x] = meshgrid(1:n, 1:m);
它在示例中有效,因为图像是方形的,但对于矩形图像,它无法正常工作。
[我本来想发表评论,但我还没有名气。]