Ellipse Fit

Given a set of points (x, y) this function returns the best fit ellipse.

Description

Usage:

[semimajor_axis, semiminor_axis, x0, y0, phi] = ellipse_fit(x, y)

Input:
                    x - a vector of x measurements
                    y - a vector of y measurements

Output:

                  semimajor_axis - Magnitude of ellipse longer axis
                  semiminor_axis - Magnitude of ellipse shorter axis
                  x0 - x coordinate of ellipse center
                  y0- y coordinate of ellipse center
                  phi - Angle of rotation in radians with respect to
                  the x-axis

Algorithm used:

Given the quadratic form of an ellipse:
a*x^2 + 2*b*x*y + c*y^2 + 2*d*x + 2*f*y + g = 0 (1)
we need to find the best (in the Least Square sense) parameters a,b,c,d,f,g.
To transform this into the usual way in which such estimation problems are presented,
divide both sides of equation (1) by a and then move x^2 to the other side. This gives us:
2*b'*x*y + c'*y^2 + 2*d'*x + 2*f'*y + g' = -x^2 (2)
where the primed parametes are the original ones divided by a. Now the usual estimation technique is used where the problem is presented as:
M * p = b, where M = [2*x*y y^2 2*x 2*y ones(size(x))],
p = [b c d e f g], and b = -x^2. We seek the vector p, given by:
p = pseudoinverse(M) * b.
From here on I used formulas (19) - (24) in Wolfram Mathworld:
http://mathworld.wolfram.com/Ellipse.html

Acknowledgements

The author wishes to acknowledge the following in the creation of this submission:
Circle fit
This submission has inspired the following:
Ellipse Fit (Taubin method), Ellipse Fit (Direct method)

• MATLAB release: MATLAB 7.3 (R2006b)

[Nikolai Chernov]

This is an unreliable ellipse fit, as it depends on the choice of the coordinate system (rotating or translating the xy frame would give you a different ellipse). Better fits should be invariant under transformations of coordinates. Besides, this fit may return a hyperbola, instead of ellipse. I have posted two wel-known and well-designed ellipse fits, Taubin Fit and Direct Fit, see files #22683 and #22684.

[Laura Suad]

Hello I am using this ellipse fit program but I need to have the: semimajor_axis, semiminor_axis, x0, y0, phi with dispersion (the fitting errors). How can I do to obtain these parameters?

[Raphaël Cautain]

I tested the function with a true ellipse, successively oriented from 0 to 2pi. All parameters were found
correct, except one :
The orientation was found with a wrong sign each time the axis direction is nearer from Ox than from Oy (that is, 1 time over 2).
Actually, recommendations from Wolfram's pages (results §23) should be - at least - implemented like this :
phi = 0.5 * acot((a-c)/(2*b));
if (a > c)
    phi = phi + pi/2;
end
Rather than like the existing line 77 :
% phi = 0.5 * acot((c-a)/(2*b));
Thus the last lines of ellipse_fit are not adapted and must be removed :
% if (a_prime < b_prime)
% phi = pi/2 - phi;

[Sreedu]

Hello,
At least how many points are required to get a reliable result?