2-D Haar 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 explores 2-D multiresolution analysis with the Haar transform. It was introduced in 1910 by Haar Haar1910 and is arguably the first example of wavelet basis.

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

Forward 2-D Haar transform

The Haar transform is the simplest orthogonal wavelet transform. It is computed by iterating difference and averaging between odd and even samples of the signal. Since we are in 2-D, we need to compute the average and difference in the horizontal and then in the vertical direction (or in the reverse order, it does not mind).

Size $N=n \times n$ of the image.

In [3]:
n = 256;
N = n*n;

First we load an image.

In [4]:
name = 'hibiscus';
f = load_image(name,n);
f = rescale( sum(f,3) );

The Haar transform operates over $J = \log_2(n)-1$ scales. It computes a series of coarse scale and fine scale coefficients $a_j, d_j^H, d_j^V, d_j^D \in \RR^{n_j \times n_j}$ where $N_j=2^j$.

In [5]:
J = log2(n)-1;

The forward Haar transform computes $$ \Hh(f) = (d_j^\om)_{j=0,\ldots,J}^{\om=H,V,D} \cup \{a_0\} . $$ Note that the set of coarse scale coefficients $(a_j)_{j>0}$ are not stored.

This transform is orthogonal, meaning $ \Hh \circ \Hh^* = \text{Id} $, and that there is the following conservation of energy $$ \sum_i \abs{f_i}^2 = \norm{f}^2 = \norm{\Hh f}^2 = \sum_{j,\om} \norm{d_j^\om}^2 + \abs{a_0}^2. $$

One initilizes the algorithm with $a_{J+1}=f$.

The first step apply a vertical transformtion, which corresponds to applying a 1-D Haar transform on each column, i.e. it computes $\tilde d_{j},\tilde a_j$ from $a_{j+1}$ as, for all $\ell=0,\ldots,2^{j+1}-1$ and $k=0,\ldots,2^j-1$, $$ \tilde a_j[k,\ell] = \frac{a_{j+1}[2k,\ell] + a_{j+1}[2k+1,\ell]}{\sqrt{2}}$ \qandq \tilde d_j[k,\ell] = \frac{a_{j+1}[2k,\ell] - a_{j+1}[2k+1,\ell]}{\sqrt{2}}. $$

Shortcut to compute $\tilde a_j, \tilde d_j$ from $a_{j+1}$.

In [6]:
haarV = @(a)[   a(1:2:length(a),:) + a(2:2:length(a),:);
                    a(1:2:length(a),:) - a(2:2:length(a),:) ]/sqrt(2);

Display the result of the vertical transform.

In [7]:
clf;
imageplot(f,'Original image',1,2,1);
imageplot(haarV(f),'Vertical transform',1,2,2);

One then apply the same 1-D horizontal transform to both $\tilde a_j, \tilde d_j$ to obtain the resulting matrices at scale $j$, i.e. to compute $a_j,d_j^H,d_j^V,d_j^D$ as, for all $\ell, k =0,\ldots,2^j-1$, $$ a_j[k,\ell] = \frac{\tilde a_{j}[k,2\ell] + \tilde a_{j}[k,2\ell+1]}{\sqrt{2}}$ \qandq d_j^H[k,\ell] = \frac{\tilde a_{j}[k,2\ell] - \tilde a_{j}[k,2\ell+1]}{\sqrt{2}}. $$ $$ d_j^V[k,\ell] = \frac{\tilde d_{j}[k,2\ell] + \tilde d_{j}[k,2\ell+1]}{\sqrt{2}}$ \qandq d_j^D[k,\ell] = \frac{\tilde d_{j}[k,2\ell] - \tilde d_{j}[k,2\ell+1]}{\sqrt{2}}. $$

Shortcut for the vertical transform.

In [8]:
haarH = @(a)haarV(a')';

Shortcut for one step of Haar transform.

In [9]:
haar = @(a)haarH(haarV(a));

Display the result of the first step of the algorithm.

In [10]:
clf;
imageplot(f,'Original image',1,2,1);
subplot(1,2,2);
plot_wavelet(haar(f),log2(n)-1); title('Transformed');

The output of the forward transform is stored in the variable |fw|. At a given iteration indexed by a scale |j|, it will store in |fw(1:2^(j+1),1:2^(j+1))| the variable $a_{j+1}$.

Initializes the algorithm at scale $j=J$.

In [11]:
j = J;

Initialize |fw| to the full signal.

In [12]:
fw = f;

At iteration indexed by $j$, select the sub-part of the signal containing $a_{j+1}$, and apply it the Haar operator.

In [13]:
fw(1:2^(j+1),1:2^(j+1)) = haar(fw(1:2^(j+1),1:2^(j+1)));

Exercise 1

Implement a full wavelet transform that extract iteratively wavelet coefficients, by repeating these steps. Take care of choosing the correct number of steps.

In [14]:
exo1()
In [15]:
%% Insert your code here.

Check that the transform is orthogonal, which means that the energy of the coefficient is the same as the energy of the signal.

In [16]:
fprintf('Should be 0: %.3f.\n', (norm(f(:))-norm(fw(:)))/norm(f(:)));
Should be 0: 0.000.

Display the wavelet coefficients.

In [17]:
clf;
subplot(1,2,1);
imageplot(f); title('Original');
subplot(1,2,2);
plot_wavelet(fw, 1); title('Transformed');

Inverse 2-D Haar transform.

The backward transform computes a signal $f_1 = \Hh^*(h)$ from a set of coefficients $ h = (d_j^\om)_{j,\om} \cup \{a_0\} $

If $h = \Hh(f)$ are the coefficients of an image $f$, one recovers $ f_1 = f $.

The inverse transform starts from $j=0$, and computes $a_{j+1}$ from the knowledge of $a_j,d_j$ by first inverting the horizontal transform, for all $k,\ell = 0,\ldots,2^j-1$, $$ \tilde a_{j}[k,2\ell] = \frac{ a_j[k,\ell] + a_j^H[k,\ell] }{\sqrt{2}}$ \qandq \tilde a_{j}[k,2\ell+1] = \frac{ a_j[k,\ell] - a_j^H[k,\ell] }{\sqrt{2}}$ $$ $$ \tilde d_{j}[k,2\ell] = \frac{ d_j^V[k,\ell] + d_j^D[k,\ell] }{\sqrt{2}}$ \qandq \tilde d_{j}[k,2\ell+1] = \frac{ d_j^V[k,\ell] - d_j^D[k,\ell] }{\sqrt{2}}. $$

In [18]:
ihaarV = @(a,d)assign( zeros(2*size(a,1),size(a,2)), ...
        [1:2:2*size(a,1), 2:2:2*size(a,1)], [a+d; a-d]/sqrt(2), 1 );

Then one invert the vertical transform, for all $k = 0,\ldots,2^j-1$ and for all $\ell=0,\ldots,2^{j+1}-1$ $$ a_{j}[2k,\ell] = \frac{ \tilde a_j[k,\ell] + \tilde d_j[k,\ell] }{\sqrt{2}}$ \qandq a_{j}[2k+1,\ell] = \frac{ \tilde a_j[k,\ell] - \tilde d_j[k,\ell] }{\sqrt{2}}$ $$

In [19]:
ihaarH = @(a,d)ihaarV(a',d')';

Shortcut to invert horizontal and then vertical transforms. This defines the invert of the Haar transform.

In [20]:
ihaar = @(a,dH,dV,dD)ihaarV( ihaarH(a,dH), ihaarH(dV,dD) );

Initialize the signal to recover $f_1$ as the transformed coefficient, and select the smallest possible scale $j=0$.

In [21]:
f1 = fw;
j = 0;

Apply one step of inverse transform.

In [22]:
s = 1:2^j; t = 2^j+1:2^(j+1); u = 1:2^(j+1);
f1(u,u) = ihaar(f1(s,s),f1(s,t),f1(t,s),f1(t,t));

Exercise 2

Write the inverse wavelet transform that computes $f_1$ from coefficients |fW|.

In [23]:
exo2()
In [24]:
%% Insert your code here.

Check that we recover exactly the original image.

In [25]:
fprintf('Should be 0: %.3f.\n', (norm(f-f1))/norm(f));
Should be 0: 0.000.

Haar Linear Approximation

Linear approximation is obtained by setting to zeros the Haar coefficient coefficients above a certain scale $j$.

Cut-off scale.

In [26]:
j = J-1;

Set to zero fine scale coefficients.

In [27]:
fw1 = zeros(n); fw1(1:2^j,1:2^j) = fw(1:2^j,1:2^j);

Computing the signal $f_1$ corresponding to these coefficients |fw1| computes the orthogonal projection on the space of images that are constant on disjoint squares of length $N/2^j \times N/2^j$.

Exercise 3

Display the reconstructed signal obtained from |fw1|, for a decreasing cut-off scale $j$.

In [28]:
exo3()
In [29]:
%% Insert your code here.

Haar Non-linear Approximation

Non-linear approximation is obtained by thresholding low amplitude wavelet coefficients. It is defined as $$ f_T = \Hh^* \circ S_T \circ \Hh(f) $$ where $S_T$ is the hard tresholding operator $$ S_T(x)_i = \choice{ x_i \qifq \abs{x_i}>T, \\ 0 \quad \text{otherwise}. }. $$

In [30]:
S = @(x,T) x .* (abs(x)>T);

Set the threshold value.

In [31]:
T = .1;

Threshold the coefficients.

In [32]:
fwT = S(fw,T);

Display the coefficients before and after thresholding.

In [33]:
clf;
subplot(1,2,1);
plot_wavelet(fw); axis('tight'); title('Original coefficients');
subplot(1,2,2);
plot_wavelet(fwT); axis('tight'); title('Thresholded coefficients');

Exercise 4

Find the threshold $T$ so that the number of remaining coefficients in |fwT| is a fixed number $m$. Use this threshold to compute |fwT| and then display the corresponding approximation $f_1$ of $f$. Try for an increasing number $m$ of coeffiients.

In [34]:
exo4()
In [35]:
%% Insert your code here.

Bibliography