\[ \newcommand{\NN}{\mathbb{N}} \newcommand{\CC}{\mathbb{C}} \newcommand{\GG}{\mathbb{G}} \newcommand{\LL}{\mathbb{L}} \newcommand{\PP}{\mathbb{P}} \newcommand{\QQ}{\mathbb{Q}} \newcommand{\RR}{\mathbb{R}} \newcommand{\VV}{\mathbb{V}} \newcommand{\ZZ}{\mathbb{Z}} \newcommand{\FF}{\mathbb{F}} \newcommand{\KK}{\mathbb{K}} \newcommand{\UU}{\mathbb{U}} \newcommand{\EE}{\mathbb{E}} \newcommand{\Aa}{\mathcal{A}} \newcommand{\Bb}{\mathcal{B}} \newcommand{\Cc}{\mathcal{C}} \newcommand{\Dd}{\mathcal{D}} \newcommand{\Ee}{\mathcal{E}} \newcommand{\Ff}{\mathcal{F}} \newcommand{\Gg}{\mathcal{G}} \newcommand{\Hh}{\mathcal{H}} \newcommand{\Ii}{\mathcal{I}} \newcommand{\Jj}{\mathcal{J}} \newcommand{\Kk}{\mathcal{K}} \newcommand{\Ll}{\mathcal{L}} \newcommand{\Mm}{\mathcal{M}} \newcommand{\Nn}{\mathcal{N}} \newcommand{\Oo}{\mathcal{O}} \newcommand{\Pp}{\mathcal{P}} \newcommand{\Qq}{\mathcal{Q}} \newcommand{\Rr}{\mathcal{R}} \newcommand{\Ss}{\mathcal{S}} \newcommand{\Tt}{\mathcal{T}} \newcommand{\Uu}{\mathcal{U}} \newcommand{\Vv}{\mathcal{V}} \newcommand{\Ww}{\mathcal{W}} \newcommand{\Xx}{\mathcal{X}} \newcommand{\Yy}{\mathcal{Y}} \newcommand{\Zz}{\mathcal{Z}} \newcommand{\al}{\alpha} \newcommand{\la}{\lambda} \newcommand{\ga}{\gamma} \newcommand{\Ga}{\Gamma} \newcommand{\La}{\Lambda} \newcommand{\Si}{\Sigma} \newcommand{\si}{\sigma} \newcommand{\be}{\beta} \newcommand{\de}{\delta} \newcommand{\De}{\Delta} \renewcommand{\phi}{\varphi} \renewcommand{\th}{\theta} \newcommand{\om}{\omega} \newcommand{\Om}{\Omega} \renewcommand{\epsilon}{\varepsilon} \newcommand{\Calpha}{\mathrm{C}^\al} \newcommand{\Cbeta}{\mathrm{C}^\be} \newcommand{\Cal}{\text{C}^\al} \newcommand{\Cdeux}{\text{C}^{2}} \newcommand{\Cun}{\text{C}^{1}} \newcommand{\Calt}[1]{\text{C}^{#1}} \newcommand{\lun}{\ell^1} \newcommand{\ldeux}{\ell^2} \newcommand{\linf}{\ell^\infty} \newcommand{\ldeuxj}{{\ldeux_j}} \newcommand{\Lun}{\text{\upshape L}^1} \newcommand{\Ldeux}{\text{\upshape L}^2} \newcommand{\Lp}{\text{\upshape L}^p} \newcommand{\Lq}{\text{\upshape L}^q} \newcommand{\Linf}{\text{\upshape L}^\infty} \newcommand{\lzero}{\ell^0} \newcommand{\lp}{\ell^p} \renewcommand{\d}{\ins{d}} \newcommand{\Grad}{\text{Grad}} \newcommand{\grad}{\text{grad}} \renewcommand{\div}{\text{div}} \newcommand{\diag}{\text{diag}} \newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} } \newcommand{\pdd}[2]{ \frac{ \partial^2 #1}{\partial #2^2} } \newcommand{\dotp}[2]{\langle #1,\,#2\rangle} \newcommand{\norm}[1]{|\!| #1 |\!|} \newcommand{\normi}[1]{\norm{#1}_{\infty}} \newcommand{\normu}[1]{\norm{#1}_{1}} \newcommand{\normz}[1]{\norm{#1}_{0}} \newcommand{\abs}[1]{\vert #1 \vert} \newcommand{\argmin}{\text{argmin}} \newcommand{\argmax}{\text{argmax}} \newcommand{\uargmin}[1]{\underset{#1}{\argmin}\;} \newcommand{\uargmax}[1]{\underset{#1}{\argmax}\;} \newcommand{\umin}[1]{\underset{#1}{\min}\;} \newcommand{\umax}[1]{\underset{#1}{\max}\;} \newcommand{\pa}[1]{\left( #1 \right)} \newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. } \newcommand{\enscond}[2]{ \left\{ #1 \;:\; #2 \right\} } \newcommand{\qandq}{ \quad \text{and} \quad } \newcommand{\qqandqq}{ \qquad \text{and} \qquad } \newcommand{\qifq}{ \quad \text{if} \quad } \newcommand{\qqifqq}{ \qquad \text{if} \qquad } \newcommand{\qwhereq}{ \quad \text{where} \quad } \newcommand{\qqwhereqq}{ \qquad \text{where} \qquad } \newcommand{\qwithq}{ \quad \text{with} \quad } \newcommand{\qqwithqq}{ \qquad \text{with} \qquad } \newcommand{\qforq}{ \quad \text{for} \quad } \newcommand{\qqforqq}{ \qquad \text{for} \qquad } \newcommand{\qqsinceqq}{ \qquad \text{since} \qquad } \newcommand{\qsinceq}{ \quad \text{since} \quad } \newcommand{\qarrq}{\quad\Longrightarrow\quad} \newcommand{\qqarrqq}{\quad\Longrightarrow\quad} \newcommand{\qiffq}{\quad\Longleftrightarrow\quad} \newcommand{\qqiffqq}{\qquad\Longleftrightarrow\qquad} \newcommand{\qsubjq}{ \quad \text{subject to} \quad } \newcommand{\qqsubjqq}{ \qquad \text{subject to} \qquad } \]

Linear Image Denoising

Linear Image Denoising

This numerical tour introduces basic image denoising methods.

Contents

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

Noisy Image Formation

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 useful to test in an oracle maner the performance of our methods.

Size \(N = n \times n\) of the image.

n = 256;
N = n^2;

We load a clean image \(x_0 \in \RR^N\).

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

Display the clean image.

clf;
imageplot(x0);

Variance of the noise.

sigma = .08;

We add some noise to it to obtain the noisy signal \(y = x_0 + w\). Here \(w\) is a realization of a Gaussian white noise of variance \(\si^2\). 

y = x0 + sigma*randn(size(x0));

Display the clean and the noisy signals.

clf;
imageplot(clamp(y));

Linear Image Denoising

We consider a noising estimator \(x \in \RR^N\) of \(x_0\) that only depends on the observation \(y\). Mathematically speaking, it is thus a random vector that depends on the noise \(w\).

A translation invariant linear denoising is necessarely a convolution with a kernel \(h\) \[ x = x_0 \star h \] where the periodic convolution between two 2-D arrays is defined as \[ (a \star b)_i = \sum_j a(j) b(i-j). \]

It can be computed over the Fourier domain as \[ \forall \om, \quad \hat x(\om) = \hat x_0(\om) \hat h(\om). \]

cconv = @(a,b)real(ifft2( fft2(a).*fft2(b) ));

We use here a Gaussian fitler \(h\) parameterized by the bandwith \(\mu\).

normalize = @(h)h/sum(h(:));
t = [0:n/2-1, -n/2:-1]';
[Y,X] = meshgrid(t,t);
h = @(mu)normalize( exp( -(X.^2+Y.^2)/(2*mu^2) ) );

Display the filter \(h\) and its Fourier transform.

mu = 10;
clf;
subplot(2,1,1);
imageplot( fftshift( h(mu) )  ); axis('tight');
title('h');
subplot(2,1,2);
imageplot( fftshift( real(fft2(h(mu))) )  ); axis('tight');
title('fft2(h)');

Shortcut for the convolution with \(h\).

denoise = @(x,mu)cconv(h(mu), x);

Display a denoised signal.

clf;
imageplot( denoise(y,mu) );

Exercice 1: (check the solution) Display a denoised signal for several values of \(\mu\). 

exo1;

Exercice 2: (check the solution) Display the evolution of the oracle denoising error \( \norm{y-x_0} \) as a function of \(\mu\). Set \(\mu\) to the value of the optimal parameter. 

exo2;

Display the results.

clf;
imageplot( denoise(y,mu) );

Wiener Filtering

We suppose here that \(x_0\) is a realization of a random vector \(x_0\), whose distribution is Gaussian with a stationary covariance \(c\), and we denote \(P_{X_0}(\om) = \hat c(\om)\) the power-spectrum of \(x_0\).

Recall that \(w\) is a realization of a random vector \(W\) distributed according to \(\Nn(0,\si^2 \text{Id})\).

The (oracle) optimal filter minimizes the risk \[ R(h) = \EE_{W,X_0}( \norm{ X_0 - h \star (X_0 + W) }^2 ). \]

One can show that the solution of this problem, the so-called Wiener filter, is defined as \[ \forall \om, \quad \hat h(\om) = \frac{ P_{X_0}(\om) }{ P_{X_0}(\om) + \si^2 }. \]

We estimate \( P_{X_0} \) using the periodogram associated to the realization \(x_0\), i.e. \[ P_{X_0} \approx \frac{1}{N} \abs{\hat x_0}^2. \] 

P = 1/N * abs(fft2(x0)).^2;

Compute the approximate Wiener filter.

h_w = real( ifft2( P ./ ( P + sigma^2 ) ) );

Note that this is a theoretical filter, because in practice one does not have access to \(x_0\).

Display it.

clf;
imageplot(crop(fftshift(h_w),n/8)); axis tight;

Display the denoising result.

clf;
imageplot(cconv(y,h_w));

Note that this denoising is not very efficient, because the hypothesis of stationarity of \(X_0\) is not realistic for such piecewise-regular signal.


Copyright (c) 2010 Gabriel Peyre