4.0

4.0 | 3 ratings Rate this file 33 Downloads (last 30 days) File Size: 2.64 KB File ID: #22423

Ellipse Fit

by

 

12 Dec 2008 (Updated )

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

| Watch this File

File Information
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

Circle Fit inspired this file.

This file inspired Ellipse Fit (Taubin Method) and Ellipse Fit (Direct Method).

MATLAB release MATLAB 7.3 (R2006b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (8)
07 Jul 2014 Aritra Das

Can anyone provide a clear idea about how to implement this practically?
For example if I want to find ellipses in an binary image say bw, how to run this code on the image to get the ellipses?
As I see there is no way to provide the matrix name as a input argument.
And if somebody could explain the input arguments "Input: x,y - a set of points in 2 column vectors. AT LEAST 5 points are needed !"
I mean this statement a bit elaborately it will be very helpful.

07 Apr 2014 duraid deikran

i have a central vacancy of point after applying certain function, how can i fit ellipse on the periphery of this vacancy using your function

11 Dec 2013 Marcello

Thank you! Great job!

How may I use this function to plot the resulting ellipse?

15 Oct 2013 dingding

thank you ! help me so much!

03 Mar 2011 Sreedu

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

15 Feb 2011 Raphael 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;

09 Apr 2010 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?

15 Jan 2009 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.

Updates
13 Dec 2008

Added input and output explanations to description part.

Contact us