%Realization of Anisotropic Binary Markov Random Field 
 clear all
 close all
%-----------------figure defaults
disp('Unisotripic Ising by Gibbs')
lw = 2;  
set(0, 'DefaultAxesFontSize', 16);
fs = 14;
msize = 5;
randn('state',3) %set the seeds (state) to have 
rand ('state',3) %the constancy of results
%Chess board
clear all
close all
%------------------------------------------
 ima = rand(128) > 0.5;
 frame = ones(128);
% 
 frame( 1,1:end) = 0;
 frame(end,1:end) = 0;
 frame(1:end, 1) = 0;
 frame(1:end,end) = 0;
 white = find(ima.*frame);
 black = find((1-ima).*frame);
 Nblack = length(black);
 Nwhite = length(white);


%Vectorised update of Ising Model
b1=0.1; b2=0.1;
% neigh = ...
%         [0 b2 0; ...
%         b1 0 b1; ...
%         0 b2 0];
neigh = ...
        [b1   0     b2; ...
         0 (b1+b2)/2 0; ...
         b2   0     b1];
alpha = - b1 - b2;
ell=0;
while 1,
for k=1:5
for pass=0:1
  field = filter2(neigh, ima, 'same') + alpha;
       if pass==0
           ima(black) = ...
            1 ./ (1 + exp(-field(black))) > rand(Nblack,1);
       else
           ima(white) = ...
            1 ./ (1 + exp(-field(white))) > rand(Nwhite,1);
       end
end
ell=ell+1;
end
figure(1); imagesc(ima); colormap(gray); title(['Iteration #  ' num2str(ell)]);
drawnow
end