
clc;
clear;
Ima=imread('5.bmp');
%Ima=imread('C7.jpg');
%Ima=imread('91.png');
%Ima=imread('86016.jpg');
%Ima=imread('105053.jpg');


Ima=double(Ima);%Ima=255*mat2gray(Ima);
[M N V]=size(Ima);
MN=M*N;
for i=1:V
    Temp=Ima(:,:,i);
    f(:,i)=Temp(:);
end


K=2;
Maxiter=1000;
epsilon=1e-9;
wsize=7;
sig=10.0;
G=fspecial('gaussian',round(2*sig)*2+1,sig);     


%{
figure;imshow(uint8(Ima));hold on;
text(6,6,'Left click to get points, right click to get end point','FontSize',[12],'Color', 'r');
L=repmat(0,[M N]);
for i=1:K-1
Mask=roipoly;
L(Mask==1)=i;
end
L=double(L)+1;
L=L(:);
%}

[L kc]=kmeans(f,K);

color=['r','y','b','c'];
resnum=0;
figure(1),imshow(uint8(Ima),[]);
hold on;
for i=1:K-1
contour(reshape(L==i,[M N]),[0.5 0.5],color(i),'Linewidth',2);
end
hold off;
resnum=resnum+1;
Liu=getframe(gca);
Liures{resnum}=frame2im(Liu);

for k=1:K
    uk=double(L==k);
    u(:,k)=uk;
    uk_area=(sum(uk)+eps);
    c(k,:)=sum(f.*repmat(uk,[1 V]))./uk_area;
    dm=f-repmat(c(k,:),[MN 1]);
    sig2(:,:,k)=(dm'*(dm.*repmat(uk,[1 V])))./uk_area+eps*eye(V);
    r(k)=1/K;
end
U=reshape(u,[M N K]);
b=repmat(1,[MN 1]);
oldJ=inf;


for k=1:K
    c1(:,:,k)=c(k,:);
    r1(:,:,k)=r(k);
end
U=reshape(u,[M N K]);
RegTerm=Reg_term(U,wsize);
for i=1:Maxiter
     
       if i==1
     
         [J1,Edata1,Ereg1]=Energycolor(f,r,c,sig2,b,RegTerm,G,u,M,N);
           beta=Edata1/(Ereg1*40);
       end
       
       u=EStep_LGMM(f,r,c,sig2,b,beta*RegTerm,G,M,N);
     U=reshape(u,[M N K]);
      [r,c,sig2,b]=MStep_LGMM(u,f,G,b,M,N); 
       RegTerm=Reg_term(U,wsize); 
       [J(i),Edata(i),Ereg(i)]=Energycolor(f,r,c,sig2,b,beta*RegTerm,G,u,M,N);
       newJ=J(i);
       if mod(i,1)==0
       figure(1)
       %subplot(221),
       imshow(uint8(Ima),[]); title(['iter=' num2str(i)]);
       hold on
       [maxv,segres]=max(U,[],3);
       for icor=1:K-1
       contour(segres==icor,[0.5,0.5],color(icor),'Linewidth',2);
       end
       hold off;
       resnum=resnum+1;
       Liu=getframe(gca);
       Liures{resnum}=frame2im(Liu);
       end
        
       if (oldJ-newJ).^2<epsilon*oldJ.^2
       [Jend,Edataend,Eregend]=Energycolor(f,r,c,sig2,b,RegTerm,G,u,M,N);
       beta=Edataend/Eregend;
       u=EStep_LGMM(f,r,c,sig2,b,beta*RegTerm,G,M,N);
       U=reshape(u,[M N K]);
       figure(1)
       %subplot(221),
       imshow(uint8(Ima),[]); title(['iter=' num2str(i)]);
       hold on
       [maxv,segres]=max(U,[],3);
       for icor=1:K-1
       contour(segres==icor,[0.5,0.5],color(icor),'Linewidth',2);
       end
       hold off;
       resnum=resnum+1;
       break;
       end
       oldJ=newJ;
       
end

[maxv,segres]=max(U,[],3);
figure,
subplot(221),imshow(uint8(Ima));hold on,contour(segres,[1 1],'r'),contour(segres,[3 3],'b');
title('segmentation');
subplot(222),imshow(segres,[]);
title('label');
subplot(223),imshow(reshape(b,[size(Ima,1) size(Ima,2)]),[])
title('background');
 res=repmat(0,[M N V]);
for k=1:K
    cor=repmat(reshape(c(k,:),[1 1 V]),[M N 1]);
    ind=find(segres==k);
    res1=cor.*repmat(double(segres==k),[1 1 V]);
    res=res+res1; 
end

