% This is a test program that tests the EZW encoder and decoder  
% matrix X is an example matrix same as Shapiro's paper given  
  
  
% Written by xiangang li  
%  
% This function calls:  
%    dominantpass.m  
%    subordinatepass.m  
%    checkdescents1.m  
%    checkchildren.m  
%    mapping.m  
  
  
  
X=[63 -34 49 10 7 13 -12 7;  
   -31 23 14 -13 3 4 6 -1;  
   15 14 3 -12 5 -7 3 9;  
   -9 -7 -14 8 4 -2 3 2;  
   -5 9 -1 47 4 6 -2 2;  
   3 0 -3 2 3 -2 0 4;  
   2 -3 6 -4 3 6 3 6;  
   5 11 5 6 0 3 -4 4];  
  
X0=X;  
  
Y0=max(X);  
     Y1=max(Y0); % y1是X中的最大值。  
     for i=0:20;  
        if 2^i<=Y1 & 2^i>0.5*Y1;  
           threshold=2^i;   % get initial threshold T0;  
           initialthreshold=threshold; % get initial threshold T0;  
           laststeplevel=i+1;% last step level  
           break;  
        end;  
     end;  
       
       
     sublist=[];  
     sub_list=[];  
     [xx,yy]=size(X);  
       
    
A=mapping(xx);  
[m,n]=size(A);  
global N;      % Let Morton scanorder vector as a global variable  
 N=zeros(m*n,2);  
  
     for i=1:m,  
        for j=1:n,  
              
                 N(A(i,j),1)=i;  
                 N(A(i,j),2)=j;  
                  
         end  
     end  
  
       order=1;    
while threshold ~= 0.5,  % if threshold~=1, do dominantpass and subordinatepass.  
           
   threshold  
     
   %Dominant Pass  
     
   [D,X,sublist,sub_list] = dominantpass(X,threshold,sublist,sub_list);   
     
   DD{order}=D  
   significantlist{order}=sub_list;  
       
    %Subordinate pass  
      
    threshold=threshold/2;  
    if threshold ==0.5,  
        break;  
    end  
      
      
   S = subordinatepass(sublist,threshold);  
   SS{order}=S  
   order=order+1;  
      
end  
  
    
%*********************************%     
%*********EZW decoder*************%  
%*********************************%  
  
 global N;  
 [m,n]=size(N);% the size of initial image  
               % m is the pixels of initial image  
 XX=zeros(sqrt(m)); % initialize the reconstructured image to zero;  
 threshold=initialthreshold; % initial theshold ;  
   
 sublist=[]; % sublist is the new position matrix   
             % for all significant coefficients 'p' and 'n';  
               
               
for level=1:laststeplevel,   
                  
            RR=zeros(size(XX)); % reference matrix RR;  
    
     [a,b]=size(DD{level}); % ?      
   
 % dominant pass  
   
       
        i=1; j=1;  
     while j<=b,  
       
         if RR(N(i,1),N(i,2))==0  
                        
              if DD{level}(j)=='p'  
                  if threshold==1  
                     XX(N(i,1),N(i,2))=threshold;   
                 else  
                     XX(N(i,1),N(i,2))=1.5*threshold;  
                 end  
                  
              end  
       
              if DD{level}(j)=='n'  
                 if threshold==1  
                     XX(N(i,1),N(i,2))=-threshold;   
                 else  
                     XX(N(i,1),N(i,2))=-1.5*threshold;  
                 end  
              end  
           
                  
              if DD{level}(j)=='t'& A(N(i,1),N(i,2))<=m/4  
                 RR=checkchildren(i,RR);% all zerotree's descendants are set to 1.  
              end  
               
              RR(N(i,1),N(i,2))=1; %reference matrix =1;  
              i=i+1;  
              j=j+1;  
                 
        else i=i+1;  
                 
               
        end  
           
    end  
   
   
 % subordinate pass  
   
   
 [xx,yy]=size(significantlist{level});  
 threshold=threshold/2;  
   
   for i=1:xx,  
         
       if level==laststeplevel|threshold==0.5  
           break  
       end  
         
      if SS{level}(i)==1  
         if XX(sub_list(i,1),sub_list(i,2))>0;  
             
            XX(sub_list(i,1),sub_list(i,2))= fix(XX(sub_list(i,1),sub_list(i,2))+ threshold/2);  
              
         else   
               
            XX(sub_list(i,1),sub_list(i,2))= fix(XX(sub_list(i,1),sub_list(i,2))-threshold/2);  
             
         end  
      end  
        
      if SS{level}(i)==0  
         if XX(sub_list(i,1),sub_list(i,2))>0;  
             
            XX(sub_list(i,1),sub_list(i,2))= fix(XX(sub_list(i,1),sub_list(i,2))-threshold/2);  
              
         else   
               
            XX(sub_list(i,1),sub_list(i,2))= fix(XX(sub_list(i,1),sub_list(i,2))+threshold/2);  
           
         end  
      end  
   end  
     
      
  
      threshold  
      level  
      XX  
        
end  
initialiamge=X0  
reconstructedimage=XX  
difference=XX-X0