filtering using FFT in images

Hi, I am experimenting with masking out areas on a FFT of an image to see the effect on the processed image:

After reading, I have the following.

I=I-mean(I(:));
f = fftshift(fft2(I));
fabs=abs(f);
figure
subplot(1,3,1)
imshow(fabs,[])

After viewing the FFT, I notice there are some white spots, so to attempt to remove these, I have defined a threshold of 90% of the max. My aim is to remove these from the fabs image and then inverse fft and then view this processed image.

thresh=0.9*max(fabs(:))

Im not too sure how to apply the the threshold to the frequency domain image and then iFFT back.

Answer 1

[Image Analyst]

Try this:

% Find pixels that are brighter than the threshold.
mask = fabs > thresh; 
% Erase those from the image
fabs(mask) = 0;
% Shift back and inverse fft
filteredImage = ifft2(fftshift(fabs)) + mean2(I);
imshow(filteredImage, []);

[Jason]

Thanks IA, your code makes sense, but Im not seeing what I expect. For example, when I put the threshold "thresh" at the max value (thresh is value of the FFT image above which I want to remove), I should get the original image back. Also, why do you have to subtract and then add the mean value?

%FFT and remove noise
I=getimage(handles.axes3);
figure
set(gcf, 'units','normalized','outerposition',[0 0 0.7 0.7]);

%View the Raw image
subplot(1,4,1)
[high,low]=Autoscaleimage(handles,I,2);
imshow(I,[low,high])
title('Raw Image')

%Perform the FFT
I=I-mean(I(:));   %Not sure why we need to subtract the mean???
f = fftshift(fft2(I));
fabs=abs(f);

%View the Frequency space image
subplot(1,4,2)
[high,low]=Autoscaleimage(handles,fabs,2);
imshow(fabs.^1.0,[low,high])
title('FFT')

%Spacially threshold the FFt image
%Remove brightest pixels

thresh=0.95*max(fabs(:));      %Define threshold
mask = fabs > thresh;          %Find pixels that are brighter than the threshold.
fabs(mask) = 0;                %Erase those from the image
filteredImage = abs(ifft2(fftshift(fabs)) + mean2(I)); % Shift back and inverse fft

%View the thresholded FFT Image
subplot(1,4,3)
[high,low]=Autoscaleimage(handles,fabs,2);
imshow(fabs,[low,high])
title('FFT with filtering')

%View the real space image to see effect
subplot(1,4,4)
[high,low]=Autoscaleimage(handles,filteredImage,2);
imshow(filteredImage,[low,high])
title('Inverse FFT')

[Image Analyst]

We'd need the original image by itself to run your code.

[Jason]

Here you go, its a tiff so I had to zip it up.

[Image Analyst]

Your filtering didn't even really filter it. If the spectrum can't be thresholded appropriately then you may have to run adapthisteq() on it first then threshold. Then get a mask and zero out and inverse transform

fabs2 = adapthisteq(...........
mask = fabs2 > someValue
filteredSpectrum  = fabs;
filteredSpectrum(mask) = 0;
% Then fftshift and ifft...

Use my interactive thresholding app in http://www.mathworks.com/matlabcentral/fileexchange/?term=authorid%3A31862

[Jason]

I think we were applying the filter incorrectly, the following works:-

%Perform the FFT
I=I-mean(I(:));   %Not sure why we need to subtract the mean???
f = fftshift(fft2(I));
fabs=abs(f);

%use log to plot
fLog = log(1 + abs(f));
% filter by a range based on fLog
filter = (fLog > .7*max(fLog(:)) ); 

%apply filter
B = abs(ifft2(f.*filter));

figure
colormap(gray)
subplot(2,2,1),imagesc(I); title('Original Image')
subplot(2,2,2),imagesc(fLog); title('Fourier Image')
subplot(2,2,3),imagesc(filter); title('Zeroed Fourier Image')
subplot(2,2,4),imagesc(B); title('Cleaned Image')

[Image Analyst]

Taking the log will work for certain shapes but is not good in general. I'm not sure who "we" are but I know that adapthisteq()' like I recommended' will flatten the background all over and doesn't depend on the "background" following an exponential shape, like taking the log does. I think you didn't try adapthisteq() or never figured out the best parameters to use with it. Anyway, glad it's working for you with this method.

[Jason]

Sorry, I mean "I" applied the filter incorrectly. Now i have it working I have tried to play around with adapthisteq(). Im not sure i understand properly its use, as the FFT image now is just a single graylevel

I tried replacing the

fLog=log(1+abs(f)) with 

fLog=adapthisteq(fabs);

[Jason]

Is there away to do what I have done but to leave the central region of the FFT image i.e. DC component

(i.e. don't mask this out)?

[Image Analyst]

It looks like you zeroed out the wrong things if you want to get rid of the periodic patterns.

[Jason]

hmmm. Im obviously not getting this!

[Image Analyst]

What do you want to do? The spikes and bars in the spectrum represent periodic structures in the spatial domain. If you get rid of those spikes, you should reduce the periodic patterns in the spatial domain.

They're inversely related. Fast changing spatial bars give rise to widely spaced spikes in the Fourier spectrum and vice versa. So the 4 big spikes you see are related to the sine wave pattern on the white rectangles. The other stuff is a 2D sinc due to there being a couple of bright squares in the image. The big spike at 0, the "DC spike" is due to the fact that the image has a non-zero mean gray level.

[Jason]

Thanks for your explanation. One thing why is there a dc offset as i subtracted off the mean originally. Are there any tricks to reject this?