function WLena()
close all
% wavelet analysis of Lena
Lenna = imread('Lenna.png');
lena = sum(Lenna,3);
size(lena)
cmap = 'gray'; % choose colormap gray
[n,m] = size(lena); % in lean, n = m = 512
figure(1); clf;
colormap(cmap);
imagesc(lena);
title('The original image');
colorbar;
wav = 'db1'; % choose mother wavelet

flag = 4;
%% compute single-level wavelet decomposition
if flag == 1
    [cA1,cH1,cV1,cD1] = dwt2(lena,wav); % compute wavelet coedfficients
    % map the coefficients onto n-by-n plane
    A1 = upcoef2('a',cA1,wav,1);
    H1 = upcoef2('h',cH1,wav,1);
    V1 = upcoef2('v',cV1,wav,1);
    D1 = upcoef2('d',cD1,wav,1);
    figure;
    imagesc(A1)
    colormap(cmap)
    title('Averaged');
    colorbar
    figure;
    colormap(cmap)
    imagesc(H1)
    title('Horizontal differences');
    colorbar;
    figure;
    colormap(cmap)
    imagesc(V1)
    title('Vertical differences');
    colorbar;
    figure;
    colormap(cmap)
    imagesc(D1)
    title('Diagonal differences');
    colorbar;
    % restore Lena using the inverse wavelet transform
    Li = idwt2(cA1,cH1,cV1,cD1,wav);
    fprintf('The max difference between the original and the reconstructed image = %d\n',...
        max(max(abs(lena - Li))));
end
%% compute wavelet decomposition
if flag == 2
    J = log2(min(n,m)) % the maximal level of decomposition
    [C,S] = wavedec2(lena,J,wav);
    figure; hold on; 
    colormap gray
    for j = 1 : J
        [cH,cV,cD] = detcoef2('all',C,S,j);
        x1 = [0 : 2^(J - j) - 1];
        x2 = [2^(J - j) : 2^(J - j + 1) - 1];
        imagesc(x2,x1,log10(max(abs(cH),1)));
        imagesc(x1,x2,log10(max(abs(cV),1)));
        imagesc(x2,x2,log10(max(abs(cD),1)));
    end
    axis ij  
    axis tight
    colorbar
    % reconstruct Lena 
    Li = waverec2(C,S,wav);
    fprintf('The max difference between the original and the reconstructed image = %d\n',...
        max(max(abs(lena - Li))));
end
%% 
if flag == 3
    J = log2(min(n,m)); % the maximal level of decomposition
    [C,S] = wavedec2(lena,J,wav);
    figure; hold on;
    k = 1;
    for j = J : -1 : 2
        subplot(2,4,k);
        A = appcoef2(C,S,wav,j);
        colormap(cmap)
        imagesc(A);
        k = k + 1;
    end
end
%% compress Lena
if flag == 4
    cfactor = 10;
    % cfactor = factor by which the default threshold is multiplied
    J = log2(min(n,m)); % the maximal level of decomposition
    [C,S] = wavedec2(lena,J,wav);
    % the larger cfactor is th more coefficients are zeroed
    [threshold,sorh,keepapp] = ddencmp('cmp','wv',lena);
    [lena1,CC,LC,P0,PL2] = wdencmp('gbl',C,S,wav,J,cfactor*threshold,sorh,keepapp);
    fprintf('P0 = %d\tPL2 = %d\n',P0,PL2);
    % P0  = percentage of zeroed coefficients
    % PL2 = percentage of preserved image energy
    figure;
    colormap(cmap)
    imagesc(lena1)
    title('Restored after compression')
end