% Demo with the 2D WP.

% Read image;
img = ReadImage('Lenna');
[N,J] = quadlength(img);
subplot(311);imagesc(img);axis image;axis off;colormap('gray')
set(gca,'FontSize',14);
title('(a)');

% Deepest decomposition scale.
deepest = 5;
coarsest = J - deepest;

% Generate Symmlet 4 CMF Filter.
qmf = MakeONFilter('Symmlet',4);

% Compute stat tree.
stats = Calc2dStatTree('WP',img - mean(img(:)),deepest,qmf,'Entropy',[]);

% Coifman-Wickerhauser best orthobasis algorithm.
bob   = Best2dBasis(stats,deepest);
wpc   = FPT2_WP(bob,img,qmf);

% Sort and threshold the WP best orthobasis coeffs.
swpc = sort(abs(wpc(:)),1,'descend');
thd  = swpc(floor(N*N/10));
ll   = wpc(1:coarsest,1:coarsest);
wpct = wpc .* (abs(wpc) > thd);
wpct(1:coarsest,1:coarsest) = ll;

% Display retained coeffs in the best orthobasis with the corresponding partition.
subplot(346);imagesc(wpct~=0);axis image;axis off;
ax = axis;hold on
Plot2dPartition(bob,'b',ax,deepest);
set(gca,'FontSize',14);
title('(b)');

% Reconstruct and display.
imwpr = IPT2_WP(bob,wpct,qmf);
subplot(347);imagesc(imwpr);axis image;axis off;
set(gca,'FontSize',14);
title('(c)');
fprintf('Relative reconstruction error WP BOB: %.3f \n',norm(img(:)-imwpr(:))/norm(img(:)));

% Compute orthogonal DWT of image.
wc = FWT2_PO(img,coarsest,qmf);
% Sort, threshold and reconstruct.
swc = sort(abs(wc(:)),1,'descend');
thd = swc(floor(N*N/10));
ll  = wc(1:coarsest,1:coarsest);
wct = wc .* (abs(wc) > thd);
wct(1:coarsest,1:coarsest) = ll;

% Display retained DWT coeffs.
subplot(3,4,10);imagesc(wct~=0);axis image;axis off;
set(gca,'FontSize',14);
title('(d)');

% Reconstruct and display. 
imdwtr = IWT2_PO(wct,coarsest,qmf);
subplot(3,4,11);imagesc(imdwtr);axis image;axis off;
set(gca,'FontSize',14);
title('(e)');
fprintf('Relative reconstruction error DWT: %.3f \n',norm(img(:)-imdwtr(:))/norm(img(:)));

saveas(gcf,'wpdemo.eps','epsc');