File Exchange

image thumbnail

fit_ellipse

version 1.0.0.0 (4.19 KB) by Ohad Gal
Find the best fit for an ellipse using a given set of points (a closed contour).

107 Downloads

Updated 02 Oct 2003

View License

This function uses the Least-Squares criterion for estimation of the best fit to an ellipse from a given set of points (x,y). The LS estimation is done for the conic representation of an ellipse (with a possible tilt).

Conic Ellipse representation = a*x^2+b*x*y+c*y^2+d*x+e*y+f=0
(Tilt/orientation for the ellipse occurs when the term x*y exists (i.e. b ~= 0))

Later, after the estimation, the tilt is removed from the ellipse (using a rotation matrix) and then, the rest of the parameters which describes an ellipse are extracted from the conic representation.

For debug purposes, the estimation can be drawn on top of a given axis handle.

Note:
1) This function does not work on a three-dimensional axis system. (only 2D)
2) At least 5 points are needed in order to estimate the 5 parameters of the ellipse.
3) If the data is a hyperbola or parabula, the function return empty fields and a status indication

Comments and Ratings (99)

Matthias

In principle the function works fine. As others have pointed out, the angle phi which is returned is somewhat inconsistent. It would make much more sense to assign to angle to one of the two axis of the ellipses and count from 0 to 180 degree so the angle is unambiguous. Also, I'm not sure what the values of a and b which are returned are supposed to represent. An illustration would be helpful.

what is axis_handle? In which format we have to give axis_handle

Hongyu Zhao

Mengqi Li

Great function! works well and fast and has very practical output. Thanks!

olimpia_m

PEI ZHAO

harsh jain

Works very well,
For me I created another output which is
angle_Xaxis_to_long_ellipse_axis = 180 - (1/2 * atan2( b,(c-a) ))*180/pi; % It's not radians but degrees
It works for my data

lone Ishfaq

Hello, can anyone help me. i am not getting plot after the execution.
also not getting any error

Hello. Can someone please help me with the following error I got while running the above uploaded code by @Ohad Gal
Error using ./
Matrix dimensions must agree.

Error in uFit_ellipse (line 157)
a = numer./denom;

lucas lenne

Xiong Zhang

Can anyone help add a line to plot the fitted ellipse ?

Hello, The fit works very well. Thank you for that.

I have a question regarding the tilt of the ellipse. If, I am not wrong the tilt of the ellipse is in the clockwise direction from the negative x-axis.

Did anybody use these fits in correlation with Lissajous Curves?

Ali Buendia

Great example.

Luuk

Very useful function. Any insight on how to modify the script such that it finds a weighted least-squares solution? Currently if the input data contains one outlier (which for some reason is not detected as an outlier in MATLAB), the fitted ellipse is affected a lot. See this example: https://imgur.com/a/nwQPXp0

Any suggestions?

Edit: Added an example

Excellent function thanks ! Also the comments are great to understand how it works !
Also thanks Frank for the fix on the angle problem.
As far as I could see the two main things to be wary of here are:
- the uncertainty that atan generates (the double angle is calculated modulo pi so the angle orientation_rad is calculated modulo pi/2) -> this is resolved by using atan2 as Frank suggests below
- the fact that in your calculations you use the rotation matrix as R = [cos(phi) sin(phi); -sin(phi) cos(phi)] when (as far as I could see) usually rotations are defined as R = [cos(phi) -sin(phi); sin(phi) cos(phi)]. This is something to take into account if you wish to make your own function to plot the ellipse; otherwise you will indeed need to use theta= pi-phi as the tilt angle.
Thanks a lot for this very useful code !!!

Hi thanks for this code! but How can I plot the ellipse and the its fit in matlab? becouse i just get output data without plot

BlueEyes

Example carried out _without_ the 'mysterious' parameter (axis_handle), by commenting rows 244,262,263 where is mentioned:
----------------
x= [0;2;6;10;12;8;4];
y= [2;6;10;10;4;0;-2];
fit_ellipse(x,y)
----------------
The ellipse is properly fit. Please, author, be kind to specify. Thx

jacopo

Hi very nice code! How can I plot the ellipse and the its fit in matlab?? Please help me!

It would be nice if the code output some measure of the error in the parameters.

I don't understand what is wanted for an axis_handle to enable the ellipse to be plotted

i have no idea what is center points between X0 Y0 or X in Y in. on my task, those two points are very different. please let me know..

Note; internal calculations of angleFromX appear to be out by 90 degrees (wrong quadrant). Great script though, extremely helpful.

Excellent matlab script to fit an ellipse to a irregular shaped data.... it would be nice if the code can be extended to plotting multiple ellipse for scattered data.

Excellent matlab script to fit an ellipse to a irregular shaped data.... it would be nice if the code can be extended to plotting multiple ellipse for scattered data.

Hi Thanks for the Code it works great! I am trying to find the area of intersection between 2 ellipses. Is there anyway i can determine that with an extension of this script? Thank you!

plasmageek

I was using a different ellipse fit code, and it was getting hung up on a more complicated data set. This one worked great and is not memory intensive.

Great submission! One small issue, when working with multiple axes and trying to plot into them - the way the plot commands are set up will plot to the active axes. As the user passes axis_handle to the funciton anyway, it'd be useful to explicitly plot into that axis.

Ho King

It is fitting the ellipse but not least square method. I really want to know the principle for your method to calculate about a b c d e f.

Is there any document?

Oron Nir

Works great but is there a paper to backup these equations?
Thanks!

rcjr15

I get the following error:
In fit_ellipse at 155
I gave x=[191 192 193 194 195]'
and y=[56 57 58 59 60]

Setoshi C

About orientation problem, angle_to_x = pi-phi instead of phi.

Rholais Lii

WARNING('') does no reset the warning state. Use LASTWARN('') instead.

Frank

One more detail, with the atan2() fix, the check for orientation_tolerance must be removed.

I believe there is an error on the code comments (not the code itself):
< A'*X'*X*A + 2*f_c'*X*A + N*f^2
> A'*X'*X*A - 2*f_c'*X*A + N*f^2
(should the sign of that factor of the cost function not be negative?)

By the way, it would be cool for the ellipse fit structure to feature the MSE error (so we know how good the fit was):
N = length(x);
MSE = (a*X'*X*a'-2*sum(X*a')+N)/N;

Frank

Excellent function, but I also had the same phase issue that Rob (below) had. The fix is line 170, change from:

orientation_rad = 1/2 * atan( b/(c-a) );

to:

orientation_rad = 1/2 * atan2( b, (c-a) );

so the returned value covers +/- 90 deg instead of only +/- 45 deg

Rob Barton

% Add my code to get the true angle of rotation of the long axis
% form the quadratic matrix
Q = [a b/2; b/2 c];

% perform eigen Decomp on it
[eigVec, eigValue] = eig(Q);

% compute the angle to the long axis
angleToX = atand(abs(eigVec(2,1))/abs(eigVec(1,1)));

% check sign to get angles from 90-180
% since vector could point other way have to check all 4 quadrants
if ( sign(eigVec(1,1)) == -1 )
% long axis points to the left
if (sign(eigVec(2,1)) == 1 )
% points up so leave in first 0-90 quadrant
angleFromX = angleToX;
else
% Points down, so will treat it as an ellipse with the long
% axis in the 90-180 range since axis can point either way.
angleFromX = 90+(90-angleToX);
end
else
% long axis points to the right
if ( sign(eigVec(2,1)) == 1 )
% long axis point up
angleFromX = 90+(90-angleToX);
else
% long axis points down
angleFromX = angleToX;
end
end

Here is the final addition I would suggest to get the true angle (angleFromX) (the full 360 degrees)

then just add:

'angleToX', angleToX, ...
'angleFromX', angleFromX, ...

to your output structure.

I put this code right after your call to

% extract parameters from the conic equation
[a,b,c,d,e] = deal( a(1),a(2),a(3),a(4),a(5) );

but you can organize it however is best.

Rob Barton

After working through your well written comments, and some web searching I figured out what you are doing, and would suggest to perhaps add the following.

Perhaps add another member in the return structure say 'angleToX' and then perhaps add some code right after yourorientation_rad to compute the angle or roatation from the long axis to the X axis as:

Q = [ a b/2; b/2 c];
[eigVec, eigValue] = eig(Q);
angleToX = atand(abs(eigVec(2,1)/abs(eigVec(1,1)));

Q just puts the Ax^2+Bxy+Cy^2 into quadratic form, and the Eigen decomposition vectors are the two axis of the ellipse. by computing the tangent we effectively compute the angle from the long axis to the x axis, so we know how much the ellipse is rotated relative to the X axis. In this manner you will get values greater than 45 degrees.

The phi you compute is for the purpose of finding how to zero out the b element.

If you add this then you can use a 'draw ellipse' function to draw the ellipse you just fit, and it will know how much to rotate it to overlay the scatter points you just fit.

Thanks very much, excellent code.

Rob Barton

I would like to use this function to fit ellipses at all angles. Could you explain how the orientation is being computed? It appears after 45 degrees it is not working as I would anticipate.

Looking at your code I can see you are using a tangent...

could you explain a bit... Can this routine be used for ellipses at more than 45 degrees?

Lucy

Lucy

I can't plot the ellipse. Could please anyone give me an explanation on how to use the code? Thanks

Christian

For my tasks this function works better than the others I've tried so far. Thank you very much Ohad!!

@Aritra: To run this function on a binary image you have to run:
[X Y] = ind2sub(size(img),find(img));
E = fit_ellipse(x,y);

then you can do:
if E.long_axis > 0
[X, Y] = calcEllipse(E, 360);
end

To plot it the ellipse:
plot(Y, X);

The function calcEllipse:
function [X,Y] = calcEllipse(varargin)
% function [X,Y] = calculateEllipse(x, y, a, b, angle, steps)
%# This functions returns points to draw an ellipse
%#
%# @param x X coordinate
%# @param y Y coordinate
%# @param a Semimajor axis
%# @param b Semiminor axis
%# @param angle Angle of the ellipse (in rad)
%#
% Source: http://stackoverflow.com/questions/2153768/draw-ellipse-and-ellipsoid-in-matlab/24531259#24531259
% Modified by Christian Fässler

steps = 360;

if nargin == 1 || nargin == 2
x = varargin{1}.X0_in;
y = varargin{1}.Y0_in;
a = varargin{1}.a;
b = varargin{1}.b;
angle = varargin{1}.phi;
if nargin == 2
steps = varargin{2};
end
else if nargin == 5 || nargin == 6
x = varargin{1};
y = varargin{2};
a = varargin{3};
b = varargin{4};
angle = varargin{5};
if nargin == 6
steps = varargin{6};
end
else
error('Wrong input');
end
end

beta = -angle;
sinbeta = sin(beta);
cosbeta = cos(beta);

alpha = linspace(0, 2*pi, steps)';
sinalpha = sin(alpha);
cosalpha = cos(alpha);

X = round(x + (a * cosalpha * cosbeta - b * sinalpha * sinbeta));
Y = round(y + (a * cosalpha * sinbeta + b * sinalpha * cosbeta));

if nargout==1, X = [X Y]; end
end

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.

Marcello

Great code, very useful!
However, I found a bug related to the orientation of the ellipse and I read some people had the same problem. The author wrote which "to correct that, the test for the orientation_tolerance should be normalized".
How can I solve the problem? May anyone help me? I would be very grateful.
Thanks in advantage, and sorry for my English

dingding

thank you very much! very good work! help me so much!

Shuqing

Thank you very much! Very helpful.

Amir

many thanks !!

Great code! Works great...I'll have to spend some time with it. For my problem, I need to force the origin to 0,0 and orthogonal axes.

I'm n0t able to plot the output results..Where do i go wrong?
%%
h=plot(ydata,zdata);
ellipse_t = fit_ellipse(ydata,zdata,h)

Gilad Kapel

Amit Ruf

Nice code, however, you probably have a bug related to the orientation of the ellipse, see:

http://www.mathworks.com/matlabcentral/fileexchange/22423

I fixed it and added a computation of the residual of the fit which provides a quality measure for the fit.
If you want my version, please mail to:
amitruf@gmail.com

Paulo

How to use this? I mean, the inputs x and y? Can I have an example? Thanks! This is very useful for my project!

Works out of the box. Very useful. Thanks.

Excellent deployment, easy to use and well commented ... THANKS

Works Great,
Thanks for sharing...

Tima Tima

Thank you. It is vary helped for my work

Thanos

Hi, I am using the script and it works fine for my data. I am trying to work the math a little bit and I am struggling at one point. why we assume that the f=-1? As far as I can see from the code, the value we assign to f affects the a, b, c, d, e which in turn affect the ellipse properties. Is that something mathematically trivial? I am looking forward for your answers
Cheers
Thanos

bear tiger

Rafal

Hello. Can someone help me and give description of algorith used in this matlab code. Maybe a link to publication. If someone has also block scheme of this algorithm send it to me please. My mail :
Monter70@gmail.com

Best Regards

Rafal

Hello i use this script and it's working perfectly,but unfortunately i have one problem ...

I load signals ,it draws an ellipse , when i want to put near ellipse points to which there is aproximation they are moved. I think it could be connected with new coordinates for ellipse. cause i use cmd plot(x,y,'b'); and there is ideal bow as ellipse has,but it's moved. Could anyone helps me?

Rakesh

Thanks a lot!! nice work!

Roy

thanks for sharing.

Rafal

Hello. I get all the ellipse parameters.I load two orthogonal signals as a x=b1ch3(:,1), y=b1ch4(:,1), there is problem because it's not drawing and fitting ellipse to this points.

scaramanga

I'm using your function in a loop, is it possible to draw ellipses with a different color at every iteration of the loop ?

Mustafa

Evgeny Pr

Sieun

thank you :D

Thanks for your sharing.

Sophie

I found.. it's just a question of angle: should have done:
t = - ellipse_t.phi;

Sophie

@Samuel: did you resolve your problem?
I do not know your exact problem but maybe what you want is this:
handle= subplot(221); % or something like this
ellipse_t = fit_ellipse(x,y,handle); % the ellipse should be drawn on the subplot else your ellipse is out of range if I am not wrong

Sophie

Thanks a lot for this script which is really useful!
I have one question because I get into trouble: I would like to get all pixels that are inside the ellipse.. I was doing something like that but it does not work:

x0 = ellipse_t.X0_in;
y0 = ellipse_t.Y0_in;
a = ellipse_t.a;
b = ellipse_t.b;
t = ellipse_t.phi;
for x=1:size(I,1)
for y=1:size(I,2)
X = (x-x0)*cos(t)+(y-y0)*sin(t);
Y = -(x-x0)*sin(t)+(y-y0)*cos(t);
if (X^2/a^2+Y^2/b^2)>1 % outside ellipse
I(y,x) = 0;
end
end
end

Any idea? Thanks in advance!

Samuel

Hi, I have one question.... When a run the code the program shows the ellipse result, but don't plot the graphic with the points and ellipse fit curve (the same graphic above). I would like to see the fit curve, but I don't know whats is the problem... Please, anybody help-me.... (sorry my english)

Kevin Shaw

Love it. Works as advertised.

Dave Peake

Exactly what I was after, perfect. Does what it says, does it well and easily.

Daniel Nilsson

Awesome! I just had to modify the program to plot the non-rotated ellipse instead, otherwise this was perfect!

a s

Roy Pur

THANKS

Raimund Leitner

Incredible useful and practical m-file

Sophie Jarlier

It's grest to have your code! Just a question: I randomly selected 5points and execute your function but I get the following error:
In fit_ellipse at 155
stopped because of a warning regarding matrix inversion
But I have 4 points that have same y=320 which is the maximum size of my image.. Should these 5points have 5different y and 5different x?
Thanks for your answer
Sophie

Heikki Suhonen

Works well for me. Easy interface, and nice optional visualization.

Willem Van der Merwe

Great. The plot function was a great help as well.

Gita Adur

Execellent!!! Very well written and extremely easy to use.

Antje _

Thanks!

Eric Tittley

Chaiyanan Sompong

Thank you. In present-day, i do my research about face detection in color image.

P F

Does what it claims to. Quite helpful!

Peter Delahunt

Thanks Ohad - it works well with color discrimination data. Best wishes, Peter

Larry O'Neill

This method works excellent with all my noisy data. Also well written and extremely easy to use.

Updates

1.0.0.0

1. added a test to identify if the data is a hyperbola or parabola - returned in the "status" field
2. the routine finds now the center point of the original (tilted) ellipse as well (fields "X0_in","Y0_in")

for ellipse with large axes, parameters "a" and/or "b", the fitting might not locate the orientation of the ellipse, especially if the angle is small.

to correct that, the test for the orientation_tolerance should be normalized.

MATLAB Release Compatibility
Created with R12.1
Compatible with any release
Platform Compatibility
Windows macOS Linux