func_spiht_enc.m

4-D SPHIT 压缩Matlab 程序

function out = func_SPIHT_Enc(M,max_bits,spatial_levels,t_levels)
%
% Matlab implementation of 4D - SPIHT (without Arithmatic coding stage)
%
% Encoder
%
% input:    M : input video in wavelet domain
%           max_bits : maximum bits can be used
%           spatial_levels : wavelet decomposition level in z-dimension
%           t_levels:wabelet decomposition level in t-dimension
%
% 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.
% 
%

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

sum(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 4]) = [size(M,1) n_max spatial_levels t_levels]; %wrong???????
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

%添加z-dimension
len2 = size(LIS,1);
zLIP=LIP; 
zLIS=LIS;
for i=1:size(M,3)-1 
    LIP=[LIP;zLIP]; 
    LIP(i*len+1:(i+1)*len,3) = (i+1)*ones(len,1);   
end 
for i=1:size(M,3)/2-1 
    LIS=[LIS;zLIS]; 
    LIS(i*len2+1:(i+1)*len2,3) = (i+1)*ones(len2,1);   
end 
%%%LIS(:, 4) = zeros(size(M,3)/2*len2, 1);

%添加t-dimension
len3 = size(LIS,1);
tLIP = LIP;
tLIS = LIS;
for i = 1:size(M,4)-1
    LIP = [LIP;tLIP];
    LIP(i*len3+1 : (i+1)*len3,4) = (i+1)*ones(len3,1);
end
for i = 1:size(M,4)-1
    LIP = [LIP;tLIP];
    LIP(i*len3+1 : (i+1)*len3,4) = (i+1)*ones(len3,1);
end

LIS(:, 5) = zeros(size(M,4)/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),LIPtemp(i,4))) >= 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).LIPtemp(i,4))>=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,5) == 0
            if bitctr >= max_bits
                return
            end
            max_d = func_Descendant(LIStemp(i,1),LIStemp(i,2),LIStemp(i,3),LIStemp(i,4),LIStemp(i,5),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); t = LIStemp(i,4) ;
                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*t-1)) >= 2^n
                   LSP = [LSP; 2*x-1 2*y-1 2*z-1 2*t-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,2*t)>=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; LISx LISy LISz LISt];
                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*t-1)) >= 2^n
                    LSP = [LSP; 2*x-1 2*y 2*z-1 2*t-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,2*t-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 2*t-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*t-1)) >= 2^n
                    LSP = [LSP; 2*x 2*y-1 2*z-1 2*t-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,2*t-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 2*t-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*t-1)) >= 2^n
                    LSP = [LSP; 2*x 2*y 2*z-1 2*t-1];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;
                    sgn = M(2*x,2*y,2*z-1,2*t-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 2*t-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*t-1)) >= 2^n
                    LSP = [LSP; 2*x-1 2*y-1 2*z 2*t-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,2*t-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 2*t-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,2*t-1)) >= 2^n
                    LSP = [LSP; 2*x-1 2*y 2*z 2*t-1];
                    out(index) = 1; bitctr = bitctr + 1;
                    index = index +1; Bits_LIS = Bits_LIS + 1;