Image Processing with Wavelets¶

Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\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}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

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