6

I want to apply notch filter on an image where it suppresses the pattern in the image but leave the rest of the image as intact as possible. i do the following steps : 

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));

so i don't know what exactly to set to zero to get the proper result. enter image description here

rayryeng
  • 102,964
  • 22
  • 184
  • 193
Yazan Houshieh
  • 113
  • 1
  • 2
  • 11

2 Answers2

10

If you take a look at the fft of the image you can clearly see the strong frequencies that are causing the pattern in the image.

enter image description here

You need to create a notch filter which zeros out the region around those high peaks. I tried using Gaussian notch filters for this operation and the resulting spectrum looked something like this.

enter image description here

The ifft image (with contrast enhanced) turns out to be

enter image description here

Here's some MATLAB code used to build and apply the filter

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);

Edit: Swapped parameters for meshgrid for reason indicated by Todd Gillette.

jodag
  • 19,885
  • 5
  • 47
  • 66
  • Yes, sharing your code would be very beneficial... especially since we don't know what parameters you chose for your Gaussian functions to smooth out the high peaks in the spectrum. – rayryeng Mar 25 '15 at 15:14
3

jodag's answer works well for me, except for a swapped m and n in the meshgrid command. It should be

    [y,x] = meshgrid(1:n, 1:m);

It worked in the example because the image is square, but with a rectangular image it does not work properly.

[I would have preferred to comment, but I don't yet have the reputation.]

Todd Gillette
  • 150
  • 1
  • 11