%This uses MLE and ICM algorithms to segment (classify each pixel) a noisy X-ray image.
%MLE segmentation is noisy; MAP (ICM algorithm) has much less noise.
clear;load foot.mat;sigma=30;beta=1000;N=size(X,1);
Y=X+sigma*randn(N,N);%Add noise 
figure,imagesc(Y),colormap(gray),axis off,title('Noisy X-ray image')
figure,hist(Y(:),1500);title('Histogram of noisy image')
%subplot(212),plot([0:300],FY(1:301)),title('Zoom of histogram of noisy image')
XHAT(N,N)=0;f=[0 150 200];JMAX=length(f);%3 levels f from histogram (black,flesh,bone)
%Compute MLE segmentation:
for J=1:JMAX;E1(:,:,J)=(Y-f(J)).^2;end;
[Z1,K1]=min(E1,[],3);XHAT=f(K1);%min(E1,[],3) is along dim=3.
figure,imagesc(XHAT),colormap(gray),axis off,title('MLE classification-still noisy')
%ICM algorithm for MAP segmentation
for I=1:9;for J=1:JMAX;
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),axis off,title('MAP (ICM algorithm) segmentation')
