function compresn(y,threshold);

%       Written by Colm Mulcahy (colm@spelman.edu), last
%       updated 31st Dec 1996.
%
%       Normalized Haar wavelet compression
%       for  256  by  256  grayscale images.
%
%       Usage:  COMPRESN(Y,THRESHOLD)  where  Y  is a
%       256  by  256  image matrix, and  THRESHOLD is the
%       "percentage of details suppressed" in the compression.
%
%       Compresses in three steps:
%
%       Matrix  Y  is wavelet transformed (8  times to the
%       rows, then  8  times to the resulting columns),
%       to yield a transformed matrix  TY;
%
%       All entries of  TY  whose magnitude do not exceed
%       max(max(abs(TY)))*threshold/100  are set to zero,
%       this doctored matrix is called  DY;
%
%       The inverse wavelet transformation is applied to
%       DY  to yield the compressed image matrix  CY.
%
%       Output is pictures of  Y  and  CY,
%       and the so-called "compression rate"  cr:
%       the reciprocal of the percentage of non-zero
%       entries in the image after stage two above.
%
%	If only 2% of these entries are non-zero, we
%       have  cr=50  and we say we have a  50:1  rate.
%       This conveniently ignores the fact that the
%       original image - or the pre-compression wavelet
%       transformed version - may have already had a
%       significant number of zero elements.
%

low=min(min(y));  high=max(max(y));     % NEEDED LATER FOR DISPLAYING
 
figure(1)
imagesc(y), colormap(gray), axis off
title('original image')
 
w=wavmat;
 
ty = w'*y*w;                            % THE WAVELET TRANSFORMATION
 
tty=ty;                                 % A TRICK TO ELIMNATE THE OVERALL
tty(1,1)=0;                             % AVERAGE FROM CONSIDERATION WHEN
dead=max(max(abs(tty)))*threshold/100;  % DECIDING JUST HOW TO THRESHOLD
clear tty;
 
dy=ty;
 
index=find(abs(dy)<=dead);
dy(index)=zeros(size(index));           % SETTING LOTS OF ELEMENTS TO ZERO
 
cy = full(w*sparse(dy)*w');             % THE INVERSE WAVELET TRANSFORMATION
 
density = nnz(dy);                      % ENTRIES USED OUT OF 256^2 = 65536
disp(['Wavelet transformed and doctored matrix uses '])
disp([ num2str(density) ' entries out of 256^2 = 65536,'])
disp(['thus is ' num2str(100*density/65536) '% dense, and we get a '])
disp(['compression ratio of ' num2str(65536/nnz(dy)) ' to 1'])  % "COMPRESSION RATIO" = 65536/DENSITY
 
figure(2)
imagesc(cy,[low high]), colormap(gray), axis off
title('compresned image')