%This implements the ICM algorithm to segment a noisy 4-square image, and
%compare the result with the (simpler) MLE segmentation (to closest level).
%The algorithm works well except near boundaries, where the MRF model breaks down,
%and where several especially-noisy values cluster together.
clear;sigma=0.5;beta=300;f=[0 1];
X(1:32,1:32)=0;X(33:64,33:64)=0;X(1:32,33:64)=1;X(33:64,1:32)=1;%4-square image.
N=size(X,1);Y=X+sigma*randn(N,N);XHAT(N,N)=0;%Add noise to 4-square image.
figure,imagesc(Y),colormap(gray),title('Noisy 4-square image g[n,m]')
figure,hist(Y(:),4000),title('Histogram p(g[n,m]) of noisy 4-square image') 
XHAT(Y>(f(1)+f(2))/2)=f(2);XHAT(Y<(f(1)+f(2))/2)=f(1);%MLE classifier.
figure,imagesc(XHAT),colormap(gray),title('MLE segmentation-still noisy')
%ICM algorithm:
for I=1:9;%ICM iteration number
E1=([XHAT(2:N,:);XHAT(1,:)]~=f(2))+([XHAT(N,:);XHAT(1:N-1,:)]~=f(2));
E1=E1+([XHAT(:,2:N) XHAT(:,1)]~=f(2))+([XHAT(:,N) XHAT(:,1:N-1)]~=f(2));
E1=2*beta*E1+(Y-f(2)).^2/2/sigma/sigma;
E2=([XHAT(2:N,:);XHAT(1,:)]~=f(1))+([XHAT(N,:);XHAT(1:N-1,:)]~=f(1));
E2=E2+([XHAT(:,2:N) XHAT(:,1)]~=f(1))+([XHAT(:,N) XHAT(:,1:N-1)]~=f(1));
E2=2*beta*E2+(Y-f(1)).^2/2/sigma/sigma;
XHAT(E2>E1)=f(2);XHAT(E1>E2)=f(1);end
figure,imagesc(XHAT),colormap(gray),...
title('ICM classifier after 9 iterations-almost noiselss')

