Wavelet image analysis and compression

Written to not use wavelet toolbox, instead using functions
dwt2haar and idwt2haar.  Alternate toolbox commands shown in comments.

Contents

Load image

RGB = imread('ngc6543a.jpg');
image(RGB);
size(RGB)
ans =

   650   600     3

Crop/pad to 640 x 576 to allow 6-level binary subsampling

pic = RGB(11:650, 12:587,:);
image(pic);

maxlev = 6;


C = dwt2Haar(pic,maxlev);
% Wavelet toolbox: [C,S] = wavedec2(pic,maxlev,'haar');

Plot cumulative wavelet power fraction vs. fraction of wavelets retained

% Unroll all wavelet coefficients into a single vector and sort by size
np = length(C(:));
Csq = C(:).^2;
[mCsqsort,Isort] = sort(-Csq);  % Tricky way to sort wavelet coeffs from
                                % largest to smallest in magnitude.
respowfrac = 1 - cumsum(mCsqsort)/sum(mCsqsort); % Cum frac of Wavelet power
semilogy((1:np)/np,respowfrac)
xlabel('Fraction of wavelets retained')
ylabel('Fraction of power not captured')
% Shows 2% of wavelet coeffs capture all but 0.2% of power and 35% capture
% all the rest.

Plot image filtered by only retaining the largest wavelets

keepfrac = 0.02;  % Zero out all but largest keepfrac of wavelet coeffs
                  % Worth trying different values of this parameter.
Izero = Isort(round(keepfrac*np):np);
Cfilt = C(:);
Cfilt(Izero) = 0;
Cfilt = reshape(Cfilt,size(C));

picfilt = idwt2Haar(Cfilt,maxlev);
% Wavelet toolbox: picfilt = waverec2(Cfilt,S,'haar');
image(uint8(picfilt));