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).

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

Ohad Gal (2021). fit_ellipse (https://www.mathworks.com/matlabcentral/fileexchange/3215-fit_ellipse), MATLAB Central File Exchange. Retrieved .

Created with
R12.1

Compatible with any release

**Inspired:**
BiofilmQ, Corneal Topography: Constructing Curvature Topography from Placido Rings Image, Elliptical Scanning Algorithm, Drop shape analysis. Fit contact angle by double ellipses or polynomials, Erosion Metrics for conical landforms App

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Rebaz Aligoksenin hande bayazitwell explained, works great. thanks a lot!

Nathaniel PlesniarskiWenfang ZhangQian XuGEUNGYU JEONGHenryMatthiasIn 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.

Shashi Kumarwhat is axis_handle? In which format we have to give axis_handle

Hongyu ZhaoMengqi LiDano RoelvinkGreat function! works well and fast and has very practical output. Thanks!

olimpia_mPEI ZHAOharsh jainOscar MALLETWorks 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 IshfaqHello, can anyone help me. i am not getting plot after the execution.

also not getting any error

SAI MANOJ KONDAPALLIHello. 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 lenneXiong ZhangCan anyone help add a line to plot the fitted ellipse ?

Shashank GuptaHello, 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 BuendiaGreat example.

Narrendar RaviChandranVan Linh NgoLuukVery 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

Jeanne StranskyExcellent 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 !!!

Milad Momenisaifalden altaieHi 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

BlueEyesExample 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

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

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

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

Min-Jong Parki 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..

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

sangeetha sugumarExcellent 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.

sangeetha sugumarExcellent 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.

Ganeshkumar MHi 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!

plasmageekI 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.

Kristjan PoderGreat 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.

Abhinna beheraHo KingIt 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 NirWorks great but is there a paper to backup these equations?

Thanks!

rcjr15I 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 CAbout orientation problem, angle_to_x = pi-phi instead of phi.

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

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

Ricardo CruzI 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;

FrankExcellent 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 BartonAfter 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 BartonI 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?

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

ChristianFor 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 DasCan 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.

MarcelloGreat 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

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

ShuqingThank you very much! Very helpful.

AmirLuigi Sanguignomany thanks !!

Martina CallaghanJeff AndersonGreat 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.

Jeff AndersonNavneet ViswanI'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 KapelAmit RufNice 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

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

David ScadutoSteven DakinWorks out of the box. Very useful. Thanks.

Hassan NaseriExcellent deployment, easy to use and well commented ... THANKS

Mohammed El-SaidWorks Great,

Thanks for sharing...

Tima TimaThank you. It is vary helped for my work

ThanosHi, 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 tigerRafalHello. 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

RafalHello 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?

RakeshThanks a lot!! nice work!

Roythanks for sharing.

RafalHello. 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.

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

MustafaEvgeny PrSieunthank you :D

Raymond ChengThanks for your sharing.

SophieI 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

SophieThanks 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!

SamuelHi, 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 ShawLove it. Works as advertised.

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

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

a sRoy PurTHANKS

Raimund LeitnerIncredible useful and practical m-file

Sophie JarlierIt'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 SuhonenWorks well for me. Easy interface, and nice optional visualization.

Willem Van der MerweGreat. The plot function was a great help as well.

Gita AdurExecellent!!! Very well written and extremely easy to use.

Antje _Thanks!

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

P FDoes what it claims to. Quite helpful!

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

Larry O'NeillThis method works excellent with all my noisy data. Also well written and extremely easy to use.