From: Bruno Luong on
This works as well:

A = [bx ax; by ay];
[U S V] = svd(A);
W = U*S; % or W = A*V'

plot(fit_x, fit_y);
hold on
plot([0 W(1,1)],[0 W(2,1)],'r')
plot([0 W(1,2)],[0 W(2,2)],'r')
axis equal

% Bruno