%WAVELET BASED COMPRESSION  
               %****************************

%removes all variables, globals, functions and
%MEX links(MATLAB loads and runs a different entry point symbol for C or Fortran MEX-files)
clear all;
% CLOSE ALL  closes all the open figure windows.
close all;
%read the image
input_image1=imread('rice.tif');
%display input image
%add noise
input_image=imnoise(input_image1,'speckle',.01);
figure;
imshow(input_image);
%give the number of decomposition level which must be integer and should not exceed 3
n=input('enter the decomposition level');

%*****************************************************************************
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('haar');
% computes four filters associated with the orthogonal or biorthogonal
%  wavelet named in the string 'wname'. 
%    The four output filters are:
 %       LO_D, the decomposition low-pass filter
 %       HI_D, the decomposition high-pass filter
  %      LO_R, the reconstruction low-pass filter
   %     HI_R, the reconstruction high-pass filter
   % Available wavelet names 'wname' are:
   % Daubechies: 'db1' or 'haar', 'db2', ... ,'db45'
    %Coiflets  : 'coif1', ... ,  'coif5'
    %Symlets   : 'sym2' , ... ,  'sym8', ... ,'sym45'
    %Discrete Meyer wavelet: 'dmey'
%**********************************************************************    

%wavedec2    - Multi-level 2-D wavelet decomposition.
[c,s]=wavedec2(input_image,n,Lo_D,Hi_D);

% gives the wavelet decomposition of the matrix input_image at level n, using the
% wavelet named in string 'wname' or low pass and high pass
% Outputs are the decomposition vector C and the
%corresponding bookkeeping matrix S.
disp(' the decomposition vector Output is');
disp(c);
%*************************************************************************************
 
%Thresholds for wavelet 2-D using Birge-Massart strategy.
[thr,nkeep] = wdcbm2(c,s,1.5,3*prod(s(1,:)));

 % give level-dependent thresholds 'thr'and numbers of coefficients to be kept 'nkeep'
 % for compression. 'thr' is obtained using a wavelet coefficients 
 %selection rule based on Birge-Massart strategy.
  %disp('level-dependent thresholds');
 %disp(thr);
 %disp(' numbers of coefficients to be');
 %disp(nkeep);
%***************************************************************************
  %compression using wavelet packets.
 [compressed_image,TREED,comp_ratio,PERFL2] =WPDENCMP(thr,'s',n,'haar','threshold',5,1);
 disp('compression ratio in percentage');
 disp(comp_ratio);
 % returns a compressed version compressed_image of input
 % signal 'thr' (2-D) obtained by wavelet packet  coefficients thresholding.
  %  The additional output argument TREED is the
  %  wavelet packet best tree decomposition of compressed_image.
  %  PERFL2 and PERF0 are L^2 recovery and compression scores in percentages.
  
  %Multi-level 2-D wavelet reconstruction.
 re_ima1 = waverec2(c,s,'haar'); 
 re_ima=uint8(re_ima1);
  subplot(1,3,1);
 imshow(input_image);
 title('i/p image');
 subplot(1,3,2);
 imshow(compressed_image);
 title('compressed image');
 subplot(1,3,3);
 imshow(re_ima);
 title('reconstructed image');