% This is an example of using the 2-d wavelet decomposition
% software in the Matlab wavelet toolbox.

% Some of the commands and comments are extracted
% verbatim from Mathworks' documentation.

% Code and documentation modifications by DPO 02/2010
% Image elk.jpg derived from a photo by Timothy O'Leary

x8 = imread('elk.jpg');
X = double(x8);
colormap('gray');
figure(1)
imagesc(X);
title('Original image, in double prec.')
figure(2)

% Perform a level 2 decomposition of the image 
% using the bior3.7 wavelet

[C,S] = wavedec2(X,2,'bior3.7');


% The coefficients of all the components of a 
% second-level decomposition (that is, the 
% second-level approximation and the first 
% two levels of detail) are returned 
% concatenated into one vector, C. 
% Argument S is a bookkeeping matrix that 
% keeps track of the sizes of each component.

% Extract approximation and detail coefficients.

% To extract the level 2 approximation coefficients from C:
%
% cA2 = appcoef2(C,S,'bior3.7',2);

% To extract the first- and second-level detail 
% coefficients from C:
%
% [cH2,cV2,cD2] = detcoef2('all',C,S,2); 
% [cH1,cV1,cD1] = detcoef2('all',C,S,1); 

% The first argument ('h', 'v', 'd', or 'all') 
% determines the type of detail (horizontal, 
% vertical, diagonal) extracted, and the last 
% argument determines the level.

% Reconstruct the Level 2 approximation.

A1 = wrcoef2('a',C,S,'bior3.7',1);
A2 = wrcoef2('a',C,S,'bior3.7',2);

% Reconstruct the level 1 and 2 details.

            H1 = wrcoef2('h',C,S,'bior3.7',1); 
            V1 = wrcoef2('v',C,S,'bior3.7',1); 
            D1 = wrcoef2('d',C,S,'bior3.7',1); 
            H2 = wrcoef2('h',C,S,'bior3.7',2);
            V2 = wrcoef2('v',C,S,'bior3.7',2); 
            D2 = wrcoef2('d',C,S,'bior3.7',2);

disp(sprintf('Norm of difference between A1+H1+V1+D1 and original image is %e', norm(X-A1-H1-V1-D1,'fro') ))

% Display the results of the level-2 decomposition.

% Advice from Mathworks:
% With all the details involved in a 
% multilevel image decomposition, it makes sense 
% to import the decomposition into the Wavelet 2-D 
% graphical tool in order to more easily display it. 
% For information on how to do this, see Loading 
% Decompositions in 
% http://www.mathworks.com/access/helpdesk/help/toolbox/wavelet/ch02_us6.html#1010498. 

% Originally the display commands were i
% image(wcodemat(..,192)) instead of imagesc(..)

colormap('gray');
subplot(2,4,1);imagesc(A1);         
title('Approximation A1')
subplot(2,4,2);imagesc(H1);         
title('Horizontal Detail H1')
subplot(2,4,3);imagesc(V1);         
title('Vertical Detail V1')
subplot(2,4,4);imagesc(D1);         
title('Diagonal Detail D1')
subplot(2,4,5);imagesc(A2);         
title('Approximation A2')
subplot(2,4,6);imagesc(H2);         
title('Horizontal Detail H2')
subplot(2,4,7);imagesc(V2);         
title('Vertical Detail V2')
subplot(2,4,8);imagesc(D2);         
title('Diagonal Detail D2')

% Added example:
% Now construct and display a 5-level approximation
% and compare with the true one.

[C,S] = wavedec2(X,5,'bior3.7');

clim = [0 256];
normX = norm(X);

figure
title('5 approximations')
subplot(2,3,1);imagesc(X,clim);         
colormap('gray')
    title('Original Image')

disp(' ')
disp('Wavelet          Relative error (in percent)')
disp('Approximation    Frobenius          Max        Median')
for j=1:5,
    A.approx{j} = wrcoef2('a',C,S,'bior3.7',j);
    pixerr = X - A.approx{j};
    pixrel = abs(pixerr ./ X);
    subplot(2,3,j+1);imagesc(A.approx{j},clim);         
    colormap('gray')
    title(sprintf('Level %d approximation',j))
    disp(sprintf('Level %d     %12.1f %16.1f %10.1f',...
        j, 100*norm(pixerr,'fro')/normX, ...
           100*max(max(pixrel)), ...
           100*median(pixrel(:)) ))
end
disp(' ')

[m,n] = size(X);
mw = length(C);

disp(sprintf('There are %d elements in the original image.', m*n))

disp(sprintf('There are %d elements in the wavelet image.', mw))

disp(' ')
disp('Pausing. Strike any key to continue.')
pause

factor = 50;
Ctrunc = C.* (abs(C)> max(abs(C))/factor); 

disp(' ')
disp('But in the wavelet images, only')
disp(sprintf('     %d of these have magnitudes within a factor',...
       sum(Ctrunc ~= 0) ))
disp(sprintf('     of %d of the largest magnitude element.',factor))
disp('We will form images using just these large magnitude elements.')

figure
title('5 approximations')
subplot(2,3,1);imagesc(X,clim);
colormap('gray')
    title('Original Image')

disp(' ')
disp('Trunc. Wavelet   Relative error (in percent)')
disp('Approximation    Frobenius          Max        Median')
for j=1:5,
    A.approx{j} = wrcoef2('a',Ctrunc,S,'bior3.7',j);
    pixerr = X - A.approx{j};
    pixrel = abs(pixerr ./ X);
    subplot(2,3,j+1);imagesc(A.approx{j},clim);
    colormap('gray')
    title(sprintf('Level %d trunc. approx.',j))
    disp(sprintf('Level %d     %12.1f %16.1f %10.1f',...
        j, 100*norm(pixerr,'fro')/normX, ...
           100*max(max(pixrel)), ...
           100*median(pixrel(:)) ))
end
disp(' ')