%For experiment 1 in the paper.
clc;
clear;
[Ima, groundtruth]=creatdataforexp1();
Ima=mat2gray(Ima);
randn('seed',100);
f=Ima+(40/255)*randn(size(Ima));
f=f*255;
[M N]=size(f);
K=4;
Maxiter=100;
epsilon=1e-9;
wsize=3;

%{
r(:,:,1:K)=1/K;
L=repmat(1,size(f));
L(47:60,86:99)=2;
for k=1:K
    uk=L==k;
    uk_area=(sum(uk(:))+eps);
    c(:,:,k)=sum(sum(f.*uk))/uk_area;
    sigma(:,:,k)=sum(sum((f-repmat(c(k),[M N])).^2.*uk))/uk_area;
 %c(:,:,k)=1/k; 
% sigma(:,:,k)=1;
end
u=repmat(1/K,[M N K]);
%}

%{
Maxf=max(f(:));
minf=min(f(:));
r=reshape(ones(1,K)/K,[1 1 K]);
c=reshape((1:K)/(1+K)*(Maxf-minf)+minf,[1 1 K]);
sigma=reshape(ones(1,K)*Maxf,[1 1 K]);
%}

%%Kmean initial.
[L kc]=kmeans(f(:),K,'start',[0;255*1/3;255*2/3;255]);
L=reshape(L,size(f));
[kc ind]=sort(kc);
KL=repmat(0,size(L));
for k=1:K
    uk=double(L==ind(k));
    KL(L==ind(k))=k;
    u(:,:,k)=uk;
    uk_area=(sum(uk(:))+eps);
    c(:,:,k)=sum(sum(f.*uk))/uk_area;
    sigma(:,:,k)=sum(sum((f-repmat(c(k),[M N])).^2.*uk))/uk_area;
    r(:,:,k)=1/K;
end

p3=SA(reshape(KL,size(groundtruth)),groundtruth);
%[fy,fx]=gradient(imfilter(f,fspecial('gaussian',[15 15],2.0),'symmetric'));
%g1=repmat(1./(1+sqrt(fx.^2+fy.^2)),[1 1 K]);
%g1=repmat(1,[M N K]);
tic;
oldE=inf;
for i=1:Maxiter
       RegTerm=Reg_term(u,wsize);
       %Set the regularization parameter \lambda.
       if i==1
       [E(i),Edata(i),Ereg(i)]=Energy_Total(f,r,c,sigma,RegTerm,u);
       beta=abs(Edata(i))/abs(Ereg(i)+eps);
       end
       u=EStep_smooth(f,r,c,sigma,beta*RegTerm);
       [r,c,sigma]=MStep(u,f);
       [E(i) Edata(i) Ereg(i)]=Energy_Total(f,r,c,sigma,beta*RegTerm,u);
       newE=E(i);
       if (oldE-newE).^2<epsilon*oldE.^2
           break;
       end
       oldE=newE;
end
t1=toc;

tic;
oldLL=-inf;
for i=1:Maxiter
        u1=EStep(f,r,c,sigma);
       [r,c,sigma]=MStep(u1,f);
       newLL=logLikehood(f,r,c,sigma);
%       if mod(i,5)==0
%       fprintf(['iter=' num2str(i) ', LL=%f\n'],newLL);
 %      end
       if (newLL-oldLL).^2<epsilon*oldLL.^2;
           break;
       end
       oldLL=newLL;
end
t2=toc;
[maxv,segres1]=max(u,[],3);
[maxv,segres2]=max(u1,[],3);
p1=SA(segres1,groundtruth);
p2=SA(segres2,groundtruth);
figure,
subplot(221),imshow(f,[]);
subplot(222),imshow(segres2,[]);title('EM');
subplot(223),imshow(segres1,[]);title('Regularized EM');

xmin=1;
xmax=length(E);
%ymin=min(min(Edata),min(Ereg));
%ymax=;max(max(Edata),max(Ereg));
figure('position',[100 200, 600 800]),plot(E,'r*-','linewidth',4);%hold on;
set(gca,'fontsize',36);
axis([xmin xmax min(E)-(max(E)-min(E))/4 max(E)]);
% Z=getframe(gcf);
% z=frame2im(Z);
% imwrite(z,'exp1_J.jpg');
figure('position',[100 200, 600 800]),plot(Edata,'b*-','linewidth',4);
set(gca,'fontsize',36);
axis([xmin xmax min(Edata) max(Edata)+(max(Edata)-min(Edata))/4]);
% Z=getframe(gca);
% z=frame2im(Z);
% imwrite(z,'exp1_E.jpg');
figure('position',[100 200, 600 800]),plot(Ereg,'g*-','linewidth',4);
set(gca,'fontsize',36);
axis([xmin xmax min(Ereg)-(max(Ereg)-min(Ereg))/4 max(Ereg)]);
% Z=getframe(gca);
% z=frame2im(Z);
% imwrite(z,'exp1_R.jpg');



