# Image Denoising with Linear Methods

This numerical tour introduces some basics about image denoising.

## Contents

- 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);