EMSeg.m

function [mask,mu,v,p]=EMSeg(ima,k)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Expectation Maximization image segmentation
%
%   Input:
%          ima: grey color image
%          k: Number of classes
%   Output:
%          mask: clasification image mask
%          mu: vector of class means 
%          v: vector of class variances
%          p: vector of class proportions   
%
%   Example: [mask,mu,v,p]=EMSeg(image,3);
%
%   Author: Prof. Jose Vicente Manjon Herrera
%    Email: jmanjon@fis.upv.es
%     Date: 02-05-2006
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% check image
ima=double(ima);
copy=ima;           % make a copy
ima=ima(:);         % vectorize ima
mi=min(ima);        % deal with negative 
ima=ima-mi+1;       % and zero values
m=max(ima);
s=length(ima);

% create image histogram

h=histogram(ima);
x=find(h);
h=h(x);
x=x(:);h=h(:);

% initiate parameters

mu=(1:k)*m/(k+1);
v=ones(1,k)*m;
p=ones(1,k)*1/k;

% start process

sml = mean(diff(x))/1000;
while(1)
        % Expectation
        prb = distribution(mu,v,p,x);
        scal = sum(prb,2)+eps;
        loglik=sum(h.*log(scal));
        
        %Maximizarion
        for j=1:k
                pp=h.*prb(:,j)./scal;
                p(j) = sum(pp);
                mu(j) = sum(x.*pp)/p(j);
                vr = (x-mu(j));
                v(j)=sum(vr.*vr.*pp)/p(j)+sml;
        end
        p = p + 1e-3;
        p = p/sum(p);

        % Exit condition
        prb = distribution(mu,v,p,x);
        scal = sum(prb,2)+eps;
        nloglik=sum(h.*log(scal));                
        if((nloglik-loglik)<0.0001) break; end;        

        clf
        plot(x,h);
        hold on
        plot(x,prb,'g--')
        plot(x,sum(prb,2),'r')
        drawnow
end

% calculate mask
mu=mu+mi-1;   % recover real range
s=size(copy);
mask=zeros(s);

for i=1:s(1),
for j=1:s(2),
  for n=1:k
    c(n)=distribution(mu(n),v(n),p(n),copy(i,j)); 
  end
  a=find(c==max(c));  
  mask(i,j)=a(1);
end
end


function y=distribution(m,v,g,x)
x=x(:);
m=m(:);
v=v(:);
g=g(:);
for i=1:size(m,1)
   d = x-m(i);
   amp = g(i)/sqrt(2*pi*v(i));
   y(:,i) = amp*exp(-0.5 * (d.*d)/v(i));
end


function[h]=histogram(datos)
datos=datos(:);
ind=find(isnan(datos)==1);
datos(ind)=0;
ind=find(isinf(datos)==1);
datos(ind)=0;
tam=length(datos);
m=ceil(max(datos))+1;
h=zeros(1,m);
for i=1:tam,
    f=floor(datos(i));    
    if(f>0 & f<(m-1))        
        a2=datos(i)-f;
        a1=1-a2;
        h(f)  =h(f)  + a1;      
        h(f+1)=h(f+1)+ a2;                          
    end;
end;
h=conv(h,[1,2,3,2,1]);
h=h(3:(length(h)-2));
h=h/sum(h);