% line=linefit(points);
%
% points = [x1 x2 ... xn; y1 y2 ... yn];

function line=linefit(points);

n=size(points,2);

x=points(1,:);
y=points(2,:);

xm=mean(x); x2m=mean(x.*x);
ym=mean(y); y2m=mean(y.*y);
xym=mean(x.*y);

A=[x2m-xm*xm    xym - xm*ym; 
    xym-xm*ym   y2m-ym*ym];
[V,D]=eig(A);
ab1=V(:,1);
c1=-ab1(1)*xm-ab1(2)*ym;
ab2=V(:,2);
c2=-ab2(1)*xm-ab2(2)*ym;

line1=[ab1;c1];
line2=[ab2;c2];

res1=line1'*[points; ones(1,n)];
res2=line2'*[points; ones(1,n)];

if (res1*res1'<res2*res2') 
    line=line1;
    residual=sqrt(res1*res1'/n)
else
    line=line2;
    residual=sqrt(res2*res2'/n)
end