function [Arec,bitstream,PSNR,MSE,D,Drec,s,p_stream]=spiht_wpackets(A,bpp,wavelet,no_decomp,pkt_depth,dec_type,varargin)
%[Arec,bitstream,PSNR,MSE,D,Drec,s,p_stream]=spiht_wpackets(A,bpp,wavelet,no_decomp,pkt_depth,dec_type,varargin)
%Version: 3.13, Date: 2006/04/10, author: Nikola Sprljan
%SPIHT image compression using Wavelet Packets (WP) decomposition
%
%Input: 
% A - array containing the original image or its filename
% bpp - vector of target decoding bitrates (bpp=bits per pixel) 
% wavelet - wavelet identification string
% no_decomp - [optional, default = max] the number of levels of dyadic 
%             decomposition. Default value is the maximum possible, e.g. down
%             to the LL subband of size 1x1 if the dimensions of an image are
%             powers of 2.
% pkt_depth - [optional, default = 0] maximum number of additional decomposition
%             levels for subbands; "packet decomposition depth"
% dec_type - type of the wavelet packets deconposition
%            if (dec_type == 'full') builds the full packets tree
%            if (dec_type == 'greedy') sub-optimal in terms of finding the 
%            minimal entropy, but faster
% varargin - parameters for entropy computation (see in decomp_packets.m)
%
%Output: 
% Arec - reconstructed image (if more than one decode bitrate is specified than
%        Arec contains the result of decoding at the highest bitrate)
% bitstream - output bitstream and side information; is of size (3,): 
%             (1,bitstream size) the bitstream 
%             (2,bitstream size) pass number of the SPIHT algorithm  
%             (3,bitstream size) bit class: 
%               {0 - header,1 - position bit,2 - sign bit,3 - refine bit}
% PSNR - PSNR of the reconstructed images (for all bpps)
% MSE - MSE of the reconstructed images
% D - wavelet coefficients after decomposition
% Drec - reconstructed wavelet coefficients 
% s - structure containing info on parent-children relationship between subbands
%     given by wavelet packets decomposition
% p_stream - stream of bits representing information on splitting decisions of
%            wavelet packets decomposition 
%
%Note: 
% spiht_wpackets.m calls the following external function (in this order):          
%   1. decomp_packets2D.m (wavelet packets decomposition),
%   2. pdf_opt.m (probability density function optimisation of the quantised 
%                 wavelet coefficient values; Laplace pdf),
%   3. recon_packets2D.m (wavelet packets reconstruction).    
% Output file SPIHTlog.txt contains statistics collected during the encoding 
% process.
%
%Uses: 
% pdf_opt.m 
% load_wavelet.m (Wavelet Toolbox)
% decomp_packets2D.m (Wavelet Toolbox)
% recon_packets2D.m (Wavelet Toolbox)
% iq_measures.m (Quality Assessment Toolbox)
%
%Example:
% [Arec,bitstream,PSNR]=spiht_wpackets('Lena512.png',1,'CDF_9x7',6); %as original binary-coded SPIHT
% [Arec,bitstream,PSNR]=spiht_wpackets(A,1,'haar');
% [Arec,bitstream,PSNR]=spiht_wpackets('Lena.bmp',[0.1 0.2 0.5 1],'CDF_9x7',5,4,'full','shannon');
% [Arec,bitstream,PSNR,MSE,D,Drec,s,p_stream]=spiht_wpackets('lena256.png',0.1,'CDF_9x7',5,4,'greedy','shannon',1.5);

if isstr(A)
    A=imread(A);
end;
if (ndims(A)>2) 
    error('Only indexed images are supported!');
end;
if (exist('decomp_packets2D.m') ~= 2)
    disp(sprintf('Warning: The function decomp_packets2D.m is not on the path, wavelet decomposition cannot be performed!'));
    disp(sprintf('Function decomp_packets2D.m is a part of Wavelet Toolbox'));
    return;
end;
if (exist('iq_measures.m') ~= 2)
    disp(sprintf('Warning: The function iq_measures.m is not on the path, the PSNR computation will not be performed!'));
    disp(sprintf('Function iq_measures.m is a part of Quality Assessment Toolbox'));
end;

if nargin<6
    dec_type='';
    if nargin<5
        pkt_depth=0;
        if nargin<4
            [Drows,Dcols]=size(A);
            rpow2=length(find(factor(Drows)==2));
            cpow2=length(find(factor(Dcols)==2));
            no_decomp=min(rpow2,cpow2);
            if 2^no_decomp==min(Drows,Dcols)
                no_decomp=no_decomp-1;   
            end;    
        end;
    end;
end;
fprintf('%d decompositions, %d packet depth, %s wavelet, highest bitrate %f bpp\n',...
    no_decomp,pkt_depth,wavelet,max(bpp));     

param = struct('N',no_decomp,'pdep',pkt_depth,'wvf',wavelet,'dec',dec_type);

ent_param = struct('ent',[],'opt',[]);
%loading wavelet
%load the wavelet here
if (isstr(wavelet))
    if (exist('load_wavelet.m') ~= 2)
        disp(sprintf('Warning: The function load_wavelet.m is not on the path, the specified wavelet cannot be loaded!'));
        disp(sprintf('Function load_wavelet.m is a part of Wavelet Toolbox'));
        return;
    end;
    param.wvf = load_wavelet(wavelet);
else
    param.wvf = wavelet;
end;

if nargin>6
    ent_param.ent=varargin{1};
    if nargin==8
        ent_param.opt=varargin{2};
    end;
else
    ent_param.ent='shannon';
end; 
param.meanvalue=round(mean(A(:)));
Am=double(A)-param.meanvalue;
%Am=A;param.meanvalue=0; %mean value shifting as in original SPIHT
tic;
[D,p_stream,s,E]=decomp_packets2D(Am,param,ent_param);
%global Dcmap; %enable this to load a decomposition matrix from the main environment
%D = Dcmap;

%mean value shifting as in original SPIHT /begin
%DLLr=size(A,1)/2^param.N;DLLc=size(A,2)/2^param.N;
%Dmean=mean(mean(D(1:DLLr,1:DLLc)));
%shift=0;
%while Dmean>1024 Dmean=Dmean/4;shift=shift+1;end;
%param.Dmean=ceil(Dmean);param.shift=shift;
%D(1:DLLr,1:DLLc)=D(1:DLLr,1:DLLc)-param.Dmean*4^param.shift;
%mean value shifting as in original SPIHT /end

fprintf('Image transformed in %.2f seconds\n',toc);
tic;
%fprintf('Entropy-  %f\n',E);
[M,s]=maximum_magnitudes(D,s);
fprintf('Maximum magnitudes computed in %.2f seconds\n',toc);
numbits=round(prod(size(A)).*bpp);
%SPIHT algorithm, coding&decoding at once (without output to a file!)
if pkt_depth>0
    %4 bits for depth of packet decomposition, 16 for the length of p_stream; not implemented
    header_size=51+4+16+length(p_stream);  
else
    header_size=51; %51 bits as in original SPIHT
end;
header=zeros(1,header_size);
header(1:8)=str2num(dec2bin(param.meanvalue,8)')';
header(9:22)=str2num(dec2bin(size(A,1),14)')';
header(23:36)=str2num(dec2bin(size(A,1),14)')';
header(37:41)=str2num(dec2bin(floor(log2(max(max(M(:,:,1))))),5)')';
header(42:45)=str2num(dec2bin(no_decomp,4)')';
if pkt_depth>0
    header(46)=1; %next five bits for the type of wavelet
    header(72:header_size)=p_stream;   
end;
fprintf('\nSPIHT algorithm...\n');
tic;
[Drec,Arec,bitstream,PSNR,MSE]=spihtalg(A,D,s,M,numbits,header_size,param,p_stream);%

bitstream(1,1:header_size)=header;
fprintf(' - completed in %.2f seconds\n',toc); 

function [Drec,Arec,bitstream,PSNR,MSE]=spihtalg(A,D,s,M,numbits,header_size,param,p_stream)
%Input: 
% A - array containing the original image or its filename
% D - wavelet coefficients of the original image
% s - structure containing info on parent-children relationship between
%     subbands (see in decomp_packets2D.m and in the function
%     maximum_magnitudes defined below) 
% M - sorted coefficients trees (see in maximum_magnitudes below)
% numbits - vector containing bit budgets for each of the specified bits
%           per pixel values
% header_size - header size in bits to be taken into account
% param - structure specifying the decomposition parametres (see
%         decomp_packets2D.m)
% p_stream - information necessary for reconstruction of wavelet packets
%            subbands (see decomp_packets2D.m)

N=param.N;
[Drows,Dcols]=size(D);
%fileout=fopen('fileout.bin','w');%f
maxcoeff=max(max(abs(D))); %maximum absolute coefficient 
if maxcoeff>0 
    nbit=floor(log2(maxcoeff)); %bits for first threshold
else %in case that image contains only DC component (which is previously substracted) 
    quantstep=0;
    return %the whole coding process is skipped
end;
%Energy=sum(sum(D.^2));
%RecEnergy=0;
%initialize all lists & variables
[LIP,endLIP,LIS,endLIS,LSP,endLSP]=initialize_lists(D,N);
LISnew=zeros(size(LIS));
LIStemp=zeros(64,3);
oldendLSP=0;
duzina=numel(D);
sizrow=Drows/2^N;
sizcol=Dcols/2^N;
Dqrows=Drows/4;
Dqcols=Dcols/4;
MagnL=M(:,:,1);
MagnD=M(:,:,2);
L_S=M(:,:,3); %matrix that links coefficient and the index of its subband
Drec=zeros(Drows,Dcols); %initialize matrix of reconstructed coefficients
bitstream=zeros(3,numbits(end),'uint8');
quantstep=2^nbit; %first threshold (threshold is equal to quantization step)
Dtype=1; %indicates D type of LIS element
pass=1; %first coding pass
%**STATS**
nowcd=fileparts(which('spiht_wpackets'));
namelog=fullfile(nowcd,'SPIHTlog.txt');
fid=fopen(namelog,'w');
count1=header_size; %bits reserved for header
[stats,fid]=show_stats(fid,'init',pass,count1);
cntbitout=1;
while numbits(cntbitout)<=header_size
    [MSE(cntbitout),PSNR(cntbitout)]=iq_measures(A,127*ones(size(A)));
    cntbitout=cntbitout+1;
end;
%2.SORTING PASS 
while 1
    %fprintf(fid,'Threshold = %.2f\n   ',quantstep);  %**STATS**
    fprintf(fid,'%2d.(%8d)   ',pass,count1+1);  %**STATS**
    %fprintf(fid,'%2d LIP\n',pass); %NS!
    %2.1. testing the significance of elements in LIP
    lipcnt=0;
    cntbit=count1;
    init=quantstep*3/2; 
    %RecEadd=init^2;
    for mlip=1:endLIP %LIP loop
        i=LIP(mlip,1);
        j=LIP(mlip,2);
        count1=count1+1;
        bitstream(1,count1)=(abs(D(i,j))>=quantstep);
        bitstream(2,count1)=pass;
        bitstream(3,count1)=1; %position bit
        if bitstream(1,count1) %significant goes to LSP
            %fprintf(fid,'%5.2f\n',D(i,j));  %**STATS**
            count1=count1+1;
            bitstream(1,count1)=(sign(D(i,j))==1);
            bitstream(2,count1)=pass;
            bitstream(3,count1)=2; %sign bit
            Drec(i,j)=sign(double(bitstream(1,count1))-0.5)*init;
            %RecEnergy=RecEnergy+RecEadd;
            endLSP=endLSP+1;
            LSP(endLSP,:)=[i j];
            %fwrite(fileout,1,'ubit1');%f
            %fwrite(fileout,sign(D(i,j)),'ubit1');%f
        else %otherwise stays in LIP
            lipcnt=lipcnt+1;
            LIP(lipcnt,:)=LIP(mlip,:);
            %fwrite(fileout,0,'ubit1');%f
        end;
        if count1+1>=numbits(cntbitout) %this limit is equivalent in terms of distortion (1 addidional bit wouldn't help)
            [Arec,PSNR(cntbitout),MSE(cntbitout)]=reconstruct(A,Drec,cntbitout,count1,quantstep,pass,param,p_stream);
            cntbitout=cntbitout+1;
        end;
        if cntbitout>length(numbits) %| abs(RecEnergy-Energy)/Energy<10^-7
            [stats,fid]=show_stats(fid,'LIP',pass,count1,stats,endLSP); %**STATS**
            [stats,fid]=show_stats(fid,'LIS',pass,count1,stats,endLSP);
            [stats,fid]=show_stats(fid,'LSP',pass,stats(pass,1).LSP-1,stats,endLSP);
            [stats,fid]=show_stats(fid,'total',pass,count1,stats,endLSP);
            fprintf('\n...LIP loop out');
            %fclose(fileout);
            return
        end;    
    end; %end of LIP loop
    endLIP=lipcnt;
    [stats,fid]=show_stats(fid,'LIP',pass,count1,stats,endLSP); %**STATS**
    %fprintf(fid,'%2d LIS\n',pass); %NS!
    %2.2. testing the significance of elements in LIS 
    mlis=0;
    mlistemp=0;
    endLIStemp=0;
    liscnt=0; 
    while mlis<endLIS
        if endLIStemp
            mlistemp=mlistemp+1;
            LISel=LIStemp(mlistemp,:);
            if mlistemp==endLIStemp
                mlistemp=0;
                endLIStemp=0;   
            end;    
        else
            mlis=mlis+1;
            LISel=LIS(mlis,:);
        end;
        i=LISel(1);
        j=LISel(2);
        subband_index=L_S(i,j);
        row=s(subband_index).band_abs(2);
        col=s(subband_index).band_abs(1);
        scalediff=s(subband_index).scalediff;
        rowrel=i-row;
        colrel=j-col;
        id=row+s(subband_index).addrows+(rowrel-1)*2.^scalediff+1;
        jd=col+s(subband_index).addcols+(colrel-1)*2.^scalediff+1;
        %2.2.1. for L type elements
        if ~LISel(3)
            polje=MagnL(id+(jd-1)*Drows); %children
            count1=count1+1;