2D Haar Wavelets
This numerical tour explores 2D multiresolution analysis with the Haar transform.
Contents
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;
