2D Daubechies Wavelets
This numerical tour explores 2D multiresolution analysis with Daubchies wavelet 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/');
Wavelets Filters
Daubechies wavelets extends the haar wavelets by using more complicated average and difference operators. These operator are computed by filtering the signal and then sub-sampling by the resulting low/high pass filtered signal.
Each wavelet is uniquely determined by its low pass filter. The larger the size p=2*k of the filter, the higher is the number k of vanishing moment.
A high number of vanishing moments allows to better compress regular parts of the signal. However, increasing the number of vanishing moments also inceases the size of the support of the wavelets, wich can be problematic in part where the signal is singular (for instance discontinuous).
Choosing the best wavelet, and thus choosing k, that is adapted to a given class of signals, thus corresponds to a tradeoff between efficiency in regular and singular parts.
- The filter with k=1 vanishing moments corresponds to the Haar filter.
- The filter with k=2 vanishing moments corresponds to the famous D4 wavelet, which compresses perfectly linear signals.
- The filter with k=3 vanishing moments compresses perfectly quadratic signals.
Set the support size. To begin, we select the D4 filter.
p = 4;
Create the low pass filter h and the high pass g. We add a zero to ensure that it has a odd length. Note that the central value of h corresponds to the 0 position.
[h,g] = compute_wavelet_filter('Daubechies',p);
Note that the high pass filter g is computed directly from the low pass filter as.
g = [0 h(length(h):-1:2)] .* (-1).^(1:length(h))
Display.
disp(['h filter = [' num2str(h) ']']); disp(['g filter = [' num2str(g) ']']);
h filter = [0 0.48296 0.83652 0.22414 -0.12941] g filter = [-0 -0.12941 -0.22414 0.83652 -0.48296]
Up and Down Filtering
The basic wavelet operation is low/high filtering, followed by down sampling. For orthogonal filter, the reverse of this process is its dual (aka its transpose), which is upsampling followed by low/high pass filtering with the reversed filters and summing.
Create a dummy signal.
n = 256; f = rand(n,1);
Low/High pass filtering followed by sub-sampling.
a = subsampling( cconv(f,h) ); d = subsampling( cconv(f,g) );
Up-sampling followed by filtering.
f1 = cconv(upsampling(a),reverse(h)) + cconv(upsampling(d),reverse(g));
Check that we really recover the same signal.
disp(strcat((['Error |f-f1|/|f| = ' num2str(norm(f-f1)/norm(f))])));
Error |f-f1|/|f| = 5.4237e-13
Forward 2D Wavelet transform
A wavelet transform is computed by iterating high pass and loss pass filterings with h and g, followed by sub-samplings. Since we are in 2D, we need to compute these filterings+subsamplings 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;
M = rescale( load_image('lena',n) );
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));
Apply high and low filtering+subsampling in the vertical direction (1st ooordinate), to get coarse and details.
Coarse = subsampling(cconv(A,h,1),1); Detail = subsampling(cconv(A,g,1),1);
Note: subsamplling(A,1) is equivalent to A(1:2:end,:) and subsamplling(A,2) is equivalent to A(:,1:2:end).
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 imge',1,2,1); imageplot(A,'Vertical transform',1,2,2);
Apply high and low filtering+subsampling in the horizontal direction (2nd ooordinate), to get coarse and details.
Coarse = subsampling(cconv(A,h,2),2); Detail = subsampling(cconv(A,g,2),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 = 17020.7492 Energy of the coefficients = 17020.7492
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 Wavelet transform.
Inversing the wavelet 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 up-sampling and then dual filtering.
Coarse = cconv(upsampling(Coarse,1),reverse(h),1); Detail = cconv(upsampling(Detail,1),reverse(g),1);
Recover the coefficient by summing.
A = Coarse + Detail;
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 up-sampling and then dual filtering.
Coarse = cconv(upsampling(Coarse,2),reverse(h),2); Detail = cconv(upsampling(Detail,2),reverse(g),2);
Recover the coefficient by summing.
A = Coarse + Detail;
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| = 8.5496e-12
Linear 2D Wavelet Approximation
Linear approximation is performed by setting to zero the fine scale wawelets coefficients and then performing the inverse wavelet transform.
Here we keep only 1/16 of the wavelet coefficient, thus calculating an m term approximation with m=n^2/16.
eta = 4; MWlin = zeros(n,n); MWlin(1:n/eta,1:n/eta) = MW(1:n/eta,1:n/eta);
Exercice 3: (the solution is exo3.m) Compute and display the linear approximation Mlin obtained from the coefficients MWlin by inverse wavelet transform.
exo3;
Non-Linear 2D Wavelet Approximation
A non-linear m-term approximation is obtained by keeping only the m largest coefficients, which creates the smallest possible error.
Removing the smallest coefficient, to keep the m-largest, is equivalently obtainedby thresholding the coefficients 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 4: (the solution is exo4.m) Find the thresholds T so that the number m of remaining coefficients in MWT are m=n^2/16. Use this threshold to compute MWT and then display the corresponding approximation Mnlin of M. Compare this result with the linear approximation.
exo4;
Exercice 5: (the solution is exo5.m) Try with Different kind of wavelets, with an increasing number of vanishing moments.
exo5;
Separable 2D Wavelet Transform
A forward wavelet transform is obtained by applying the 1D wavelet transform to
Exercice 6: (the solution is exo6.m) Implement the foward separable transform. Wavelet transformm in 1D each column M(:,i) to obtain coefficients MWsep(:,i). Then re-transform each row MWsep(i,:)', and store the result in MW(i,:)'.
exo6;
Display the result.
clf; subplot(1,2,1); opt.separable = 0; plot_wavelet(MW,1,opt); title('Isotropic wavelets'); subplot(1,2,2); opt.separable = 1; plot_wavelet(MWsep,1,opt); title('Separable wavelets');
Exercice 7: (the solution is exo7.m) Implement the backward separable transform to recover an image M1 from the coefficients MWsep, which backward transform each row and then each columns.
exo7;
Check that we recover exactly the original image.
disp(strcat((['Error |M-M1|/|M| = ' num2str(norm(M(:)-M1(:))/norm(M(:)))])));
Error |M-M1|/|M| = 7.5175e-12
Exercice 8: (the solution is exo8.m) Perform m=n^2/16-terms best approximation with separable wavelet transform. Compare the result with isotropic wavelet approximation.
exo8;
The Shape of a Wavelet
An isotropic wavelet coefficient MW[k] corresponds to an inner product M,psi_{j,p}^s where k depends on the scale j and the position p and the orientation s (horizontal, vertical or diagonal).
The wavelet image M1 = psi_{j,p} is computed by applying the inverse wavelet transform to MW where MW[k]=1 and MW[l]=0 for l \neq k.
Exercice 9: (the solution is exo9.m) Compute wavelets at several scales and orientation. Here we show only horizontal wavelets, in 2D.
exo9;
Exercice 10: (the solution is exo10.m) Display Daubechies wavelets with different orientation, for different number of VM.
exo10;
