function [e a b x0 y0] = smc(rho, xm, ym)
rho = rho/sum(rho(:));

if nargin < 3
    [ym,xm] = ndgrid(1:size(rho,1),1:size(rho,2));
end

%s_rho = size(rho);
%[x,y] = meshgrid(1:s_rho(2),1:s_rho(1));

x0 = sum(sum(rho.*xm));
y0 = sum(sum(rho.*ym));
 
%x_avg = sum(sum(x.*rho));
%y_avg = sum(sum(y.*rho));

%xx = x - x_avg;
%yy = y - y_avg;

xx = xm-x0;
yy = ym-y0;

x_rms = sqrt(sum(sum(xx.^2.*rho)));
y_rms = sqrt(sum(sum(yy.^2.*rho)));
 
xy_cor = sum(sum(xx.*yy.*rho));
xy = xy_cor/x_rms^2;
 
e = sqrt(x_rms^2*y_rms^2-xy_cor^2);
a = -xy_cor/e;
b = x_rms^2/e;
end