function output_image = decode(xdim, ydim, bitstring);
%DECODE   Arithmetic decoding of binary images.
%   I = DECODE(XDIM, YDIM, BS) returns the XDIM-by-YDIM image decoded from
%   the bitstring BS.  Predicts the pixels using a simple Bernoulli-model.

output_image = zeros(xdim, ydim);
bufsize = 53;
bs2 = 2^bufsize;
lo = 0;
hi = 1;
bitpos = 1;
bitstring=[bitstring zeros(1,bufsize)];
value = 0;

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

% Read in the first 'bufsize' bits
for i=1:bufsize 
  value = 2*value + bitstring(bitpos)/bs2; 
  bitpos = bitpos + 1;
end

% 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 (value >= thresh)
      lo = thresh;
      output_image(y, x) = 1; 
    else
      hi = thresh;
      output_image(y, x) = 0;
    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);
	value = 2*(value-0.5) + bitstring(bitpos)/bs2;
      else % (hi <= 0.5)
	% Zoom in the lower half [0 0.5] -> [0 1]
	lo = 2*lo;
	hi = 2*hi;
	value = 2*value + bitstring(bitpos)/bs2;
      end
      bitpos = bitpos + 1;
    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);
      value = 2*(value-0.25) + bitstring(bitpos)/bs2;
      bitpos = bitpos + 1;
    end

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