Wavelet Block Thresholding
This numerical tour presents block thresholding methods, that makes use of the structure of wavelet coefficients of natural images to perform denoising.
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/');
Generating a Noisy Image
Here we use an additive Gaussian noise.
First we load an image.
name = 'boat';
n = 256;
M0 = rescale( load_image(name,n) );
Noise level.
sigma = .08;
Then we add some Gaussian noise to it.
M = M0 + sigma*randn(size(M0));
Display.
clf; imageplot(M0, 'Original', 1,2,1); imageplot(clamp(M), 'Noisy', 1,2,2);
Orthogonal Wavelet Block Thresholding
Wavelet coefficients of natural images are not independant one from each other. One can thus improve the denoising results by thresholding block of coefficients togethers. Block thresholding is only efficient when used as a soft thresholder. Here we use a Stein soft thresholder.
Parameters for the orthogonal wavelet transform.
Jmin = 4; options.ti = 0;
First compute a wavelet transform.
MW = perform_wavelet_transf(M,Jmin,+1,options);
The block size.
w = 4;
Compute indexing of the blocks.
[dX,dY,X,Y] = ndgrid(0:w-1,0:w-1,1:w:n-w+1,1:w:n-w+1); I = X+dX + (Y+dY-1)*n;
Extract the blocks.
H = reshape(MW(I(:)),size(I));
Set the threshold value.
T = 1.25*sigma;
Compute the average energy of each block, and duplicate.
v = mean(mean(abs(H).^2,1),2); v = repmat( max3(v,1e-15), [w w]);
Threshold the blocks. We use here a Stein block thresholding. All values within a block are atenuated by the same factor.
HT = max3(1-T^2*v.^(-1),0) .* H;
Reconstruct the thresholded coefficients.
MWT = MW; MWT(I(:)) = HT(:);
For soft thresholding, we restore back the low frequencies.
MWT(1:2^Jmin,1:2^Jmin) = MW(1:2^Jmin,1:2^Jmin);
Display the thresholded coefficients.
clf; plot_wavelet(MWT,Jmin);
Exercice 1: (the solution is exo1.m) Display the evolution of the denoising SNR when T varies for both Stein soft thresholding, soft thresholding, and hard thresholding. Store in Mblock the best block thresholding denoising result.
exo1;
We compute the classical soft thresholding for comparison/
T = 1.3*sigma;
MWT = perform_thresholding(MW,T,'soft',options);
Msoft = perform_wavelet_transf(MWT,Jmin,-1,options);
Display a comparison of the two.
clf; imageplot(clamp(Msoft), strcat(['Soft orthogonal, SNR=' num2str(snr(M0,Msoft))]), 1,2,1); imageplot(clamp(Mblock), strcat(['Block thresholding, SNR=' num2str(snr(M0,Mblock))]), 1,2,2);
Exercice 2: (the solution is exo2.m) For the Stein block thresholding, study the influence of the block size.
exo2;
Translation invariant Block Thresholding
Block thresholding can also be applied to a translation invariant wavelet transform. It gives state of the art denoising results.
For Scilab user, you must increase the memory.
extend_stack_size(4);
Translation invariant transform.
options.ti = 1; MW = perform_wavelet_transf(M,Jmin,+1,options);
The block size.
w = 4;
Compute indexing of the blocks.
[dX,dY,X,Y,J] = ndgrid(0:w-1,0:w-1,1:w:n-w+1,1:w:n-w+1, 1:size(MW,3)); I = X+dX + (Y+dY-1)*n + (J-1)*n^2;
Extract the blocks.
H = reshape(MW(I(:)),size(I));
Then the thresholding proceeds as for the orthogonal denoising.
Exercice 3: (the solution is exo3.m) Compute the optimal threshold for TI-block thresholding. Do not forget to restore the low frequency coefficients after threshold. Store the best denoising result in Mti.
exo3;
Display.
clf; imageplot(clamp(Mblock), strcat(['Orthogonal, SNR=' num2str(snr(M0,Mblock))]), 1,2,1); imageplot(clamp(Mti), strcat(['Invariant, SNR=' num2str(snr(M0,Mti))]), 1,2,2);
