# 混合高斯模型用于图像分割

http://download.csdn.net/detail/dayenglish/8384547

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=h/sum(h);

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

prb = distribution(mu,v,p,x);%prb等于每一个像素值在三个类别上的分布概率
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;