% isingminimum.m tries to minimize the Ising energy pixel by pixel

clf

subplot(3,2,1)
G = imread('grapple.gif');           % read the .gif file
s = size(G);                         % image dimensions
imshow(G,[0 1]);                     % display the image with 2 gray levels
title('Original image')

subplot(3,2,2)
p = 0.1;                             % probability of flipping
L = log((1-p)/p);                    % constant in D(x,y)
Y = (G == uint8(rand(s)>p));         % change pixels with probability p
imshow(Y,[0 1]);
title('Image with p=0.1 multiplicative noise');

subplot(3,2,3)
M = min(s,511);                      % upper left corner of pixels
imshow(Y(1:M(1),1:M(2)),[0 1]);
title('Upper left corner');

drawnow
X = Y;                               % original estimate we will refine
beta = 2*0.5;                        % scale beta for our way of computing k,d

k = sum(sum(X(1:s(1)-1,:)~=X(2:s(1),:)));
k = k + sum(sum(X(:,1:s(2)-1)~=X(:,2:s(2)))); % count # differing pairs
d = sum(sum(X~=Y));                  % differences between X and Y
ho= beta*k+L*d;

fprintf('Energy of original estimate (data) is %8.1f\n', ho);

for sw=1:2,                                  % total of two sweeps
 for a=2:M(1)-1,                             % go through pixels, avoid edges
  for b=2:M(2)-1,
    n = ones(4,1)*[a b] + [[1 0]; [0 1]; [-1 0]; [0 -1]];    % neighbors
    h = L*double(X(a,b) ~= Y(a,b));                          % current energy
    hc= L*double(X(a,b) == Y(a,b));                          % changed energy?
    for j=1:4,
            h  = h  + beta*double(X(n(j,1),n(j,2))~=X(a,b));
            hc = hc + beta*double(X(n(j,1),n(j,2))==X(a,b));
    end 
    if hc < h,                               % if energy would go down,
      X(a,b) = (X(a,b)==0);                  % flip this pixel
    end
  end
 end
 if sw==1,
   subplot(3,2,4)
   imshow(X(1:M(1),1:M(2)),[0 1]);
   title('Corner after one sweep');
   drawnow
 else
   subplot(3,2,5)
   imshow(X(1:M(1),1:M(2)),[0 1]);
   title('Corner after two sweeps');
 end
end

k = sum(sum(X(1:s(1)-1,:)~=X(2:s(1),:)));
k = k + sum(sum(X(:,1:s(2)-1)~=X(:,2:s(2))));  % count # differing pairs
d = sum(sum(X~=Y));                            % differences between X and Y
hf= beta*k+L*d;

fprintf('Final energy of estimate is           %8.1f\n', hf);
fprintf('At that rate, final energy would be   %8.1f\n', ho-(prod(s)/prod(M))*(ho-hf));

subplot(3,2,6)
Z = medfilt2(X,[5 5],'symmetric');
imshow(Z,[0 1]);
title('5 by 5 median filtered');

k = sum(sum(Z(1:s(1)-1,:)~=Z(2:s(1),:)));
k = k + sum(sum(Z(:,1:s(2)-1)~=Z(:,2:s(2))));  % count # differing pairs
d = sum(sum(Z~=Y));                            % differences between Z and Y
hm= beta*k+L*d;

fprintf('Energy of median filtered estimate    %8.1f\n', hm);