Image Denoising with Linear Methods

This numerical tour introduces some basics about image denoising.

Installing toolboxes and setting up the path.
Noise in Signal and Image
Linear Signal Denoising
Linear Image Denoising
Wiener filtering

Installing toolboxes and setting up the path.

You need to download the following files: signal toolbox and general toolbox.

You need to unzip these toolboxes in your working directory, so that you have toolbox_signal and toolbox_general in your directory.

For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.

Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.

Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');

Noise in Signal and Image

In these numerical tour, we simulate noisy acquisition by adding some white noise (each pixel is corrupted by adding an independant Gaussian variable). This is helpful since we know the original, clean, image, and can test and compare several algorihtms by computing the recovery error.

We load a signal.

name = 'piece-regular';
n = 1024;
x0 = load_signal(name,n);
x0 = rescale(x0);

We add some noise to it.

sigma = .04; % noise level
x = x0 + sigma*randn(size(x0));
clf;
subplot(2,1,1);
plot(x0); axis([1 n -.05 1.05]);
subplot(2,1,2);
plot(x); axis([1 n -.05 1.05]);

We load an image.

n = 256;
name = 'hibiscus';
M0 = load_image(name,n);
M0 = rescale( sum(M0,3) );

Then we add some gaussian noise to it.

sigma = .08; % noise level
M = M0 + sigma*randn(size(M0));
clf;
imageplot(M0, 'Original', 1,2,1);
imageplot(clamp(M), 'Noisy', 1,2,2);

Linear Signal Denoising

A translation invariant linear denoising is necessarely a convolution with a kernel h. It correspond to a linear diagonal operation over the Fourier domain that multiply each noisy Fourier coefficient by the Fourier transform of h.

In practice, one uses a Gaussian fitler h, and the only parameter is the width (variance) of the filter.

% width of the filter
mu = 4;
% compute a Gaussian filter of width mu
t = (-length(x)/2:length(x)/2-1)';
h = exp( -(t.^2)/(2*mu^2) );
h = h/sum(h(:));

The Fourier transform of a Gaussian discrete filter is nearly a Gaussian whose width is proportional to 1/mu.

% Fourier transform of the (centered) filter
hf = real( fft(fftshift(h)) );
hf = fftshift(hf);
% display
clf;
subplot(2,1,1);
plot( t,h  ); axis('tight');
title('Filter h');
subplot(2,1,2);
plot( t,hf  ); axis('tight');
title('Fourier transform');

Since we use periodic boundary condition, the convolution of x with h can be computer over the Fourier domain.

% Fourier coefficients of the noisy signal
xf = fft(x);
% Fourier coefficients of the denoised signal
xhf = xf .* fft(fftshift(h));
% Denoised signal
xh = real( ifft(xhf) );

We display the denoised signal. Although most of the noise is removed, the singularity have been blurred.

clf;
subplot(2,1,1);
plot( t,x  ); axis('tight');
title('Noisy');
subplot(2,1,2);
plot( t,xh  ); axis('tight');
title('Denoised');

We display the noisy and denoised Fourier coefficients. One can see that the denoising remove the high frequency coefficients.

% log of Fourier transforms
epsilon = 1e-10;
L0 = log10(epsilon+abs(fftshift(fft(x0))));
L = log10(epsilon+abs(fftshift(xf)));
Lh = log10(epsilon+abs(fftshift(xhf)));
% display Fourier transforms
clf;
subplot(2,1,1);
plot( t, L, '-' );
axis([-length(x)/2 length(x)/2 -4 max(L)]);
title('log of noisy Fourier coefs.');
subplot(2,1,2);
plot( t, Lh, '-' );
axis([-length(x)/2 length(x)/2 -4 max(L)]);
title('log of denoised Fourier coefs.');

It is non-trivial to select the width parameter mu to minimize the denoising error. It should account for both the variance of the noise and the power spectrum of the image.

We display the blurring for an increasing value of mu.

mulist = linspace(.5,4,4);
clf;
for i=1:length(mulist)
    mu = mulist(i);
    % compute the filter
    h = exp( -(t.^2)/(2*mu^2) );
    h = h/sum(h(:));
    % perform the blurring
    xh = real( ifft( fft(x) .* fft(fftshift(h)) ));
    subplot( 4,1,i );
    plot(t, clamp(xh) ); axis('tight');
    title(strcat(['\mu=' num2str(mu)]));
end

Exercice 1: (the solution is exo1.m) Try for various Gaussian variance mu to compute the denoising xh. Compute, in an oracle manner, the best variance muopt by computing the residual error snr(x0,xh).

exo1;

The optimal smoothing width is 1.2333 pixels, SNR=24.8113dB.

Display the results.

% compute the optimal filter
h = exp( -(t.^2)/(2*muopt^2) );
h = h/sum(h(:));
% perform blurring
xh = real( ifft( fft(x) .* fft(fftshift(h)) ));
% display
clf;
subplot(2,1,1);
plot(t, clamp(x)); axis('tight');
title('Noisy');
subplot(2,1,2);
plot(t, clamp(xh)); axis('tight');
title('Denoised');

Linear Image Denoising

We denoise similarely a 2D image using a 2D Gaussian filter whose width mu is optimized to match the noise level and the regularity of the signal.

We use a simple gaussian blur to denoise an image.

% we use cyclic boundary condition since it is quite faster
options.bound = 'per';
% number of pixel of the filter
mu = 10;
Mh = perform_blurring(M,mu,options);
clf;
imageplot(clamp(M), 'Noisy', 1,2,1);
imageplot(clamp(Mh), 'Blurred', 1,2,2);

We display the blurring for an increasing value of mu.

mulist = linspace(3,15,6);
clf;
for i=1:length(mulist)
    mu = mulist(i);
    Mh = perform_blurring(M,mu,options);
    imageplot(clamp(Mh), strcat(['\mu=' num2str(mu)]), 2,3,i);
end

Exercice 2: (the solution is exo2.m) Try for various Gaussian variance to compute the denoising Mh. Compute, in an oracle manner, the best variance muopt by computing the residual error snr(M0,Mh).

exo2;

The optimal smoothing width is 3.53 pixels, SNR=21.4382dB.

Display the results

% optimal filter
Mgauss = perform_blurring(M,muopt,options);
% display
clf;
imageplot(M, strcat(['Noisy, SNR=' num2str(snr(M0,M)) 'dB']), 1,2,1);
imageplot(Mgauss, strcat(['Gaussian denoise, SNR=' num2str(snr(M0,Mgauss)) 'dB']), 1,2,2);

Wiener filtering

In a probabilistic setting, for translation invariant signal distributions, the Wiener filtering is the optimal filtering.

Perform the wiener filtering

[Mwien,Hwien] = peform_wiener_filtering(M0,M,sigma);

display the filter

k = 5;
clf;
imageplot(Hwien(n/2-k+2:n/2+k,n/2-k+2:n/2+k), 'Wiener filter (zoom)');

display the result

% display
clf;
imageplot( clamp(Mgauss), strcat(['Gaussian denoise, SNR=' num2str(snr(M0,Mgauss)) 'dB']), 1,2,1);
imageplot( clamp(Mwien), strcat(['Wiener denoise, SNR=' num2str(snr(M0,Mwien)) 'dB']), 1,2,2);