Andrew W. Fitzgibbon, Maurizio Pilu, and Robert B. Fisher
Direct least-squares fitting of ellipses,
IEEE Transactions on Pattern Analysis and Machine Intelligence, 21(5), 476--480, May 1999
It's important to scale the image coordinates to be near one before running the algorithm or you'll get numerical instability. See the improved matlab code for an implementation which includes this step.
If you're implementing it in a language other than matlab, there's a better way to solve the eigenvalue problem than we knew when we wrote the paper, which involves dividing the matrix into blocks.
[A B] [x] = l [D 0] [x] [B' C] [y] [0 0] [y]is the same as the pair of equations
A x + B y = l D x B' x + C y = 0 ==> y = -inv(C) B' xso the first eqn is
A x + B (-inv(C) B') x = l D xor
inv(D) * (A - B*inv(C)*B') x = l xWhat this means is that you can solve the system by breaking the original matrices up into the blocks A,B,C,D and writing
[allx, allD] = eig(inv(D) * (A - B*inv(C)*B')) x = "column of allx corresponding to the positive eigenvalue"; y = -inv(C) * B' * x;