% Code started 30th 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 indices
format short
% Prepare for RASTER scan [Sha93]
xm=mapping_128;
xm=xm(:);
% ####################
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
%mat=dwt2d(x,fopt,4);
mat=trans(x,'bior3.7');
mat=fix(mat);
seqt(xm)=mat;
T=2^fix(log2(max(max(abs(mat)))));
global maxy
maxy=T;
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
output=[];

% Initializaton complete, Starting processing
bit_number=1;

for xx=1:8

    sendlsp=LSP;
    %------------------------------------------------------------------------------------------------
for ii=1:length(LIP) % i.e for each entry in the LIP, do
    if abs(seqt(LIP(ii)))>=T
        output=[output '1'];
        LSP=[LSP LIP(ii)]; % Moving ii to LSP, Removing ii from LSP done later
        % Output sign
        if seqt(LIP(ii))>=0
            output=[output '0'];
        else
            output=[output '1'];
        end
    else
        output=[output '0'];
    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
        %-----------------------------------------------------------------------------
        if ~isempty(find(abs(seqt(descendants(LIS(ij))))>=T)) % Check for offsprings of ii
            output=[output '1'];
            out=1;
            
        else
            output=[output '0'];
        end
        
        if out==1
            % Star 1
            var1=offspring(LIS(ij));
            for kl=1:4
                if abs(seqt(var1(kl)))>=T
                    output=[output '1'];
                    out1=1;
                else
                    output=[output '0'];
                end
                
                if out1==1
                    LSP=[LSP var1(kl)];
                    
                    % sign
                    if seqt(var1(kl))>=0
                        output=[output '0'];
                    else
                        output=[output '1'];
                    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)));
        if ~isempty(find(abs(seqt(lij))>=T))
            output=[output '1'];
            out=1;
        else
            output=[output '0'];
        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
% ----------------------------------------------------
% Remove repeating elements
if ~isempty(LIS)
    
    for z=1:length(track)
        LIS(track(z):track(z)+1)=9999;
    end
    % ----------------------------------------------------
    LIS=LIS(find(LIS~=9999));

end

if ~isempty(sendlsp)
    output=refinement(output,sendlsp,seqt,bit_number);
end
bit_number=bit_number+1;

[  T  length(LSP) length(find(abs(seqt)>=T))]

T=T/2;

end

disp(' OVER ')
save filename output
%output
%LSP
%LIP
%LIS