function [xdim, ydim, bitstring, optlength] = encode(input_image);
%ENCODE   Arithmetic coding of binary images.
%   [XDIM, YDIM, BS, OL] = ENCODE(I) returns the bitstring BS containing
%   the compressed image, the dimensions XDIM and YDIM of the image I and
%   the shortest theoretical coding length OL.  Predicts the pixels using
%   a simple Bernoulli-model.

[xdim, ydim] = size(input_image);
bufsize = 53;
bs2 = 2^bufsize;
lo = 0;
hi = 1;
memorybits = 0;
bitstring = [];
optlength = 0;

% Initialize the frequency of zeros
freq0 = 1;
total = 2*freq0;

% Loop through all pixels of image
for y = 1:ydim
  for x = 1:xdim
    p0 = freq0/total;
    % The threshold has to be rounded to the same accuracy as the value
    thresh = round((lo+(hi-lo)*p0)*bs2)/bs2;
    % Make sure that the threshold is strictly between lo and hi
    if (thresh == lo)
      thresh = thresh + 1/bs2;
    elseif (thresh == hi)
      thresh = thresh - 1/bs2;
    end
    if (input_image(y,x))
      lo = thresh;
      optlength = optlength - log(1-p0)/log(2);
    else
      hi = thresh;
      optlength = optlength - log(p0)/log(2);
    end

    while (lo >= 0.5 | hi <= 0.5)
      if (lo >= 0.5)
	% Zoom in the upper half [0.5 1] -> [0 1]
	lo = 2*(lo-0.5);
	hi = 2*(hi-0.5);
	bitstring = [bitstring 1 zeros(1,memorybits)];
      else % (hi <= 0.5)
	% Zoom in the lower half [0 0.5] -> [0 1]
	lo = 2*lo;
	hi = 2*hi;
	bitstring = [bitstring 0 ones(1,memorybits)];
      end
      memorybits = 0;
    end
    while (lo >= 0.25 & hi <= 0.75)
      % Zoom in the middle [0.25 0.75] -> [0 1]
      lo = 2*(lo-0.25);
      hi = 2*(hi-0.25);
      % We don't know yet which bits to write, but keep the number in memory
      memorybits = memorybits + 1;
    end

    % Update the frequencies
    if (input_image(y, x) == 0)
      freq0 = freq0 + 1;
    end
    total = total + 1;
  end
end

% Make sure the value we write is between the limits
if (lo > 0 | hi < 1 | memorybits > 0)
  if (lo == 0)        % (lo == 0 & (hi < 1 | memorybits > 0))
    bitstring = [bitstring 0 ones(1,memorybits)];
  elseif (hi == 1)    % (hi == 1 & lo > 0)
    bitstring = [bitstring 1 zeros(1,memorybits)];
  elseif (lo <= 0.25) % (lo <= 0.25 & hi > 0.5)
    bitstring = [bitstring 0 ones(1,memorybits+1)];
  else                % (lo < 0.5 & hi > 0.75)
    bitstring = [bitstring 1 zeros(1,memorybits+1)];
  end
end
