function func_SPIHT_Main
% Matlab implementation of SPIHT (without Arithmatic coding stage)
%
% Main function 
%
% input:    Orig_I : the original image.
%           rate : bits per pixel
% output:   img_spiht 
%
% Jing Tian
% Contact me : scuteejtian@hotmail.com
% This program is part of my undergraduate project in GuangZhou, P. R. China.
% April - July 1999

%-----------   Input   ----------------
%I=imread('3668.tif');I1=I(1:512,1:512);
%X=mat2gray(I1);
load woman
Orig_I = X;
rate =0.5;

%-----------   Pre-processing   ----------------

OrigSize = size(Orig_I, 1);
max_bits = floor(rate * OrigSize^2);

OutSize = OrigSize;
image_spiht = zeros(size(Orig_I));
% "image " is the input of codec
[nRow, nColumn] = size(Orig_I);


%-----------   Wavelet Decomposition   ----------------
n = size(Orig_I,1);
n_log = log2(n); 
level = n_log;
% wavelet decomposition level can be defined by users manually.

type = 'db4';
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(type);

[I_W, S] = func_DWT(Orig_I,3, Lo_D, Hi_D);

%-----------   Coding   ----------------
img_enc = func_SPIHT_Enc(I_W, max_bits, nRow*nColumn,5);  
L=length(img_enc);
%-----------   Decoding   ----------------
img_dec = func_SPIHT_Dec(img_enc);

%-----------   Wavelet Reconstruction   ----------------
img_spiht = func_InvDWT(img_dec, S, Lo_R, Hi_R, 3);
X0=img_spiht;
%-----------   PSNR analysis   ----------------
Q = 255;
%r=(16*256*256)/L
MSE = sum(sum((img_spiht - Orig_I).^2))/( nRow * nColumn);
psnr = 10*log10(Q*Q/MSE)
%figure(1);imshow(X)
figure(2);subplot(221);image(X0);colormap(map);