im = +(rand(128, 128) > 0.5)*2 - 1; 
imagesc((im+1)/2);colormap(gray);

% beta negative: low energy configuration has more chance to appear
beta = -1;

for n=1:1e3
    im = gibbs(im,beta);
    imagesc((im+1)/2); colormap(gray)
    title(sprintf('gibbs iter %d',n));
    drawnow;
end



