function y = mat2huff(x)
%MAT2HUFF Huffman encodes a matrix
% Y = MAT2HUFF(X) Huffman encodes matrix X using symbol probabilities,
% unit-width histogram bins between X's minimum and maximum values.   
% The encoded data is returned as a structure Y:
%   Y.code         The Huffman encoded values of X, stored in 
%                  uint16 vector.
%   The other fields of Y contain additional decoding information,
%   including:
%   Y.min          The minimum value of X plus 32768
%   Y.size         The size of X
%   Y.hist         The histogram of X
%
% If X is logical, uint8, uint16, uint32, int8, int16, or double with
% integer values, it can be input directly to MAT2HUFF. The minimum value
% of X must be representable as an int16
%
% If X is double with non-integer values, ---for example, an image with
% values between 0 and 1 --- first scale X to an appropriate integer range
% before the call. For example, use Y = mat2huff(round(255*X)) for 256
% gray level encoding
%
% The number of Huffman code words is round(max(X(:)))-round(min(X(:)))+1
% You may need to scale input X to generate codes of reasonable length. The maximum
% row or column dimension of X is 65535
%
%   See also HUFF2MAT, HUFFMAN

if (ndims(x)~=2) | ~isreal(x) | (~isnumeric(x) & ~islogical(x))
    error('X must be 2-D real numeric or logical matrix.');
end

% Store the size of input x
y.size = uint32(size(x));

% Find the range of x values and store its minimum value biased by +32768 
% as UINT16
x = round(double(x));
xmin = min(x(:));
xmax = max(x(:));
pmin = double(int16(xmin));
pmin = uint16(pmin+32768);
y.min = pmin;

% Compute the input histogram between xmin and xmax with unit width bins,
% scale to UINT16, and store 
x = x(:)';
h = histc(x, xmin:xmax);

if max(h)>65535
    h = 65535*h/max(h);
end
h = uint16(h);
y.hist = h;

% Code the input matrix and store the result.
map = huffman(double(h));           % Make Huffman code map
hx = map(x(:) - xmin + 1);          % Map image
hx = char(hx)';                     % Convert to char array
hx = hx(:)';
hx(hx ==' ') = [];                  % Remove blanks
ysize = ceil(length(hx)/16);        % Compute encoded size
hx16 = repmat('0', 1, ysize*16);    % Pre-allocate modulo-16 vector
hx16(1:length(hx)) = hx;            % make hx modulo 16 in length
hx16 = reshape(hx16,16,ysize);      % Reshape to 16 characters words
hx16 = hx16' - '0';                 % Convert binary string to decimal
twos = pow2(15:-1:0);
y.code = uint16(sum(hx16.*twos(ones(ysize,1),:),2))';