function [E,Edata,Ereg]=Energycolor(f,r,c,sigma,b,RegTerm,G,u,M,N)
K=size(sigma,3);
MN=size(f,1);
V=size(f,2);

B=reshape(b,[M N]);
B2=B.^2;
GB=imfilter(1./B,G,'symmetric');
GB2=imfilter(1./B2,G,'symmetric');
GB=GB(:);
GB2=GB2(:);

for k=1:K
  %   dm=f-repmat(c(k,:),[MN 1]);
  %   E1(:,k)=(0.5*sum(dm.*(dm/sigma(:,:,k)),2)-log(r(k)+eps)+log(sqrt(det(sigma(:,:,k)))+eps)+log(u(:,k)+eps)).*u(:,k);
  Ck=repmat(c(k,:),[MN 1]);
  Temp1=sum(f.*(f/sigma(:,:,k)),2).*GB2-2*sum(Ck.*(f/sigma(:,:,k)),2).*GB+sum(Ck.*(Ck/sigma(:,:,k)),2);
  Temp2=-log(r(k)+eps)+0.5*log(det(sigma(:,:,k))+eps)+log(u(:,k)+eps);
 % assignin('base','Temp1',Temp1);
 % assignin('base','Temp2',Temp2);
  E1(:,k)=(0.5*Temp1+Temp2).*u(:,k);
  RR=RegTerm(:,:,k);
  E2(:,k)=RR(:).*u(:,k);
end
Edata=sum(E1(:))+sum(log(b.^V+eps));
Ereg=sum(E2(:));
E=Edata+Ereg;
%{
K=size(r,3);
[M N]=size(f);
SIGMA=repmat(sigma,[M N 1])+eps;
R=repmat(r,[M N 1])+eps;
F=repmat(f,[1 1 K]);
b=b+eps;
C=repmat(c,[M N 1]);

G_B=repmat(imfilter(1./b,G,'symmetric'),[1 1 K]);
G_B2=repmat(imfilter(1./(b.^2),G,'symmetric'),[1 1 K]);

Temp1=F.^2.*G_B2-2*G_B.*C.*F+C.^2;
Temp1=Temp1./(2*SIGMA);
Temp2=-log(R)+1/2*log(SIGMA)+log(u+eps);
Temp3=(Temp1+Temp2).*u;
Edata=sum(Temp3(:))+sum(log(b(:)));
Temp4=RegTerm.*u;
Ereg=sum(Temp4(:));
E=Edata+Ereg;
%}