function out = func_SPIHT_Enc(M, max_bits, level)
%
% Matlab implementation of 3D - SPIHT (without Arithmatic coding stage)
%
% Encoder
%
% input:    M : input video in wavelet domain
%           max_bits : maximum bits can be used
%           level : wavelet decomposition level
%
% output:   out : bit stream
%
% This code is based on the implementation of the 2D-SPIHT by Jing Tian. It
% extends the code in order to handle video sequences (3D-SPIHT). Small changes have
% been also made to the 2D-SPIHT code.
% 
%
% Athanasopoulos Dionysios 
% Postgraduate Student
% Computer Engineering and Informatics Dept.
% University of Patras, Greece
%

%-----------   Initialization  -----------------
bitctr = 0;
out = 2*ones(1,max_bits - 14);
n_max = floor(log2(abs(max(max(max(M))))));

sum(sum(sum(floor(log2(abs(M(M~=0)))))))
Bits_Header = 0;
Bits_LSP = 0;
Bits_LIP = 0;
Bits_LIS = 0;

%-----------   output bit stream header   ----------------
% image size, number of bit plane, wavelet decomposition level should be
% written as bit stream header.
out(1,[1 2 3]) = [size(M,1) n_max level]; bitctr = bitctr + 24;
index = 4;
Bits_Header = Bits_Header + 24;

%-----------   Initialize LIP, LSP, LIS   ----------------

temp = [];
% bandsize1 = round (2.^(log2(size(M, 1)) - level + 1));
% bandsize2 = round (2.^(log2(size(M, 2)) - level + 1));
bandsize1 = 2;
bandsize2 = 2;
temp1 = 1 : bandsize1;
for i = 1 : bandsize2
    temp = [temp; temp1];
end
LIP(:, 1) = temp(:);
temp = [];
temp1 = 1 : bandsize2;
for i = 1 : bandsize1
    temp = [temp; temp1];
end
temp = temp';
LIP(:, 2) = temp(:);
len = size(LIP,1);
LIP(:,3) = ones(len,1);

LIS = LIP;
pstart = 1;
pend = bandsize2 / 2;
for i = 1 : bandsize1 / 2
    LIS(pstart : pend, :) = [];
    pdel = pend - pstart + 1;
    pstart = pstart + bandsize2 - pdel;
    pend = pend + bandsize2 - pdel;
end


len2 = size(LIS,1);
tLIP=LIP; 
tLIS=LIS;
for i=1:size(M,3)-1 
    LIP=[LIP;tLIP]; 
    LIP(i*len+1:(i+1)*len,3) = (i+1)*ones(len,1);   
end 
for i=1:size(M,3)/2-1 
    LIS=[LIS;tLIS]; 
    LIS(i*len2+1:(i+1)*len2,3) = (i+1)*ones(len2,1);   
end 
LIS(:, 4) = zeros(size(M,3)/2*len2, 1);

LSP = [];

n = n_max;

%-----------   CODING   ----------------
while(bitctr < max_bits)
        
    % Sorting Pass
     LIPtemp = LIP;temp = 0;
    for i = 1:size(LIPtemp,1)
        temp = temp+1;
        if (bitctr + 1) >= max_bits
            if (bitctr < max_bits)
                out(length(out))=[];
            end
            return
        end
    
        if abs(M(LIPtemp(i,1),LIPtemp(i,2),LIPtemp(i,3))) >= 2^n 
            out(index) = 1; bitctr = bitctr + 1;
            index = index +1; Bits_LIP = Bits_LIP + 1;
            sgn = M(LIPtemp(i,1),LIPtemp(i,2),LIPtemp(i,3))>=0;
            out(index) = sgn; bitctr = bitctr + 1;
            index = index +1; Bits_LIP = Bits_LIP + 1;
            LSP = [LSP; LIPtemp(i,:)];
            LIP(temp,:) = []; temp = temp - 1;
        else
            out(index) = 0; bitctr = bitctr + 1;
            index = index +1;
            Bits_LIP = Bits_LIP + 1;
        end
    end
    
    LIStemp = LIS; temp = 0; i = 1;
    while ( i <= size(LIStemp,1))
        temp = temp + 1;
        if LIStemp(i,4) == 0
            if bitctr >= max_bits
                return
            end
            max_d = func_Descendant(LIStemp(i,1),LIStemp(i,2),LIStemp(i,3),LIStemp(i,4),M);
            if max_d >= 2^n
                out(index) = 1; bitctr = bitctr + 1;
                index = index +1; Bits_LIS = Bits_LIS + 1;
                x = LIStemp(i,1); y = LIStemp(i,2);
                z = LIStemp(i,3); 
                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(M(2*x-1,2*y-1,2*z-1)) >= 2^n
                    LSP = [LSP; 2*x-1 2*y-1 2*z-1];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = M(2*x-1,2*y-1,2*z-1)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x-1 2*y-1 2*z-1];
                end
                
                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(M(2*x-1,2*y,2*z-1)) >= 2^n
                    LSP = [LSP; 2*x-1 2*y 2*z-1];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = M(2*x-1,2*y,2*z-1)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x-1 2*y 2*z-1];
                end
                
                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(M(2*x,2*y-1,2*z-1)) >= 2^n
                    LSP = [LSP; 2*x 2*y-1 2*z-1];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = M(2*x,2*y-1,2*z-1)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x 2*y-1 2*z-1];
                end
                
                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(M(2*x,2*y,2*z-1)) >= 2^n
                    LSP = [LSP; 2*x 2*y 2*z-1];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = M(2*x,2*y,2*z-1)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x 2*y 2*z-1];
                end

                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(M(2*x-1,2*y-1,2*z)) >= 2^n
                    LSP = [LSP; 2*x-1 2*y-1 2*z];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = M(2*x-1,2*y-1,2*z)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x-1 2*y-1 2*z];
                end
                
                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(M(2*x-1,2*y,2*z)) >= 2^n
                    LSP = [LSP; 2*x-1 2*y 2*z];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = M(2*x-1,2*y,2*z)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x-1 2*y 2*z];
                end
                
                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(M(2*x,2*y-1,2*z)) >= 2^n
                    LSP = [LSP; 2*x 2*y-1 2*z];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = M(2*x,2*y-1,2*z)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x 2*y-1 2*z];
                end
                
                if (bitctr + 1) >= max_bits
                    if (bitctr < max_bits)
                        out(length(out))=[];
                    end
                    return
                end
                if abs(M(2*x,2*y,2*z)) >= 2^n
                    LSP = [LSP; 2*x 2*y 2*z];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = M(2*x,2*y,2*z)>=0;
                    out(index) = sgn; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                else
                    out(index) = 0; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    LIP = [LIP; 2*x 2*y 2*z];
                end
  
                if ((2*(2*x)-1) < size(M,1) & (2*(2*y)-1) < size(M,2) & (2*(2*z)-1) < size(M,3))
                    LIS = [LIS; LIStemp(i,1) LIStemp(i,2) LIStemp(i,3) 1];
                    LIStemp = [LIStemp; LIStemp(i,1) LIStemp(i,2) LIStemp(i,3) 1];
                end
                LIS(temp,:) = []; temp = temp-1;
                
            else
                out(index) = 0; bitctr = bitctr + 1;
                index = index +1; Bits_LIS = Bits_LIS + 1;
            end
        else
            if bitctr >= max_bits
                return
            end
            max_d = func_Descendant(LIStemp(i,1),LIStemp(i,2),LIStemp(i,3),LIStemp(i,4),M);
            if max_d >= 2^n
                out(index) = 1; bitctr = bitctr + 1;
                index = index +1;
                x = LIStemp(i,1); y = LIStemp(i,2); z = LIStemp(i,3);
                LIS = [LIS; 2*x-1 2*y-1 2*z-1 0; 2*x-1 2*y 2*z-1 0; 2*x 2*y-1 2*z-1 0; 2*x 2*y 2*z-1 0;2*x-1 2*y-1 2*z 0; 2*x-1 2*y 2*z 0; 2*x 2*y-1 2*z 0; 2*x 2*y 2*z 0];
                LIStemp = [LIStemp; 2*x-1 2*y-1 2*z-1 0; 2*x-1 2*y 2*z-1 0; 2*x 2*y-1 2*z-1 0; 2*x 2*y 2*z-1 0;2*x-1 2*y-1 2*z 0; 2*x-1 2*y 2*z 0; 2*x 2*y-1 2*z 0; 2*x 2*y 2*z 0];
                LIS(temp,:) = []; temp = temp - 1;
            else
                out(index) = 0; bitctr = bitctr + 1;
                index = index +1; Bits_LIS = Bits_LIS + 1;
            end
        end
        i = i+1;
    end
    
    % Refinement Pass
    temp = 1;
    value = floor(abs(2^(n_max-n+1)*M(LSP(temp,1),LSP(temp,2),LSP(temp,3))));
    while (value >= 2^(n_max+2) & (temp <= size(LSP,1)))
        if bitctr >= max_bits
            return
        end
        s = bitget(value,n_max+2);
        out(index) = s; bitctr = bitctr + 1;
        index = index +1; Bits_LSP = Bits_LSP + 1;
        temp = temp + 1;
        if temp <= size(LSP,1)
            value = floor(abs(2^(n_max-n+1)*M(LSP(temp,1),LSP(temp,2),LSP(temp,3))));
        end
    end
    
    n = n - 1;
end