% Usage: halftone = floyd(original)
%
% Created by Ramesh Neelamani, Rice University on 11/9/99
%
% This gives the Floyd-Stienberg halftoned image of the input image
% \bibitem{Floyd}
% R.~W. Floyd and L.~Stienberg, ``An adaptive algoritm for spatial grayscale,''
% {\em Proc. SID}, vol.~17, no.~2, pp.~75--77, 1976.
%

function[halftoned]  = floyd(x)
 x = double(x);
 
[M,N] = size(x);
T = 127.5;
y = x;
error = 0;

for rows = 1:M-1  

  z(rows,1) =255*(y(rows,1)>=T);
  error = -z(rows,1) + y(rows,1);
  y(rows,1+1) = 7/16 * error + y(rows,1+1);
  y(rows+1,1+1) = 1/16 * error + y(rows+1,1+1);
  y(rows+1,1) = 5/16 * error + y(rows+1,1);

  for cols = 2:N-1
    z(rows,cols) =255*(y(rows,cols)>=T);
    error = -z(rows,cols) + y(rows,cols);
    y(rows,cols+1) = 7/16 * error + y(rows,cols+1);
    y(rows+1,cols+1) = 1/16 * error + y(rows+1,cols+1);
    y(rows+1,cols) = 5/16 * error + y(rows+1,cols);
    y(rows+1,cols-1) = 3/16 * error + y(rows+1,cols-1);
  end

  z(rows,N) =255*(y(rows,N)>=T);
  error = -z(rows,N) + y(rows,N);
  y(rows+1,N) = 5/16 * error + y(rows+1,N);
  y(rows+1,N-1) = 3/16 * error + y(rows+1,N-1);

end
rows = M;
if rows == M
  z(rows,1) =255*(y(rows,1)>=T);
  error = -z(rows,1) + y(rows,1);
  y(rows,1+1) = 7/16 * error + y(rows,1+1);
  
  for cols = 2:N-1
    z(rows,cols) =255*(y(rows,cols)>=T);
    error = -z(rows,cols) + y(rows,cols);
    y(rows,cols+1) = 7/16 * error + y(rows,cols+1);
  end

  z(rows,N) =255*(y(rows,N)>=T);
  error = -z(rows,N) + y(rows,N);
end

halftoned = z;

