huff03.m

图像压缩的Huffman编码与解码,matlab源代码

function varargout = huff03(arg1, arg2, arg3)
% Huffman encoder/decoder, one or two vectors of non-negative
% integers are Huffman coded.
% Note that number of arguments decide whether it is encoding og decoding
% If it is more than one input, encoding is done
%
% [y, br, bre] = huff03(x1, x2, Speed);       % encoding
% [y, br] = huff03(x1, x2);                   % encoding
% y = huff03(x1,0);                           % encoding
% x1 = huff03(y);                             % decoding
% [x1, x2] = huff03(y);                       % decoding
% ------------------------------------------------------------------
% Arguments:
%  x1, x2   a column vector of non-negative integers representing the
%           symbol sequence.
%           Best results are achieved if x1, x2 are mostly small integers
%  y        a column vector of non-negative integers (bytes) representing 
%           the code, 0 <= y(i) <= 255. Note that the first part of y will
%           be the Huffman-table/tree information.
%  Speed    Set this to 1 to cheat during encoding.
%           y will then be a sequence of zeros only, but it will be
%           of correct length and the other output arguments will 
%           be correct. If last argumet is scalar, then it is assumed
%           to be 'Speed'.
% Optional output arguments
%  br       Actual bit rate, br(1) for x1, br(2) for x2, br(3) overall
%  bre      0th order entropy (minimum bit rate) of the symbols. 
%           bre(1) is for x1, bre(2) is for x2, bre(3) is overall
% ------------------------------------------------------------------

% SOME NOTES ON THE FUNCTION
% huff03 depens on other functions for Huffman code, and the functions in this file
%  hufflen  - find length of codewords
%  huffcode - find huffman codewords
%  hufftree - find huffman tree

%----------------------------------------------------------------------
%----------------------------------------------------------------------

global y Byte BitPos Speed Level TryAll

Level=0;     % start level
TryAll=1;    % 1 - best compression   0 - fastest compression

% check input and output arguments
if (nargin < 1); 
   error('huff03: function must have input arguments, see help.'); 
end
if (nargout < 1); 
   error('huff03: function must have output arguments, see help.'); 
end
if (nargin == 1)
   Encode=0;Decode=1;
else
   Encode=1;Decode=0;
end;   

% assigne values to arguments
if Decode
   y=arg1(:);
end
if Encode
   x1=arg1(:);
   if ((nargin==2) & (length(arg2(:))==1))
      Speed=arg2;x2=[];
   elseif ((nargin==2) & (length(arg2(:))>1))
      x2=arg2(:);
      Speed=0;
   elseif ((nargin==3) & (length(arg3(:))==1))
      x2=arg2(:);
      Speed=arg3;
   else
      error('huff03: wrong input arguments, see help.'); 
   end
end
clear arg1 arg2 arg3
   
if Encode
   L1=length(x1);L2=length(x2);
   br=[0,0,0];bre=[0,0,0];
   % initalize the global variables
   y=zeros(2*(L1+L2),1);
   Byte=0;BitPos=1;  % ready to write into first position
   % start encoding, first bit indicate one or two sequences
   if (L2 == 0) 
      if Speed
         BitPos=BitPos-1;
         if (~BitPos); Byte=Byte+1; BitPos=8; end; 
      else
         PutBit(0); 
      end;
      [bits1, ent1]=EncodeVector(x1);
      ent2=0;bits2=0;
   else
      if Speed
         BitPos=BitPos-1;
         if (~BitPos); Byte=Byte+1; BitPos=8; end; 
      else
         PutBit(1); 
      end;
      [bits1, ent1]=EncodeVector(x1);
      [bits2, ent2]=EncodeVector(x2);
   end
   y=y(1:Byte);   
   varargout(1) = {y};
   if (nargout >= 2)
      if (L1>0); br(1)=bits1/L1; end;
      if (L2>0); br(2)=bits2/L2; end;
      br(3)=(1+bits1+bits2)/(L1+L2);
      varargout(2) = {br};
   end
   if (nargout >= 3)
      bre(1)=ent1;
      bre(2)=ent2;
      bre(3)=(L1*ent1+L2*ent2)/(L1+L2);
      varargout(3) = {bre};
   end
end
if Decode
   % initalize the global variables, y is set earlier
   Byte=0;BitPos=1;  % ready to read from first position
   % start decoding, is it one or two sequences
   if GetBit         
      x1=DecodeVector;
      x2=DecodeVector;
   else
      x1=DecodeVector;
      x2=[];
   end
   if (nargout == 1)
      if (length(x2)>0)
         x1=[x1;x2];
         warning('huff03: two vectors are concatenated.');
      end
      varargout(1) = {x1};
   end
   if (nargout == 2)
      varargout(1) = {x1};
      varargout(2) = {x2};
   end
end

return     % end of main function, huff03

% the EncodeVector and DecodeVector functions are the ones
% where actual coding is going on.
function [bits, ent] = EncodeVector(x, bits, L, HL, Method, Maxx, Meanx)
global y Byte BitPos Speed Level TryAll
Level = Level + 1;
if (nargin==1)            % this is the same as (Level==1)
   L=length(x);
   Maxx=max(x);
   Meanx=mean(x);
   % find the histogram, hist function is slow if many bins
   if (Maxx < 72)
      Hi=hist(x,0:Maxx);
   else
      Hi=zeros(Maxx+1,1);
      for l=1:L; Hi(x(l)+1)=Hi(x(l)+1)+1; end;
   end
   Hinz=nonzeros(Hi);
   ent=log2(L)-sum(Hinz.*log2(Hinz))/L;  % find entropy
   HL=hufflen(Hi);
   if ((Maxx>250) | TryAll)
      HLlen=HuffTabLen(HL,0);
   else
      HLlen=HuffTabLen(HL,1);
   end   
   I=find(HLlen==min(HLlen));Method=I(1);
   % find number of bits to use
   bits=6;
   bits=bits+HLlen(Method);
   bits=bits+sum(HL.*Hi);
   if ((Maxx<250) & TryAll)
      % test if 'non-optimal' Huffman Code would give fewer bits 
      % (due to shorte table), this makes the function a little bit slower
      HL1=hufflen(Hi+1);
      HLlen1=HuffTabLen(HL1,0);
      I1=find(HLlen1==min(HLlen1));Method1=I1(1);
      % find number of bits to use
      bits1=6;
      bits1=bits1+HLlen1(Method1);
      bits1=bits1+sum(HL1.*Hi);
      if (bits1 < bits)               % Use these codeword lengths instead
         bits=bits1;Method=Method1;
         HL=HL1;HLlen=HLlen1;
      end
   end
   if (L>=16); bits=bits+4; end;
   if (L>=272); bits=bits+4; end;
   if (L>=4368); bits=bits+4; end;
else                % arguments are given, do not need to be calculated
   ent=0;
end
%
% Here we have: x, bits, L, HL, Method, Maxx, Meanx, ent
if ((Level < 8) & TryAll & (L>10))           % maximum level is set to 8
   xm=median(x);
   % xm=mean(x);
   x1=zeros(L,1);x2=zeros(L,1);
   x2(1)=x(1);i1=0;i2=1;
   for i=2:L
      if (x(i-1) <= xm) 
         i1=i1+1; x1(i1)=x(i);
      else
         i2=i2+1; x2(i2)=x(i);
      end
   end
   x1=x1(1:i1);x2=x2(1:i2);
   % find bits1 and bits2 for x1 and x2
   L1=length(x1);L2=length(x2);
   Maxx1=max(x1);Maxx2=max(x2);
   Meanx1=mean(x1);Meanx2=mean(x2);
   if (Maxx1 < 72)
      Hi1=hist(x1,0:Maxx1);
   else
      Hi1=zeros(Maxx1+1,1);
      for l=1:L1; Hi1(x1(l)+1)=Hi1(x1(l)+1)+1; end;
   end
   if (Maxx2 < 72)
      Hi2=hist(x2,0:Maxx2);
   else
      Hi2=zeros(Maxx2+1,1);
      for l=1:L2; Hi2(x2(l)+1)=Hi2(x2(l)+1)+1; end;
   end
   HL1=hufflen(Hi1);HL2=hufflen(Hi2);
   if ((Maxx1>250) | TryAll); HLlen1=HuffTabLen(HL1,0); else; HLlen1=HuffTabLen(HL1,1); end;
   if ((Maxx2>250) | TryAll); HLlen2=HuffTabLen(HL2,0); else; HLlen2=HuffTabLen(HL2,1); end;
   I=find(HLlen1==min(HLlen1));Method1=I(1);
   I=find(HLlen2==min(HLlen2));Method2=I(1);
   bits1=6+HLlen1(Method1)+sum(HL1.*Hi1);
   bits2=6+HLlen2(Method2)+sum(HL2.*Hi2);
   if ((Maxx1<250) & TryAll)
      HL3=hufflen(Hi1+1);
      HLlen3=HuffTabLen(HL3,0);
      I=find(HLlen3==min(HLlen3));Method3=I(1);
      bits3=6+HLlen3(Method3)+sum(HL3.*Hi1);
      if (bits3 < bits1); bits1=bits3;Method1=Method3;HL1=HL3;HLlen1=HLlen3; end;
   end
   if ((Maxx2<250) & TryAll)
      HL3=hufflen(Hi2+1);
      HLlen3=HuffTabLen(HL3,0);
      I=find(HLlen3==min(HLlen3));Method3=I(1);
      bits3=6+HLlen3(Method3)+sum(HL3.*Hi2);
      if (bits3 < bits2); bits2=bits3;Method2=Method3;HL2=HL3;HLlen2=HLlen3; end;
   end
   if (L1>=16); bits1=bits1+4; end;
   if (L1>=272); bits1=bits1+4; end;
   if (L1>=4368); bits1=bits1+4; end;
   if (L2>=16); bits2=bits2+4; end;
   if (L2>=272); bits2=bits2+4; end;
   if (L2>=4368); bits2=bits2+4; end;
else
   bits1=bits;bits2=bits;
end
% Here we may have: x1, bits1, L1, HL1, Method1, Maxx1, Meanx1
% and               x2, bits2, L2, HL2, Method2, Maxx2, Meanx2
% but at least we have bits1 and bits2
if ((bits1+bits2) < bits)
   if Speed
      BitPos=BitPos-1;
      if (~BitPos); Byte=Byte+1; BitPos=8; end; 
   else
      PutBit(1); 
   end;
   [bits1, temp] = EncodeVector(x1, bits1, L1, HL1, Method1, Maxx1, Meanx1);
   [bits1, temp] = EncodeVector(x2, bits2, L2, HL2, Method2, Maxx2, Meanx2);
   bits=bits1+bits2+1;
else
   bits=bits+1; 
   % disp(['huff03-EncodeVector: Level=',int2str(Level),'  ',int2str(L),...
   %       ' sybols stored in ',int2str(bits),' bits.']);
   if Speed
      % advance Byte and BitPos without writing to y
      Byte=Byte+floor(bits/8);
      BitPos=BitPos-mod(bits,8);
      if (BitPos<=0); BitPos=BitPos+8; Byte=Byte+1; end;
   else
      % put the bits into y
      StartPos=Byte*8-BitPos;
      PutBit(0);
      PutVLIC(L);        % max for L is 69903
      PutHuffTab(HL,Method);
      HK=huffcode(HL);
      for i=1:L;
         n=x(i)+1;    % symbol number (value 0 is first symbol, symbol 1)
         for k=1:HL(n)
            PutBit(HK(n,k));
         end
      end
      % check if one has used as many bits as calculated
      BitsUsed=Byte*8-BitPos-StartPos;
      if (BitsUsed~=bits)
         disp(['BitsUsed=',int2str(BitsUsed),'  bits=',int2str(bits)]);
         error('huff03-EncodeVector: Logical error, (BitsUsed~=bits).'); 
      end
   end
end
Level = Level - 1;
return    % end of EncodeVector

function x = DecodeVector
global y Byte BitPos
if GetBit
   x1=DecodeVector;
   x2=DecodeVector;
   L=length(x1)+length(x2);
   xm=median([x1;x2]);
   % xm=mean([x1;x2]);
   x=zeros(L,1);
   x(1)=x2(1);
   i1=0;i2=1;
   for i=2:L
      if (x(i-1) <= xm) 
         i1=i1+1; x(i)=x1(i1);
      else
         i2=i2+1; x(i)=x2(i2);
      end
   end
else
   L=GetVLIC;
   x=zeros(L,1);
   HL=GetHuffTab;
   Htree=hufftree(HL);
   root=1;pos=root;     
   l=0;  % number of symbols decoded so far
   while l<L
      if GetBit
         pos=Htree(pos,3);
      else
         pos=Htree(pos,2);
      end
      if Htree(pos,1)           % we have arrived at a leaf
         l=l+1;
         x(l)=Htree(pos,2)-1;   % value is one less than symbol number
         pos=root;              % start at root again
      end
   end
end
return    % end of DecodeVector

% Functions to write and read the Huffman Table Information
% Several possible schemes may be used to store the 
% Huffman Codeword Lengths, normally the best one is
% chosen. For decoding we do not need to know the method