The Wavelet Transform - Sparse Image and Signal Processing - page 43

% 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 coefficients.
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)');

% Compute orthogonal DWT of image.
wc = FWT2_PO(img,coarsest,qmf);