How can I remove the additive periodical noise from the image?

I got an image (Project1.png) that had been corrupted by periodical noise A= 20 * sin(0.1*pi*x+0.1*pi*y). Try to remove the noise by apply the DFT on the image and multiple with the design filter (FilterImage.jpg) that I been created then inverse back to get the filtered image. But its seem like the final image not remove the noise at all. Any idea how to remove it or suggestion on how to do it? The design filter I created was basically based on the DFT of the noise (F=fft2((A)); Fs=(fftshift(F));) Thanks

clear all;

close all;

clc;

[filename, pathname]=uigetfile('*.*','Select Grey Scale Image');

filewithpath=strcat(pathname,filename);

img=imread(filewithpath);

[row,col]=size(img)

F=fft2(double(img)); %DFT of Image

Fs=fftshift(F); %Shifting Spectrum to Centre

% Getting H(u,v) - Filter

H=imread('FilterImage.jpg');

H=double(H/255);

Fsf=Fs.*H;

%Inverse DFT

fimg=ifft2(fftshift(Fsf));

imgr=uint8(real(fimg));

subplot(221);

imshow(img);

subplot(222);

imshow(log(1+abs(Fs)),[]);

subplot(223);

imshow(imgr, []);

subplot(224);

imshow(log(1+abs(Fsf)),[]);

Answer 1

[Bjorn Gustavsson]

Your Filterimage is not the notch-filter you need, that is only the 1-D amplitude (possibly power) of the 1-D fft of the filter you need. Have a look at chapter 5.7 Interactive Restoration (I've got an old edition), there you should find the information you need to produce a 2-D notch-filter removing the relevant "interference" peaks from the 2-D fft of your image.

HTH

[AP]

Thank you so much. I need to have a look back again. But wondering what book that you mean for this? 'chapter 5.7 Interactive Restoration'.

[Bjorn Gustavsson]

I instantly assumed that the project1.png image was from the front-page of Digital Image Processing by Gonzalez and Woods (and Eddins depending on variant). In my edition this type of notch-filter for supression of harmonic interference signal is described in detail in section 5.7. Since you are trying to learn this technique I thought it would be best to direct you to one of the best descriptions of that technique.

The important steps are:

1, to realise that the interference-signal will have a distinctive 2-D pattern that has a few spikes in the 2-D Fourier-domain.

2, that it is reasonably straightforwar to set those Fourier-components to zero.

3, that it is imperative to get the correct positive and negative wave-number-components

3b, to get them right it most likely best not to shift the zero-frequency component to the centre - instead it is better to count up from 1 and down from end when identifying the interfering components.

That's it.

[AP]

Awesome, thank you so much for the suggestion. Will have a look slowly.

[Bjorn Gustavsson]

Before you do anything you first have to look at the image. In detail. That also means look at the intensity in a couple of selected rows and columns. You should first do something like this:

Im = imread('project1.png');
plot(Im(12,:),'.-')

and take a good look at what the signal looks like. Then have a think about what that tells you about the added obfuscating signal. To me it looks like someone addded a signal something like this:

[x,y] = meshgrid(1:size(Im,2),1:size(Im,1));
kx = 2*pi/20.5;
ky = 2*pi/20.5;
imNoisy = Im + A*max(0,sin(kx*x+ky*y));

That type of noise will be a pain to remove, but you'll figure it out.

[AP]

Thank you for that. But the noise added had been given, its a A= 20 * sin(0.1*pi*x+0.1*pi*y). As i didnt have the original image, so its a bit hard to see what the 'spike' different on the DFT of the original image and the corrupted image. Any oder suggestion that I can have a look?

Now on process to read the section that you suggest before & this problem of Weiner Filter.

[Bjorn Gustavsson]

Well that is claimed to be the noise that's added. You clearly did not do what I adviced you to and plot the 12th row and look at the intensity-variation there. If you had done so you would've seen what I saw, and I am perfectly comfortable writing that the claim of the noise being 20 * sin(0.1*pi*x+0.1*pi*y) is not true. This is very clear to see. The reason for this is that in the addition there was an impact due to type-casting, intended or not. Your teacher have given you a task that is unsuitable for the intended learning objective. But hey! you'll have an opportunity to learn this in another way and something else!

Do this:

x = 1:size(Im,2);
y = 1:size(Im,1);
[x,y] = meshgrid(x,y);
A = 20 * sin(0.1*pi*x+0.1*pi*y);
B = max(0,A);
subplot(2,2,1)
imagesc(A)
subplot(2,2,2)
imagesc(B)
subplot(2,2,3)
imagesc(log10(abs(fftshift(fft2(A)))))
subplot(2,2,4)
imagesc(log10(abs(fftshift(fft2(B)))))

This should give you some insight into what the difference in the two noise-images have on the fft of the noisy signal, and the complication of removing the second noise-type - it has to do with the decrease of the amplitudes of Fourier-coefficients of signals with continuous gradients and with discontinuous gradients.

For this type of noise (when you've been given the (incorrect-but-close) shape of the noise you can easily remove it. But that requires very little manual work for one image, but might be tricky to automate.

HTH

[AP]

Noted, thank you so much!! Manage to remove it by now.

[Bjorn Gustavsson]

My pleasure.

The lesson is (as it always seems to be for me): Really look at the data before going to work...

...I tried the automatic notch-filtering first too...

Answer 2

[Image Analyst]

See attached demo.

Because your bars are so sharp I'd probably rather fix the image in the spatial domain rather than the frequency domain.

[AP]

Thank you for the suggestion, try to work on it now. But based on the coding given wondering where the amplitudeThreshold = 10.9; get from and brightSpikes(115:143, :) = 0;? I mean the value 10.9 and 115:143 from? Thanks

[Image Analyst]

I'd have to look at it again, but that zeros out those rows in the brightSpikes matrix. I probably identified locations by looking at them in the display. If the image starts as 256x256, then those rows would be either the highest or lowest spatial frequencies, depending on if fftshift() was used.

[AP]

Good good, figure it out now. But this is what I get after doing filter. Still didn't enhance the image. Any others suggestion for that?

[Image Analyst]

I warned you that spatial filtering in the Fourier domain would not work very well for that type of noise. The reason is that the bars noise has lots of frequencies so using a notch/band pass filter, or even filtering out spikes by intensity would not get all the energy of the bars. So you'll filter out some of the bars, but not all, and you'll filter out some of the desired content, but not all. So it won't be great. But if you do it in the spatial domain you can do a better job. You can identify the angle of the lines with something like the radon transform, and then examine the profile perpendicular to the lines to determine if the intensity was a fraction of the unaltered intensity, or if it's an offset subtracted from the unaltered intensity (multiplicative vs. additive noise).