% 2D image analysis.

% 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);

close all
clf
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)
pause
close(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 branches
ra2 = 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);

clf
subplot(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');

clf
subplot(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')