% Code started 31th Oct 2003, Ex taken from J.M.Shapiro's paper on EZW
% This code is not intented to be optimized, just to make things work
% and to understand the principles of SPHIT
% Reference : A New, Fast, Efficient Image Codec using Set Partitioning of Hierarchical Trees
%                                           :- Amir. Said, W. A Pearlman
clear,clc
load filename output
format long
load indices
% #####################
choice=1;
% #####################
if choice==1
    load lena.mat
    x=x(80:80+127,70:70+127);
elseif choice==2
    load kitten
    x=kitten(50:50+127,50:50+127);
elseif choice==3
    load penny
    x=P;
elseif choice==4
    x=imread('barbara.tif');
    x=im2double(x);
    x=220*x./x(1);
    x=x(50:50+127,300:300+127);
end
orig=x;

T=2048;  % For barbara 2048 % else 1024
mmx=T;
xm=mapping_128;
xm=xm(:);

%==========================================
mat=trans(x,'bior3.7');
mat=fix(mat);
dect(xm)=mat;  % Original
%==========================================
clear x

seqt(128*128)=0;

global iii
iii=0;
LIS=[];
% ==============================================================================================
% Initial lists
for jh=65:256
    LIS=[LIS jh 0];%[2 0 3 0 4 0];      % List of Insignificant Sets ,A-0 & B-1, Co-ordinates
end
LIP=[1:256];%[1 2 3 4];          % List of Insignificant Pixels, Co-ordinates
LSP=[];                 % List of Significant Pixels, Co-ordinates

% Initializaton complete, Starting processing

for xx=1:8

    getlsp=LSP;
%------------------------------------------------------------------------------------------------
for ii=1:length(LIP) % i.e for each entry in the LIP, do
    bit=input_bit;
    if bit=='1'
        LSP=[LSP LIP(ii)];       
        bit=input_bit;
        if bit=='1'
            seqt(LIP(ii))=-mean([T 2*T]);
        else
            seqt(LIP(ii))=mean([T 2*T]);
        end
    end
end
%------------------------------------------------------------------------------------------------

% Now remove the common elements i.e perform LIP-LSP
% D(i,j) means all descendants of (i,j), Function Descendants does this
% O(i,j) means offsprings of (i,j), Function offspring does this
% L(i,j) = D(i,j) - O(i,j), Done by mark_proper

LIP=mark_proper(LIP,LSP);

track=[];
ij=1;
while ij<=length(LIS)
%for ij=1:2:length(LIS) % For each entry in LIS

    if LIS(ij+1)==0 % A type        
      
        out=0;out1=0;% LIS(ij)  % problem
        %-----------------------------------------------------------------------------
        bit=input_bit;
        if bit=='1' % Check for offsprings of ii
            out=1;
        end
        
        if out==1
            % Star 1
            var1=offspring(LIS(ij));
            for kl=1:4
                bit=input_bit;
                if bit=='1'
                    out1=1;
                end
                
                if out1==1
                    LSP=[LSP var1(kl)];
                    
                    % sign
                    bit=input_bit;
                    if bit=='0'
                        seqt(var1(kl))=mean([T 2*T]);
                    else
                        seqt(var1(kl))=-mean([T 2*T]);
                    end
                    
                else
                    LIP=[LIP var1(kl)];
                end
                out1=0;
            end
            
            % Star 2
            lij=mark_proper(descendants(LIS(ij)),offspring(LIS(ij)));
            
            if ~isempty(lij)
                % Move ij to the end of LIS as an entry of type B
                LIS=[LIS LIS(ij) 1];
                track=[track ij];
            else
                track=[track ij];
            end   
        end
    end
    
    if LIS(ij+1)==1 % B type
        out=0;
        lij=mark_proper(descendants(LIS(ij)),offspring(LIS(ij)));
        bit=input_bit;
        if bit=='1'
            out=1;
        end
        
        if out==1
            var1=offspring(LIS(ij));
            for mn=1:4
                LIS=[LIS var1(mn) 0];
            end
            track=[track ij];
        end   
        out=0;
    end
    ij=ij+2;
end

if ~isempty(LIS)
    
    % ----------------------------------------------------
    % Remove repeating elements
    for z=1:length(track)
        LIS(track(z):track(z)+1)=9999;
    end
    % ----------------------------------------------------
    LIS=LIS(find(LIS~=9999));
end

if ~isempty(getlsp)
    seqt=irefinement(seqt,T,getlsp);
end

rec=round(reshape(seqt(xm(:)),128,128));
%xr=idwt2d(rec,fopt,4);
xr=itrans(rec,'bior3.7');
figure,imshow(mat2gray(fix(xr)))
[snr , msr]=PSNR(xr,orig);

% ******************************************************

format short g
[T iii snr msr max(abs(abs(seqt)-abs(dect)))] 
format long
% ******************************************************

T=T/2;
end
%round(reshape(seqt(x(:)),128,128))


%LSP
%LIP
%LIS