%Mats K 010206
iptsetpref('ImshowTruesize','manual')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Load the original image:
%%%%%%%%%%%%%%%%%%%%%%%%
   I=imread('start64.tif');
   N=64; %The size of the image
   L=5;  %The number of Potts states
   I=double(I);  
   figure(1),imshow(I,[1 5])
   title('Original image')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Create the noise-corrupted image:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   n=.75*randn(N); 
%Add the Gaussian noise to the original image and round:
   Idegr=round(I+n);
%Truncate the values outside of {1,...5}:
   Idegr(find(Idegr>5))=5*ones(length(find(Idegr>5)),1);
   Idegr(find(Idegr<1))=1*ones(length(find(Idegr<1)),1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2),imshow(Idegr,[1 5])
title('Noise corrupted image')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Initialize parameters:
%%%%%%%%%%%%%%%%%%%%%%%
Inew=Idegr;
K=10;
C=.5;
figure(3)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Simulation
%After each sweep, there is a break so you can see present result.
%%%%%%%%%%%%%%%
for s=1:K
    %Calculate the temperature for this sweep:
    T=C/log(1+s);
    for k=1:N
        for l=1:N
            %Look up the neighbours:
            n=neighbours8(Inew,[l,k]);
            %Update using a `Metropolis' update procedure
            Inew(k,l)=geman_update(Inew(k,l),n,5,T);
        end
    end
    T,imshow(Inew,[1 5]),
    title(['Image after ' num2str(s) ' iterations of simulated annealing'])
    disp('Press a key to continue'),pause
    disp('Here we go again')
end
figure(3),imshow(Inew,[1 5])
title(['Image after ' num2str(K) 'iteration of simulated annealing'])
