% isingminmac.m tries to minimize the Ising energy pixel by pixel; for the Macintosh.

clf

subplot(3,2,1)
load grapple                         % load G from the file grapple.mat
G = G(1:300,1:300);                  % extract a subimage to use less memory
s = size(G);                         % image dimensions
imagesc(G);                          % display the image with 2 gray levels
title('Original image')
colormap(gray)                       % will display black and white

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
imagesc(Y);                          % *** LOOKS mostly black on Mac, is OK ***
title('Image with p=0.1 multiplicative noise');

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:s(1)-1,                             % go through pixels, avoid edges
  for b=2:s(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)
   imagesc(X);
   title('Image after one sweep');
   drawnow
 else
   subplot(3,2,5)
   imagesc(X);
   title('Image 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');
%imagesc(Z);
%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);