Download All

Code covered by the BSD License  

Highlights from
A Numerical Tour of Signal Processing

from A Numerical Tour of Signal Processing by Gabriel Peyre
A set of Matlab experiments that illustrates advanced computational signal and image processing.

Image Denoising with Wavelets

Image Denoising with Wavelets

This numerical tour uses wavelets to perform both linear and non-linear image denoising.

Contents

Installing toolboxes and setting up the path.

You need to download the general purpose toolbox and the signal toolbox.

You need to unzip these toolboxes in your working directory, so that you have toolbox_general/ and toolbox_signal/ 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 commands 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(path,p); % scilab users must *not* execute this

Then you can add these toolboxes to the path.

% Add some directories to the path
getd('toolbox_signal/');
getd('toolbox_general/');

Image loading and adding Gaussian noise.

A simple noise model is additive Gaussian noise.

First we load an image.

name = 'boat';
n = 256;
M0 = rescale( load_image(name,n) );

Then we add some gaussian noise to it.

sigma = .08; % noise level
M = M0 + sigma*randn(size(M0));

Display.

clf;
imageplot(M0, 'Original', 1,2,1);
imageplot(clamp(M), 'Noisy', 1,2,2);

Hard Thresholding vs. Soft Thresholding

A thresholding is a 1D non-linear function applied to each wavelet coefficients. The most important thresholding are the hard thresholding (related to L0 minimization) and the soft thresholding (related to L1 minimization).

We compute here the 1D thresholding curve.

% threshold value
T = 1;
v = -linspace(-3,3,2000);
% hard thresholding of the t values
v_hard = v.*(abs(v)>T);
% soft thresholding of the t values
v_soft = max(1-T./abs(v), 0).*v;

Display them.

clf;
hold('on');
h = plot(v, v_hard);
if using_matlab()
    set(h, 'LineWidth', 2);
end
h = plot(v, v_soft, 'r--');
if using_matlab()
    set(h, 'LineWidth', 2);
end
axis('equal'); axis('tight');
legend('Hard thresholding', 'Soft thresholding');
hold('off');

Wavelet Denoising with Hard Thesholding

It is possible to perform non linear denoising by thresholding the wavelet coefficients. This allows to better respect the sharp features of the image.

Some parameters for the orthogonal wavelet transform.

options.ti = 0;
Jmin = 4;

First we compute the wavelet coefficients of the noisy image.

MW = perform_wavelet_transf(M,Jmin,+1,options);

Select the threshold value. In practice a threshold of 3*sigma is close to optimal for natural images. 

T = 3*sigma;

Then we hard threshold the coefficients below the noise level.

MWT = MW .* (abs(MW)>T);

Display the thresholded coefficients.

clf;
subplot(1,2,1);
plot_wavelet(MW,Jmin);
title('Noisy coefficients');
set_axis(0);
subplot(1,2,2);
plot_wavelet(MWT,Jmin);
title('Thresholded coefficients');
set_axis(0);

One can then reconstruct from these noisy coefficients.

Mhard = perform_wavelet_transf(MWT,Jmin,-1,options);

Display the denoising result.

clf;
imageplot(clamp(M), strcat(['Noisy, SNR=' num2str(snr(M0,M),3)]), 1,2,1);
imageplot(clamp(Mhard), strcat(['Hard denoising, SNR=' num2str(snr(M0,Mhard),3)]), 1,2,2);

Wavelet Denoising with Soft Thesholding

The image suffers from many artifacts (wavelets poping arround). It is possible to improve the result by using soft thresholding. Two important remark should be made

  • First one must use a lower threshold because soft thresholding also lower the value of non-thresholded coeffcients. In practice, a threshold of 3/2*sigma works well.
  • The low frequency part of the coefficients must not be thresholded.

Select the threshold.

T = 3/2*sigma;

Perform the soft thresholding.

MWT = perform_thresholding(MW,T,'soft');

Re-inject the low frequencies.

MWT(1:2^Jmin,1:2^Jmin) = MW(1:2^Jmin,1:2^Jmin);

Re-construct.

Msoft = perform_wavelet_transf(MWT,Jmin,-1,options);

Display the denoising result.

clf;
imageplot(clamp(Mhard), strcat(['Hard denoising, SNR=' num2str(snr(M0,Mhard),3)]), 1,2,1);
imageplot(clamp(Msoft), strcat(['Soft denoising, SNR=' num2str(snr(M0,Msoft),3)]), 1,2,2);

Exercice 1: (the solution is exo1.m) Determine the best threshold T for both hard and soft thresholding. Test several T values in [.8*sigma, 4.5*sigma]. Do not forget to keep the low frequency wavelet coefficients. What can you conclude from these results ? Test with another image. 

exo1;

Translation Invariant Denoising with Cycle Spinning

Orthogonal wavelet transforms are not translation invariant. It means that the processing of an image and of a translated version of the image give different results.

Any denoiser can be turned into a translation invariant denoiser by performing a cycle spinning. The denoiser is applied to several shifted copies of the image, then the resulting denoised image are shifted back to the original position, and the results are averaged.

Number m of translations. A perfect invariance is obtained for m=2^Jmin, but here we enforce only approximate translation invariance. 

m = 4;

Initialize the denoised image.

Mti = zeros(n,n);

Generate a set of shift.

[dY,dX] = meshgrid(0:m-1,0:m-1);

Here is a typical application of one step of the scheme.

The index i should run in 1,...,m^2

i = 10;

Apply the shift, using circular boundary conditions.

Ms = circshift(M,[dX(i) dY(i)]);

Apply here the denoising to Ms.

After denoising, do the inverse shift.

Ms = circshift(Ms,-[dX(i) dY(i)]);

Accumulate the result to obtain at the end the denoised image that averahe the translated results.

Mti = (i-1)/i*Mti + 1/i*Ms;

Exercice 2: (the solution is exo2.m) Perform the cycle spinning denoising by iterating on i. 

exo2;

Exercice 3: (the solution is exo3.m) Study the influence of the number m of shift on the denoising quality 

exo3;

Translation Invariant Wavelet Transform

Another way to achieve translation invariance, which is equivalent to performing cycle spinning, is to replace the orthogonal wavelet transform by a redundant translation invariant wavelet transform.

The chief advantage of this approach is that it is faster to compute, since an invariance to m=2^J translations is obtained with an O((3*J+1)*n^2) algorithm.

A translation invariant wavelet transform is implemented by ommitting the sub-sampling at each stage of the transform. This correspond to the decomposition of the image in a redundant familly of n^2*(3*J+1) atoms where n^2 is the number of pixels and J is the number of scales of the transforms.

For Scilab, we need to extend a little the available memory.

extend_stack_size(4);

The invariant transform is obtained using the same function, by activating the switch options.ti=1. 

options.ti = 1;
MW = perform_wavelet_transf(M0,Jmin,+1,options);

MW(:,:,1) corresponds to the low scale residual. Each MW(:,:,3*j+k+1) for k=1:3 (orientation) corresponds to a scale of wavelet coefficient, and has the same size as the original image. 

clf;
i = 0;
for j=1:2
    for k=1:3
        i = i+1;
        imageplot(MW(:,:,i+1), strcat(['Scale=' num2str(j) ' Orientation=' num2str(k)]), 2,3,i );
    end
end

Translation Invariant Wavelet Denoising

Orthogonal wavelet denoising does not performs very well because of its lack of translation invariance. A much better result is obtained by not sub-sampling the wavelet transform, which leads to a redundant tight-frame.

First we compute the translation invariant wavelet transform

options.ti = 1;
MW = perform_wavelet_transf(M,Jmin,+1,options);

Then we threshold the set of coefficients.

T = 3.5*sigma;
MWT = perform_thresholding(MW,T,'hard');

We can display some wavelets coefficients

J = size(MW,3)-5;
clf;
imageplot(MW(:,:,J), 'Noisy coefficients', 1,2,1);
imageplot(MWT(:,:,J), 'Thresholded coefficients', 1,2,2);

We can now reconstruct

Mti = perform_wavelet_transf(MWT,Jmin,-1,options);

Display the denoising result.

clf;
imageplot(clamp(Msoft), strcat(['Soft orthogonal, SNR=' num2str(snr(M0,Msoft),3)]), 1,2,1);
imageplot(clamp(Mti), strcat(['Hard invariant, SNR=' num2str(snr(M0,Mti),3)]), 1,2,2);

Exercice 4: (the solution is exo4.m) Determine the best threshold T for both hard and soft thresholding, but now in the translation invariant case. What can you conclude ? 

exo4;


Copyright ® 2008 Gabriel Peyre

Contact us at files@mathworks.com