%This implements the ICM algorithm on a noisy 64X64 Shepp-Logan phantom, 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.
%Note: This example is not in the text, but it is a good example of classification.
clear;load shepp64;sigma=0.05;beta=100;N=size(X,1);
Y=X+sigma*randn(N,N);%Add noise to 64X64 Shepp-Logan phantom.
figure,imagesc(Y),colormap(gray),title('Noisy 64X64 Shepp-Logan phantom')
figure,hist(Y(:),4096),title('Histogram of noisy 64X64 Shepp-Logan phantom')
f=[0,0.2,0.3,1];JMAX=length(f);XHAT(N,N)=0;%4 levels;values from histogram peaks.
%There is a peak at 0.3 but it is drowned out by the noise.
%MLE classifier:
E1(:,:,1)=(Y-f(1)).^2;E1(:,:,2)=(Y-f(2)).^2;
E1(:,:,3)=(Y-f(3)).^2;E1(:,:,4)=(Y-f(4)).^2;
[Z1,K1] = min(E1,[],3);XHAT=f(K1);
figure,imagesc(XHAT),colormap(gray),title('MLE segmentation-still noisy')
%ICM algorithm for classification:
for I=1:9;for J=1:JMAX;%I=iteration number.
E(:,:,J)=([XHAT(2:N,:);XHAT(1,:)]~=f(J))+([XHAT(N,:);XHAT(1:N-1,:)]~=f(J));
E(:,:,J)=E(:,:,J)+([XHAT(:,2:N) XHAT(:,1)]~=f(J))+([XHAT(:,N) XHAT(:,1:N-1)]~=f(J));
E(:,:,J)=2*beta*E(:,:,J)+(Y-f(J)).^2/2/sigma/sigma;end;
[Z,K] = min(E,[],3);XHAT=f(K);end;
figure,imagesc(XHAT),colormap(gray),title('ICM segmentation-almost noiseless')
