function [r,c,sig2,b]=MStep_LGMM(u,f,G,b,M,N)

K=size(u,2);
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');
GBv=repmat(GB(:),[1,V]);
GB2v=repmat(GB2(:),[1 V]);

for k=1:K
    uk=u(:,k);
    UK_area=sum(uk)+eps;
    r(k)= UK_area/MN;
    ukv=repmat(uk,[1 V]);
    c(k,:)=sum(f.*GBv.*ukv)/UK_area;
    Ck=repmat(c(k,:),[MN 1]);
    sig2(:,:,k)=(f'*((f.*GB2v-Ck.*GBv).*ukv)+Ck'*((-f.*GBv+Ck).*ukv))/UK_area+eps*eye(V);
    Temp1=reshape(sum(Ck.*(f/sig2(:,:,k)),2),[M N]);
    Temp2=reshape(sum(f.*(f/sig2(:,:,k)),2),[M N]);
    UK=reshape(uk,[M N]);
    S(:,:,k)=imfilter(UK.*Temp1,G,'symmetric');
    T(:,:,k)=imfilter(UK.*Temp2,G,'symmetric');
end

s=sum(S,3);
s=s(:);
t=sum(T,3);
t=t(:);
b=(-s+sqrt(s.^2+(4*V)*t))/(2*V);
