function [r,c,sigma]=MStep(u,f)
[M N K]=size(u);
UK_area=sum(sum(u,1),2)+eps;
r=UK_area./(M*N);
F=repmat(f,[1 1 K]);
c=sum(sum(u.*F,1),2)./UK_area;
C=repmat(c,[M N 1]);
sigma=sum(sum((F-C).^2.*u,1),2)./UK_area;
