function w=simbin(beta,N,iter,w0)
% Simulates iter steps to obtain a binary image w,
% Let u=the ratio of black neighbours i to j.
%  P( w(j)=0 | w(i), N(j-i+(a+1,b+1))=1 ) =
%  exp( beta * (2*u-1) ) / (1 + exp( ... ))
%
% w=simbin(beta,N,iter,w0);
%
%  beta 1|1
%  N    (2a+1)|(2b+1), non-scaled neighbour-pattern. (0/1)
%  iter 1|1, number of iterations
%  w0   m|n, the initial image 
%

w=w0;
[m,n]=size(w);
[a,b]=size(N);
a=floor(a/2);
b=floor(b/2);

NN=conv2(ones(size(w)),N,'same');
w=[zeros(m,b) w zeros(m,b)];
w=[zeros(a,n+2*b);w;zeros(a,n+2*b)];
for k=1:iter
  for i=1:m
    e=rand(1,n);
    e=(log(e./(1-e))/beta+1)/2;
    for j=1:n
      q=sum(sum(N.*w(i:(i+2*a),j:(j+2*b))))/NN(i,j);
      w(i+a,j+b)=(q>=e(j));
    end
  end
end
w=w(a+(1:m),b+(1:n));
