function varargout = estimateBits(xC,ecm,v)
% estimateBits    Estimate bits for coding a cell-array of integer sequences
%
% use:
%  bits = estimateBits(xC, verbose)
%  bits = estimateBits(xC, ecm, verbose)
%  [bits, r] = estimateBits(...)
%  [bits, r, y] = estimateBits(xC, ecm, verbose)
%-------------------------------------------------------------------------
% arguments:
%   xC          a cell array of integer sequences
%   ecm         the entropy coding to use, may be 'none','Huff06','Arith06'
%               or 'Arith07'. If omitted 'none' is used, this estimate bits
%               based on zero-order entropy of each sequence
%               (Huff06, Arith06 and Arith07, are in catalog ...\Comp\ )
%   verbose     0 for no display, 1 for some display, 2 for very verbose
%               (entropy coding may also display some lines)
%   bits        is the (estimated) total number of bits to use
%   r           a matrix that sum up the results, size is (numel(xC)+1) x 4
%               one line for each of the input sequences, the columns are
%               r(:,1) - number of elements in the sequence
%               r(:,2) - zero-order entropy of the sequence
%               r(:,3) - bits needed to code the sequence
%               r(:,4) - bit rate for the sequence, Res(:,3)/Res(:,1)
%               Then the last line is total, bits = r(end,3)
%   y           The byte sequence
%-------------------------------------------------------------------------

%----------------------------------------------------------------------
% Copyright (c) 2010.  Karl Skretting.  All rights reserved.
% University of Stavanger (Stavanger University), Signal Processing Group
% Mail:  karl.skretting@uis.no   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  22.01.2010  KS: was an internal function in 'compressImage.m'
%----------------------------------------------------------------------

if (nargin == 2)
    v = ecm;
    ecm = 'none';
end
if ( ~strcmpi(ecm,'none') &&  (exist(ecm,'file') ~= 2) )
    ecm = 'none';
end
    
if (v > 0)
    disp(['estimateBits: input is ',int2str(numel(xC)),' sequences',...
        ', entropy coding method (ecm) = ',ecm,'.']);
end

if strcmpi(ecm,'none')
    bits = 0;
    r = zeros(numel(xC)+1,4);
    y = [];
    for i = 1:numel(xC)
        if (numel(xC{i}) > 0)
            counts = hist(xC{i}, min(xC{i}):max(xC{i}));
            a = numel(find(counts));
            e = entropy0(counts);
            b = ceil(numel(xC{i})*e) + ceil(a*log2(a));
            if (numel(xC{i}) < 5)
                b = numel(xC{i})*(3+ceil(log2(max(xC{i}))));
            end
            bits = bits+b;
            if (v > 1); 
                disp(['Seq. ',int2str(i),' has length ',int2str(numel(xC{i})),...
                    '. Entropy is ',num2str(e),' and bit estimate is ',...
                    int2str(ceil(numel(xC{i})*e)),' + ',...
                    int2str(ceil(a*log2(a))),' = ',int2str(b),' bits.']);
            end
            r(i,:) = [numel(xC{i}), e, b, b/numel(xC{i})]; 
        end
    end
    r(end,1) = sum(r(1:numel(xC),1));
    r(end,2) = sum(r(1:numel(xC),1).*r(1:numel(xC),2))/r(end,1);
    r(end,3) = bits;
    r(end,4) = bits/r(end,1);
elseif strcmpi(ecm,'Huff06')
    [y,r] = Huff06(xC, 5, 0);
    bits = numel(y)*8;
    if (v > 1); dispLines(r); end;
elseif strcmpi(ecm,'Arith06')
    [y,r] = Arith06(xC);
    bits = numel(y)*8;
    if (v > 1); dispLines(r); end;
elseif strcmpi(ecm,'Arith07')
    [y,r] = Arith07(xC);
    bits = numel(y)*8;
    if (v > 1); dispLines(r); end;
end

if (v > 0)
    disp(['Total is  ',int2str(bits),' bits.']);
end

varargout{1} = bits;
if (nargout > 1);     varargout{2} = r;   end;
if (nargout > 2);     varargout{3} = y;   end;

return


% zero-order entropy
function H = entropy0(S)
N = sum(sum(S));		% if S is probability this is 1
if ((N>0) && (min(S(:))>=0))
   Snz=nonzeros(S);
   H=log2(N)-sum(Snz.*log2(Snz))/N;
else
   H=0;
end
return

function dispLines(r)
for i=1:(size(r,1)-1);
    disp(['Seq. ',int2str(i),' has length ',int2str(r(i,1)),...
        '. Entropy is ',num2str(r(i,2)),' and bit estimate is ',...
        int2str(r(i,3)),' (bit rate ',num2str(r(i,4)),').']);
end
i = size(r,1);
disp(['Total length ',int2str(r(i,1)),...
    '. Average entropy is ',num2str(r(i,2)),' and bit estimate is ',...
    int2str(r(i,3)),' (bit rate ',num2str(r(i,4)),').']);
return