Periodic noise in images


Let us consider an image affected by a superposition of pure tones. As any periodic function can be decomposed as such superposition (the Fourier series), it is fair to call generically periodic noise to this situation although the randomness commonly associated to the idea of noise is not present here.

We are going to think that this superposition is of very few waves, modeling the case in which the Fourier coefficients have a quick decay. A natural method to de-noise the image is to examine its DFT and eliminate the isolated peaks, because they corresponds to pure tones.

In the octave code below, precisely in the function denoip, I have devised a method to detect automatically isolated peaks based on the simple idea that they have much larger value that the median of the surrounding values. There is a parameter alpha that quantifies the meaning of "much larger", namely it gives the threshold for the ratio between them. On the other hand, r embodies the scope of the "surrounding values", it indicates that we take the median on a square of size 2r+1. Once the spurious peaks in the DCT (I use it to avoid complex values) have been detected and substituted by the median, the inverse of the DCT gives the reconstruction.

The rest of the code is to apply this method to noisy versions of the following 512x512 photo. We take all the time r=2.

goose

We first consider the addition of cos(32*pi*X)/4 where X (resp. Y) is the horizontal (resp. vertical) coordinate taking the side of the photo as the unit length. The aspect of the resulting noisy image is:

pnois1

It is not needed any fine tuning of alpha, values between 4 and 5 and even in a wider range give similarly good results. For instance

pnois1_1
r=2,     alpha=4.5

Now we complicate the situation adding to the noisy image still more noise summing sin(10*pi*(X+Y))/4. The aspect of the result is:

pnois2

With alpha=5 and alpha=8 we get the following results. They are quite similar but note that in the second reconstruction the more noticeable periodic lights in the upper and lower parts are an artifact.

pnois2_1
r=2,     alpha=5

pnois2_2
r=2,     alpha=8

Finally, let us examine what happens when we put alpha to the limit. If it is very large then some peaks will survive and we will observe periodic noise in the reconstruction.

pnois2_r1
r=2,     alpha=80

On the other hand, if alpha is very small some values will be considered wrongly as peaks and replaced by other values. The outcome is a loss in the quality of the image.

pnois2_r2
r=2,     alpha=2

The code

This is the octave code that produces the images. It also works with matlab, if the signal processing toolbox is available, omitting the pkg lines.

clear all
close all


pkg load image
pkg load signal

ima = imread('../images/goose.jpg');
ima = im2double(ima);


x = linspace(0,1 ,512);
y = linspace(0,1 ,512);
[X, Y] = meshgrid(x, y);

iman1 = ima + 0.25*cos(32*pi*X);
iman2 = iman1 + 0.25*sin(10*pi*(X+Y));

% Noisy image 1 (scaled)
imwrite( im2double(mat2gray(iman1)) ,'./pnois1.jpg' );

% Noisy image 2 (scaled)
imwrite( im2double(mat2gray(iman2)) ,'./pnois2.jpg' );

% Noise 1, r=2, alpha = 4.5
res = denoip( iman1, 2, 4.5 );
imwrite( res ,'./pnois1_1.jpg' );

% Noise 2, r=2, alpha = 5
res = denoip( iman2, 2, 5 );
imwrite( res ,'./pnois2_1.jpg' );

% Noise 2, r=2, alpha = 8
res = denoip( iman2, 2, 8 );
imwrite( res ,'./pnois2_2.jpg' );

% Noise 2, r=2, alpha = 80
res = denoip( iman2, 2, 80 );
imwrite( res ,'./pnois2_r1.jpg' );

% Noise 2, r=2, alpha = 2
res = denoip( iman2, 2, 2 );
imwrite( res ,'./pnois2_r2.jpg' );

It calls to the function denoip:

% Clean periodic noise
% iman = noisy image
% r = media filter of size 2r+1 x 2r+1
% alpha = threshold to consider a peak

function res = denoip( iman, r, alpha )
  % DCT
  J = dct2(iman);

  % mask of peaks
  % It is a peak if it is not in the first corner
  % and takes a value very large in comparison
  % to the median
  aJ = abs(J);
  % mask = (aJ > alpha*medfilt2(aJ,true(2*r+1), 'symmetric' ));
  mask = (aJ > alpha*medfilt2(aJ,[2*r+1,2*r+1], 'symmetric' ));
  mask(1:r,1:r) = zeros(r);

  % Cleaned values of the DCT
  % nJ = (1-mask).*J + mask.*medfilt2(J, logical(ones(2*r+1)) );
  nJ = (1-mask).*J + mask.*medfilt2(J, [2*r+1,2*r+1] );
  
  % Invert the DCT an normalizes just in case
  res = idct2(nJ);
  res = im2double(mat2gray(res));
  
end