# 2D Haar Wavelets

This numerical tour explores 2D multiresolution analysis with the Haar transform.

## Contents

- Installing toolboxes and setting up the path.
- Forward 2D Haar transform
- Inverse 2D Haar transform.
- 2D Haar Image Approximation

## Installing toolboxes and setting up the path.

You need to download the general purpose toolbox and the signal toolbox.

You need to unzip these toolboxes in your working directory, so that you have `toolbox_general/` and `toolbox_signal/` 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 commands 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(path,p); % scilab users must *not* execute this

Then you can add these toolboxes to the path.

    % Add some directories to the path
    getd('toolbox_signal/');
    getd('toolbox_general/');

## Forward 2D 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 2D, 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).

First we load an image.

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

Initialize the transformed coefficients as the image itself and set the initial scale as the maximum one. `MW` will be iteratively transformated and will contains the coefficients.

    MW = M;
    j = log2(n)-1;

Select the sub-part of the image to transform.

    A = MW(1:2^(j+1),1:2^(j+1));

Compute average and differences along the vertical direction.

    Coarse = ( A(1:2:size(A,1),:) + A(2:2:size(A,1),:) )/sqrt(2);
    Detail = ( A(1:2:size(A,1),:) - A(2:2:size(A,1),:) )/sqrt(2);

Concatenate them in the vertical direction to get the result.

    A = cat3(1, Coarse, Detail );

Display the result of the vertical transform.

    clf;
    imageplot(M,'Original image',1,2,1);
    imageplot(A,'Vertical transform',1,2,2);

Compute average and differences along the horizontal direction.

    Coarse = ( A(:,1:2:size(A,1)) + A(:,2:2:size(A,1)) )/sqrt(2);
    Detail = ( A(:,1:2:size(A,1)) - A(:,2:2:size(A,1)) )/sqrt(2);

Concatenate them in the horizontal direction to get the result.

    A = cat3(2, Coarse, Detail );

Assign the transformed data.

    MW(1:2^(j+1),1:2^(j+1)) = A;

Display the result of the horizontal transform.

    clf;
    imageplot(M,'Original image',1,2,1);
    subplot(1,2,2);
    plot_wavelet(MW,log2(n)-1); title('Transformed')

Exercice 1: (the solution is exo1.m) Implement a full wavelet transform that extract iteratively wavelet coefficients, by repeating these steps. Take care of choosing the correct number of steps.

    exo1;

Check for orthogonality of the transform (conservation of energy).

    disp(strcat(['Energy of the signal       = ' num2str(norm(M(:)).^2)]));
    disp(strcat(['Energy of the coefficients = ' num2str(norm(MW(:)).^2)]));

    Energy of the signal       = 12708.0303
    Energy of the coefficients = 12708.0303

Display the wavelet coefficients.

    clf;
    subplot(1,2,1);
    imageplot(M); title('Original');
    subplot(1,2,2);
    plot_wavelet(MW, 1); title('Transformed');

## Inverse 2D Haar transform.

Inversing the Haar transform means retrieving a signal `M1` from the coefficients `MW`. If `MW` are exactely the coefficients of `M`, then `M=M1` up to machine precision.

Initialize the image to recover `M1` as the transformed coefficient, and select the smallest possible scale.

    M1 = MW;
    j = 0;

Select the sub-coefficient to transform.

    A = M1(1:2^(j+1),1:2^(j+1));

Retrieve coarse and detail coefficients in the vertical direction (you can begin by the other direction, this has no importance).

    Coarse = A(1:2^j,:);
    Detail = A(2^j+1:2^(j+1),:);

Undo the transform by sum and difference.

    A(1:2:size(A,1),:) = ( Coarse + Detail )/sqrt(2);
    A(2:2:size(A,2),:) = ( Coarse - Detail )/sqrt(2);

Retrieve coarse and detail coefficients in the horizontal direction.

    Coarse = A(:,1:2^j);
    Detail = A(:,2^j+1:2^(j+1));

Undo the transform by sum and difference.

    A(:,1:2:size(A,1)) = ( Coarse + Detail )/sqrt(2);
    A(:,2:2:size(A,2)) = ( Coarse - Detail )/sqrt(2);

Assign the result.

    M1(1:2^(j+1),1:2^(j+1)) = A;

Exercice 2: (the solution is exo2.m) Write the inverse wavelet transform that computes `M1` from the coefficients `MW`. Compare `M1` with `M`.

    exo2;

Check that we recover exactly the original image.

    disp(strcat((['Error |M-M1|/|M| = ' num2str(norm(M(:)-M1(:))/norm(M(:)))])));

    Error |M-M1|/|M| = 2.0389e-15

## 2D Haar Image Approximation

To perform image processing, one usually modify the coefficients `MW`. Denoising and approximation are obtained by removing low amplitude coefficients (setting them to 0).

For instance, it is possible to threshold the coefficients in order to set to 0 the smallest coefficients.

First select a threshold value (the largest the threshold, the more agressive the approximation).

    T = .2;

Then set to 0 coefficients with magnitude below the threshold.

    MWT = MW .* (abs(MW)>T);

Display thresholded coefficients.

    clf;
    subplot(1,2,1);
    plot_wavelet(MW); axis('tight'); title('Original coefficients');
    subplot(1,2,2);
    plot_wavelet(MWT); axis('tight'); title('Thresholded coefficients');

Exercice 3: (the solution is exo3.m) Find the thresholds `T` so that the number `m` of remaining coefficients in `MWT` are `.05*n^2` and `.2*n^2`. Use this threshold to compute `MWT` and then display the corresponding approximation `M1` of `M`.

    exo3;