# Image Processing with Wavelets¶


This numerical tour overviews the use of wavelets for image approximation and denoising.

In [2]:
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/introduction_5_wavelets_2d')


## Wavelet Approximation¶

Image approximation is obtained by thresholding wavelets coefficients.

First we load an image $f \in \mathbb{R}^N$ of $N = n_0 \times n_0$ pixels.

In [3]:
name = 'cortex';
n0 = 512;
f = load_image(name,n0);
f = rescale( sum(f,3) );


Display it.

In [4]:
clf;
imageplot(f);


An orthogonal wavelet basis $\mathcal{B} = \{ \psi_{j,n}^k \}_{j,n}$ of $\mathbb{R}^N$ is composed of multiscale atoms parameterized by their scale $2^j$, position $2^j n \in [0,1]^2$ and orentation $k \in \{H,V,D\}$.

A forward wavelet transform computes the set of inner products $$\Psi f = \{ \langle f,\psi_{j,n}^k\rangle \} \in \mathbb{R}^N$$ where the wavelet atoms are defined as $$\psi_{j,n}^k(x) = \psi^k\left( \frac{x - 2^j n}{2^j} \right).$$

Set the minimum scale for the transform.

In [5]:
Jmin = 0;


A short-cut for the wavelet transform $\Psi$:

In [6]:
Psi = @(f)perform_wavelet_transf(f,Jmin,+1);


A short-cut for the inverse wavelet transform $\Psi^{-1} = \Psi^*$:

In [7]:
PsiS = @(fw)perform_wavelet_transf(fw,Jmin,-1);


Perform the wavelet transform to compute $\Psi f$.

In [8]:
fW = Psi(f);


Display the transformed coefficients.

In [9]:
clf;
plot_wavelet(fW);


To perform non-linear image approximation, one remove the small amplitude coefficients. This is performed using a hard thresholding $$H_T(f,\mathcal{B}) = \Psi^{-1} \circ H_T \circ \Psi (f) = \sum_{|\langle f,\psi_{j,n}^k\rangle| > T} \langle f,\psi_{j,n}^k\rangle \psi_{j,n}^k.$$

$T$ should be adapted to ensure a given number $M$ of non-zero coefficients, and then $f_M = H_T(f,\mathcal{B})$ is the best $M$ terms approximation of $f$ in $\mathcal{B}$.

Select a threshold.

In [10]:
T = .5;


Shortcut for the thresholding operator $H_T$.

In [11]:
Thresh = @(fW,T)fW .* (abs(fW)>T);


Perform hard thresholding of the coefficients.

In [12]:
fWT = Thresh(fW,T);


Exercise 1

Compute the ratio $M/N$ of non-zero coefficients.

In [13]:
exo1()

Ratio of non-zero coefficients : M/N = 0.00481.

In [14]:
%% Insert your code here.


Display the thresholded coefficients.

In [15]:
clf;
subplot(1,2,1);
plot_wavelet(fW);
title('Original coefficients');
subplot(1,2,2);
plot_wavelet(fWT);


Perform reconstruction using the inverse wavelet transform $\Psi^*$.

In [16]:
f1 = PsiS(fWT);


Display approximation.

In [17]:
clf;
imageplot(f, 'Image', 1,2,1);
imageplot(clamp(f1), strcat(['Approximation, SNR=' num2str(snr(f,f1),3) 'dB']), 1,2,2);


Number of coefficients for the approximation.

In [18]:
M = n0^2/16;


Exercise 2

Compute a threshold $T$ to keep only $M$ coefficients.

In [19]:
exo2()

In [20]:
%% Insert your code here.


Perform hard thresholding.

In [21]:
fWT = Thresh(fW,T);


Check the number of non-zero coefficients in |fWT|.

In [22]:
disp(strcat(['      M=' num2str(M)]));
disp(strcat(['|fWT|_0=' num2str(sum(fWT(:)~=0))]));

      M=16384
|fWT|_0=16384


Exercise 3

Compute an approximation with an decreasing number of coefficients.

In [23]:
exo3()

In [24]:
%% Insert your code here.


## Orthognal Wavelet Denoising¶

Image denoising is obtained by thresholding noisy wavelets coefficielts.

Here we consider a simple setting where we intentionnaly add some noise $w$ to a clean image $f$ to obtain $y = f + w$.

A denoiser computes an estimate $\tilde f$ of $f$ from the observations $y$ alone. In the mathematical model, since $y$ is a random variable depending on $w$, so is $\tilde f$. A mathematical evaluation of the efficiency of the denoiser is the average risk $E_w( \|f-\tilde f\|^2 )$.

Here we consider a single realization of the noise, so we replace the risk by the oracle error $\|f-\tilde f\|^2$. This allows us to bench the efficiency of the denoising methods by comparing the result to $f$. But you have to keep in mind that for real application, one does not have access to $f$.

We consider a Gaussian white noise $w$ of variance $\sigma^2$.

In [25]:
sigma = .1;


Generate a noisy image.

In [26]:
y = f + randn(n0,n0)*sigma;


Display.

In [27]:
clf;
imageplot(f, 'Clean image', 1,2,1);
imageplot(clamp(y), ['Noisy image, SNR=' num2str(snr(f,y),3) 'dB'], 1,2,2);


A denoising is obtained by thresholding the wavelet coefficients $$\tilde f = H_T(f,\mathcal{B}).$$

The asymptotically optimal threshold of Donoho and Johnstone is $T = \sqrt{2 \log(N)} \sigma$. In practice, one observes that much better result are obtained using $T \approx 3 \sigma$.

Compute the noisy wavelet coefficients.

In [28]:
fW = Psi(y);


Compute the threshold value using the $3\sigma$ heuristic.

In [29]:
T = 3*sigma;


Perform hard thresholding.

In [30]:
fWT = Thresh(fW,T);


Display the thresholded coefficients.

In [31]:
clf;
subplot(1,2,1);
plot_wavelet(fW);
title('Original coefficients');
subplot(1,2,2);
plot_wavelet(fWT);


Perform reconstruction.

In [32]:
f1 = PsiS(fWT);


Display denoising.

In [33]:
clf;
imageplot(clamp(y), 'Noisy image', 1,2,1);
imageplot(clamp(f1), strcat(['Denoising, SNR=' num2str(snr(f,f1),3) 'dB']), 1,2,2);