clf;subplot(2,1,1);plot(x0); axis([1 n -1 1]);set_graphic_sizes([], 20);title('Original signal');subplot(2,1,2);plot(x); axis([1 n -1 1]);set_graphic_sizes([], 20);title('Noisy signal');%%% _Exercice 1:_ (the solution is <../private/denoising_wavelet/exo1.m exo1.m>)% What is the optimal threshold |T| to remove as much as possible of noise% ? Try several values of |T|.exo1;%%% In order to be optimal without knowing in advance the amplitude of the% coefficients of |x0|, one needs to set |T| just above the noise level.% This means that |T| should be roughly equal to the maximum value of a% Gaussian white noise of size |n|.%%% _Exercice 2:_ (the solution is <../private/denoising_wavelet/exo2.m exo2.m>)% The theory predicts that the maximum of |n| Gaussian variable of variance |sigma^2|% is smaller than |sqrt(2*log(n))| with large probability (that tends to 1% when |n| increases). This is also a sharp result. Check this numerically% by computing with Monte Carlo sampling the maximum with |n| increasing% (in power of 2). Check also the deviation of the maximum when you% perform several trial with |n| fixed.exo2;%% Image loading and adding Gaussian noise.% A simple noise model is additive Gaussian noise.%%% First we load an image.name = 'boat';n = 256;M0 = load_image(name,n);M0 = rescale( M0, .05, .95 );%%% Then we add some gaussian noise to it.sigma = .08; % noise levelM = M0 + sigma*randn(size(M0));clf;imageplot(M0, 'Original', 1,2,1);imageplot(clamp(M), 'Noisy', 1,2,2);%% Hard Thresholding vs. Soft Thresholding%%% A thresholding is a 1D non-linear function applied to each wavelet% coefficients. The most important thresholding are the hard thresholding% (related to L0 minimization) and the soft thresholding (related to L1% minimization).% threshold valueT = 1;v = -linspace(-3,3,2000);% hard thresholding of the t valuesv_hard = v.*(abs(v)>T);% soft thresholding of the t valuesv_soft = max(1-T./abs(v), 0).*v;% displayclf;hold('on');plot(v, v_hard);plot(v, v_soft, 'rREPLACE_WITH_DASH_DASH');axis('equal'); axis('tight');legend('Hard thresholding', 'Soft thresholding');hold('off');%% Orthogonal Wavelet Denoising% It is possible to perform non linear denoising by thresholding the% wavelet coefficients. This allows to better respect the sharp features of% the image.%%% First we compute the wavelet coefficients% of the noisy image.options.ti = 0;Jmin = 4;MW = perform_wavelet_transf(M,Jmin,+1,options);%%% Then we hard threshold the coefficients below the noise level.% In practice a threshold of |3*sigma| is close to optimal for natural% images.T = 3*sigma;MWT = perform_thresholding(MW,T,'hard');clf;subplot(1,2,1);plot_wavelet(MW,Jmin);title('Noisy coefficients');set_axis(0);subplot(1,2,2);plot_wavelet(MWT,Jmin);title('Thresholded coefficients');set_axis(0);%%% One can then reconstruct from these noisy coefficients.Mhard = perform_wavelet_transf(MWT,Jmin,-1,options);% displayclf;imageplot(clamp(M), 'Noisy', 1,2,1);imageplot(clamp(Mhard), strcat(['Hard denoising, SNR=' num2str(snr(M0,Mhard))]), 1,2,2);%%% The image suffers from many artifacts (wavelets poping arround). It is% possible to improve the result by using soft thresholding. Two important% remark should be made%%% * First one must use a lower threshold because soft thresholding also lower the value of non-thresholded coeffcients. In practice, a threshold of |3/2*sigma| works well.% * The low frequency part of the coefficients must not be thresholded.T = 3/2*sigma;MWT = perform_thresholding(MW,T,'soft');% re-inject the low frequenciesMWT(1:2^Jmin,1:2^Jmin) = MW(1:2^Jmin,1:2^Jmin);% re-constructMsoft = perform_wavelet_transf(MWT,Jmin,-1,options);% displayclf;imageplot(clamp(Mhard), strcat(['Hard denoising, SNR=' num2str(snr(M0,Mhard))]), 1,2,1);imageplot(clamp(Msoft), strcat(['Soft denoising, SNR=' num2str(snr(M0,Msoft))]), 1,2,2);%%% _Exercice 3:_ (the solution is <../private/denoising_wavelet/exo3.m exo3.m>)% Determine the best threshold |T| for both hard and soft thresholding.% To that end, check for |T=alpha*sigma| (for hard) and |T=alpha*sigma/2|% (for hard) and compute the denoising error.% What can you conclude from these results ?% Test with another image.exo3;%% Estimating the noise level% In practice, the noise level |sigma| is unknown. % A good estimator is given by the median of the wavelet coefficients at% the finer scale. An even simple estimator is given by the normalized% derivate along X or Y direction%% % First we extract the high frequency residual.H = (M(1:n-1,:) - M(2:n,:))/sqrt(2);% histograms[h,t] = hist(H(:), 100);h = h/sum(h);% displayclf;imageplot(H, 'High freq. coefficients', 2,1,1);subplot(2,1,2);bar(t, h);axis([-.5 .5 0 max(h)]);%%% The mad estimator (median of median) must be rescaled so that % it gives the correct variance for gaussian noise.sigma_est = mad(H(:),1)/0.6745;disp( strcat(['Estimated noise level=' num2str(sigma_est), ', true=' num2str(sigma)]) );%% Translation Invariant Wavelet Transform% Orthogonal wavelet transforms are not translation invariant.% It means that the processing of an image and of a translated version of% the image give different results. A translation invariant wavelet% transform is implemented by ommitting the sub-sampling at each stage of% the transform. This correspond to the decomposition of the image in a% redundant familly of N*(J+1) atoms where N is the number of pixel and J% is the number of scales of the transforms.%%% For Scilab, we need to extend a little the available memory.extend_stack_size(4);%%% The invariant transform is obtained using the same function, by% activating the switch |options.ti=1|.options.ti = 1;MW = perform_wavelet_transf(M0,Jmin,+1,options);%%% |MW(:,:,1)| corresponds to the low scale residual.% Each |MW(:,:,3*j+k+1)| for k=1:3 (orientation) corresponds to a scale of wavelet coefficient, and has% the same size as the original image.clf;i = 0;for j=1:2    for k=1:3        i = i+1;        imageplot(MW(:,:,i+1), strcat(['Scale=' num2str(j) ' Orientation=' num2str(k)]), 2,3,i );    endend%% Translation Invariant Wavelet Denoising% Orthogonal wavelet denoising does not performs very well because of its% lack of translation invariance.% A much better result is obtained by not sub-sampling the wavelet% transform, which leads to a redundant tight-frame.%%% First we compute the translation invariant wavelet transformoptions.ti = 1;MW = perform_wavelet_transf(M,Jmin,+1,options);%%% Then we threshold the set of coefficients.T = 3.5*sigma;MWT = perform_thresholding(MW,T,'hard');%%% We can display some wavelets coefficientsJ = size(MW,3)-5;clf;imageplot(MW(:,:,J), 'Noisy coefficients', 1,2,1);imageplot(MWT(:,:,J), 'Thresholded coefficients', 1,2,2);%%% We can now reconstructMti = perform_wavelet_transf(MWT,Jmin,-1,options);% displayclf;imageplot(clamp(Msoft), strcat(['Soft orthogonal, SNR=' num2str(snr(M0,Msoft))]), 1,2,1);imageplot(clamp(Mti), strcat(['Hard invariant, SNR=' num2str(snr(M0,Mti))]), 1,2,2);%%% _Exercice 4:_ (the solution is <../private/denoising_wavelet/exo4.m exo4.m>)% Determine the best threshold |T| for both hard and soft thresholding,% but now in the translation invariant case. What can you conclude ?exo4;%% Wavelet Block Thresholding% Wavelets 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.%%% You can perform the block thresholding for an arbitrary block size.options.ti = 0;MW = perform_wavelet_transf(M,Jmin,+1,options);% soft block thresholdingT = 2.5*sigma/2;options.block_size = 4;MWT = perform_thresholding(MW,T,'block',options);% displayplot_wavelet(MWT,Jmin);%%% You can reconstruct the image. Test with several values for |T| in order to determine the best threshold. Mblock = perform_wavelet_transf(MWT,Jmin,-1,options);% displayclf;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 5:_ (the solution is <../private/denoising_wavelet/exo5.m exo5.m>)% Try block thresholding for a variety of block size and determine% the best SNR.exo5;%%% Block thresholding can also be applied to a translation invariant wavelet% transform. It gives state of the art denoising results.% transformoptions.ti = 1;MW = perform_wavelet_transf(M,Jmin,+1,options);% thresholdT = 2.5*sigma/2;options.block_size = 5;MWT = perform_thresholding(MW,T,'block',options);% transform backMblockti = perform_wavelet_transf(MWT,Jmin,-1,options);% displayclf;imageplot(clamp(Mti), strcat(['Hard TI, SNR=' num2str(snr(M0,Mti))]), 1,2,1);imageplot(clamp(Mblockti), strcat(['Block TI, SNR=' num2str(snr(M0,Mblockti))]), 1,2,2);##### SOURCE END #####-->   </body></html>