% Demo with the 2D biorthogonal DWT.

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

% Coarsest decomposition scale.
coarsest = 3;

% Generate 7/9 Biorthonormal CMF Filter Pairs.
[qmf,dqmf] = MakeBSFilter('CDF',[4 4]);

% Compute biorthogonal DWT of image.
wc = FWT2_PB(img,coarsest,qmf,dqmf);

% Plot DWT NLA error.
swc   = sort(abs(wc(:)),1,'descend');
wcerr = flipud(cumsum(flipud(swc).^2));
subplot(222);plot(log2(1:N*N),log2(wcerr));set(gca,'XTick',[0:2:log2(N*N)]);
axis([4 14 log2(wcerr(2^14)) log2(wcerr(2^4))]);
set(gca,'FontSize',14);
xlabel('log_2(Number of retained DWT coeff)');ylabel('log_2(m-term approximation error)');
title('(b)');

% Keep only 10% of coefficients.
thd = swc(floor(N*N/10));
ll  = wc(1:coarsest,1:coarsest);
wct = wc .* (abs(wc) >= thd);
wct(1:coarsest,1:coarsest) = ll;
subplot(223);imagesc(wct~=0);axis image;axis off
set(gca,'FontSize',14);
title('(c)');

% Reconstruct from remaining 10% DWT coefficients.
imr = IWT2_PB(wct,coarsest,qmf,dqmf);
subplot(224);imagesc(imr);axis image;axis off
set(gca,'FontSize',14);
title('(d)');

fprintf('Relative reconstruction error: %.4f \n',norm(img(:)-imr(:))/norm(img(:)));

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