% file mcmc7.m: Ising model with uniform prior
%
clear all; close all
% input matrix consisting of letter A. A is made of 1's 
% while the background is made of -1's.
F = imread('C:/lettera.bmp'); %or other path...
[M,N] = size(F);
sigma = 2;
d = double(F); d=2*((d-29)./198)-1; %d either -1 1
d = d + sigma*randn(size(d)); %noisy letter A.
%-----------------------------------------------
figure(1);
subplot(1,2,1);imagesc(F);set(gca,'Visible','off');colormap gray; axis square;
subplot(1,2,2);imagesc(d);set(gca,'Visible','off');colormap gray; axis square;
drawnow
J = 0.3; %clumping parameter
f = ones(M,N);
figure(2);
subplot(1,2,1); hf=imagesc(f); set(gca,'Visible','off');
colormap gray; axis square; drawnow;
mf = zeros(M,N);
subplot(1,2,2); hm=imagesc(mf); set(gca,'Visible','off');
colormap gray; axis square; drawnow;
SS = 10000;
misfit = [];
adj = [-1 1 0 0; 0 0 -1 1];
iter = 0;
while 1
   ix = ceil(N*rand(1)); iy = ceil(M*rand(1));
   pos = iy + M*(ix-1); fp = -f(pos);
   LkdRat = exp(d(pos)*(fp - f(pos))/sigma.^2);
   alpha = LkdRat/10;
   if rand<alpha % Prob of acceptance = min(1,alpha)
        f(pos) = fp;
   end
   iter = iter + 1;
   if rem(iter,SS) == 0,
        mf = mf+f; NS = iter/SS; iter
        set(hf,'CData',f);
        set(hm,'CData',mf); drawnow
        if iter/SS > length(misfit)
               misfit = [misfit,zeros(100,1)];
               misfit(iter/SS) = sum(sum((d-f).^2))/sigma;
        end
    end
end