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.
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.
n = 256;
N = n*n;
First we load an image.
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$.
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}$.
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.
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.
haarH = @(a)haarV(a')';
Shortcut for one step of Haar transform.
haar = @(a)haarH(haarV(a));
Display the result of the first step of the algorithm.
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$.
j = J;
Initialize |fw| to the full signal.
fw = f;
At iteration indexed by $j$, select the sub-part of the signal containing $a_{j+1}$, and apply it the Haar operator.
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.
exo1()
%% 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.
fprintf('Should be 0: %.3f.\n', (norm(f(:))-norm(fw(:)))/norm(f(:)));
Should be 0: 0.000.
Display the wavelet coefficients.
clf;
subplot(1,2,1);
imageplot(f); title('Original');
subplot(1,2,2);
plot_wavelet(fw, 1); title('Transformed');