function example4()% EXAMPLE4 2D image analysis. Demonstrates image compression using% biorthogonal 9/7 filters.% This file is part of the examples provided for the MIT course% 18.327 / 1.130 "Wavelets, Filter Banks and Applications"% (http://web.mit.edu/1.130).	% Author : Kevin Amaratunga, MIT% Version: 1.0	% Load a test image.  Matlab test images consist of a matrix, X,% color palette, map, which maps each value of the matrix to a% color.  Here, we will apply the Discrete Wavelet Transform to X.load woman2%load detfingr; X = X(1:200,51:250); image(X)colormap(map);axis image; set(gca,'XTick',[],'YTick',[]); title('Original')pause% We will use the 9/7 filters with symmetric extension at the% boundaries.dwtmode('sym')wname = 'bior4.4'% Plot the structure of a two stage filter bank.t = wtree(X,2,'bior4.4');plot(t)pauseclose(2)% Compute a 2-level decomposition of the image using the 9/7 filters.[wc,s] = wavedec2(X,2,wname);% Extract the level 1 coefficients.a1 = appcoef2(wc,s,wname,1);         h1 = detcoef2('h',wc,s,1);           v1 = detcoef2('v',wc,s,1);           d1 = detcoef2('d',wc,s,1);           % Extract the level 2 coefficients.a2 = appcoef2(wc,s,wname,2);h2 = detcoef2('h',wc,s,2);v2 = detcoef2('v',wc,s,2);d2 = detcoef2('d',wc,s,2);% Display the decomposition up to level 1 only.ncolors = size(map,1);              % Number of colors.sz = size(X);cod_a1 = wcodemat(a1,ncolors); cod_a1 = wkeep(cod_a1, sz/2);cod_h1 = wcodemat(h1,ncolors); cod_h1 = wkeep(cod_h1, sz/2);cod_v1 = wcodemat(v1,ncolors); cod_v1 = wkeep(cod_v1, sz/2);cod_d1 = wcodemat(d1,ncolors); cod_d1 = wkeep(cod_d1, sz/2);image([cod_a1,cod_h1;cod_v1,cod_d1]);axis image; set(gca,'XTick',[],'YTick',[]); title('Single stage decomposition')colormap(map)pause% Display the entire decomposition upto level 2.cod_a2 = wcodemat(a2,ncolors); cod_a2 = wkeep(cod_a2, sz/4);cod_h2 = wcodemat(h2,ncolors); cod_h2 = wkeep(cod_h2, sz/4);cod_v2 = wcodemat(v2,ncolors); cod_v2 = wkeep(cod_v2, sz/4);cod_d2 = wcodemat(d2,ncolors); cod_d2 = wkeep(cod_d2, sz/4);image([[cod_a2,cod_h2;cod_v2,cod_d2],cod_h1;cod_v1,cod_d1]);axis image; set(gca,'XTick',[],'YTick',[]); title('Two stage decomposition')colormap(map)pause% Here are the reconstructed branchesra2 = wrcoef2('a',wc,s,wname,2);rh2 = wrcoef2('h',wc,s,wname,2);rv2 = wrcoef2('v',wc,s,wname,2);rd2 = wrcoef2('d',wc,s,wname,2);ra1 = wrcoef2('a',wc,s,wname,1);rh1 = wrcoef2('h',wc,s,wname,1);rv1 = wrcoef2('v',wc,s,wname,1);rd1 = wrcoef2('d',wc,s,wname,1);cod_ra2 = wcodemat(ra2,ncolors);cod_rh2 = wcodemat(rh2,ncolors);cod_rv2 = wcodemat(rv2,ncolors);cod_rd2 = wcodemat(rd2,ncolors);cod_ra1 = wcodemat(ra1,ncolors);cod_rh1 = wcodemat(rh1,ncolors);cod_rv1 = wcodemat(rv1,ncolors);cod_rd1 = wcodemat(rd1,ncolors);subplot(3,4,1); image(X); axis image; set(gca,'XTick',[],'YTick',[]); title('Original')subplot(3,4,5); image(cod_ra1); axis image; set(gca,'XTick',[],'YTick',[]); title('ra1')subplot(3,4,6); image(cod_rh1); axis image; set(gca,'XTick',[],'YTick',[]); title('rh1')subplot(3,4,7); image(cod_rv1); axis image; set(gca,'XTick',[],'YTick',[]); title('rv1')subplot(3,4,8); image(cod_rd1); axis image; set(gca,'XTick',[],'YTick',[]); title('rd1')subplot(3,4,9); image(cod_ra2); axis image; set(gca,'XTick',[],'YTick',[]); title('ra2')subplot(3,4,10); image(cod_rh2); axis image; set(gca,'XTick',[],'YTick',[]); title('rh2')subplot(3,4,11); image(cod_rv2); axis image; set(gca,'XTick',[],'YTick',[]); title('rv2')subplot(3,4,12); image(cod_rd2); axis image; set(gca,'XTick',[],'YTick',[]); title('rd2')pause% Adding together the reconstructed average at level 2 and all of% the reconstructed details gives the full reconstructed image.Xhat = ra2 + rh2 + rv2 + rd2 + rh1 + rv1 + rd1;sprintf('Reconstruction error (using wrcoef2) = %g', max(max(abs(X-Xhat))))% Another way to reconstruct the image.XXhat = waverec2(wc,s,wname);sprintf('Reconstruction error (using waverec2) = %g', max(max(abs(X-XXhat))))% Compression can be accomplished by applying a threshold to the% wavelet coefficients.  wdencmp is the function that does this.% 'h' means use hard thresholding. Last argument = 1 means do not% threshold the approximation coefficients.%    perfL2 = energy recovery = 100 * ||wc_comp||^2 / ||wc||^2.%             ||.|| is the L2 vector norm.%    perf0 = compression performance = Percentage of zeros in wc_comp.thr = 20;                                                    [X_comp,wc_comp,s_comp,perf0,perfL2] = wdencmp('gbl',wc,s,wname,2,thr,'h',1);clfsubplot(1,2,1); image(X); axis image; set(gca,'XTick',[],'YTick',[]);title('Original')cod_X_comp = wcodemat(X_comp,ncolors);subplot(1,2,2); image(cod_X_comp); axis image; set(gca,'XTick',[],'YTick',[]);title('Compressed using global hard threshold')xlabel(sprintf('Energy retained = %2.1f%% \nNull coefficients = %2.1f%%',perfL2,perf0))pause% Better compression can be often be obtained if different thresholds% are allowed for different subbands.thr_h = [21 17];        % horizontal thresholds.              thr_d = [23 19];        % diagonal thresholds.                thr_v = [21 17];        % vertical thresholds.  thr = [thr_h; thr_d; thr_v];[X_comp,wc_comp,s_comp,perf0,perfL2] = wdencmp('lvd',X,wname,2,thr,'h');clfsubplot(1,2,1); image(X); axis image; set(gca,'XTick',[],'YTick',[]);title('Original')cod_X_comp = wcodemat(X_comp,ncolors);subplot(1,2,2); image(cod_X_comp); axis image; set(gca,'XTick',[],'YTick',[]);title('Compressed using variable hard thresholds')xlabel(sprintf('Energy retained = %2.1f%% \nNull coefficients = %2.1f%%',perfL2,perf0))% Return to default settings.dwtmode('zpd')