This function uses the LeastSquares 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 threedimensional 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
1.0.0.0  1. added a test to identify if the data is a hyperbola or parabola  returned in the "status" field


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. 
Inspired: Corneal Topography: Constructing Curvature Topography from Placido Rings Image, Elliptical Scanning Algorithm, Drop shape analysis. Fit contact angle by double ellipses or polynomials
Create scripts with code, output, and formatted text in a single executable document.
Matthias (view profile)
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.
Shashi Kumar (view profile)
what is axis_handle? In which format we have to give axis_handle
Hongyu Zhao (view profile)
Mengqi Li (view profile)
Dano Roelvink (view profile)
Great function! works well and fast and has very practical output. Thanks!
olimpia_m (view profile)
PEI ZHAO (view profile)
harsh jain (view profile)
Oscar MALLET (view profile)
Works very well,
For me I created another output which is
angle_Xaxis_to_long_ellipse_axis = 180  (1/2 * atan2( b,(ca) ))*180/pi; % It's not radians but degrees
It works for my data
lone Ishfaq (view profile)
Hello, can anyone help me. i am not getting plot after the execution.
also not getting any error
SAI MANOJ KONDAPALLI (view profile)
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 (view profile)
Xiong Zhang (view profile)
Can anyone help add a line to plot the fitted ellipse ?
Shashank Gupta (view profile)
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 xaxis.
Did anybody use these fits in correlation with Lissajous Curves?
Ali Buendia (view profile)
Great example.
Narrendar RaviChandran (view profile)
Van Linh Ngo (view profile)
Luuk (view profile)
Very useful function. Any insight on how to modify the script such that it finds a weighted leastsquares 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 Stransky (view profile)
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= piphi as the tilt angle.
Thanks a lot for this very useful code !!!
Milad Momeni (view profile)
saifalden altaie (view profile)
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 (view profile)
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
Wojciech Matwij (view profile)
jacopo (view profile)
Hi very nice code! How can I plot the ellipse and the its fit in matlab?? Please help me!
Samuel Lazerson (view profile)
It would be nice if the code output some measure of the error in the parameters.
Heather Morgan (view profile)
I don't understand what is wanted for an axis_handle to enable the ellipse to be plotted
MinJong Park (view profile)
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..
Frederic Moisy (view profile)
Jessica Hiscocks (view profile)
Note; internal calculations of angleFromX appear to be out by 90 degrees (wrong quadrant). Great script though, extremely helpful.
sangeetha sugumar (view profile)
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.
sangeetha sugumar (view profile)
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.
Ganeshkumar M (view profile)
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 (view profile)
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.
Kristjan Poder (view profile)
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.
Abhinna behera (view profile)
Ho King (view profile)
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 (view profile)
Works great but is there a paper to backup these equations?
Thanks!
rcjr15 (view profile)
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 (view profile)
About orientation problem, angle_to_x = piphi instead of phi.
Rholais Lii (view profile)
WARNING('') does no reset the warning state. Use LASTWARN('') instead.
Frank (view profile)
One more detail, with the atan2() fix, the check for orientation_tolerance must be removed.
Ricardo Cruz (view profile)
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 (view profile)
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/(ca) );
to:
orientation_rad = 1/2 * atan2( b, (ca) );
so the returned value covers +/ 90 deg instead of only +/ 45 deg
Rob Barton (view profile)
% 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 90180
% 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 090 quadrant
angleFromX = angleToX;
else
% Points down, so will treat it as an ellipse with the long
% axis in the 90180 range since axis can point either way.
angleFromX = 90+(90angleToX);
end
else
% long axis points to the right
if ( sign(eigVec(2,1)) == 1 )
% long axis point up
angleFromX = 90+(90angleToX);
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 (view profile)
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 (view profile)
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 (view profile)
Lucy (view profile)
I can't plot the ellipse. Could please anyone give me an explanation on how to use the code? Thanks
Christian (view profile)
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/drawellipseandellipsoidinmatlab/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 (view profile)
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 (view profile)
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 (view profile)
thank you very much! very good work! help me so much!
Shuqing (view profile)
Thank you very much! Very helpful.
Amir (view profile)
Luigi Sanguigno (view profile)
many thanks !!
Martina Callaghan (view profile)
Jeff Anderson (view profile)
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.
Jeff Anderson (view profile)
Navneet Viswan (view profile)
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 (view profile)
Amit Ruf (view profile)
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 (view profile)
How to use this? I mean, the inputs x and y? Can I have an example? Thanks! This is very useful for my project!
David Scaduto (view profile)
Steven Dakin (view profile)
Works out of the box. Very useful. Thanks.
Hassan Naseri (view profile)
Excellent deployment, easy to use and well commented ... THANKS
Mohammed ElSaid (view profile)
Works Great,
Thanks for sharing...
Tima Tima (view profile)
Thank you. It is vary helped for my work
Thanos (view profile)
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 (view profile)
Rafal (view profile)
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 (view profile)
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 (view profile)
Thanks a lot!! nice work!
Roy (view profile)
thanks for sharing.
Rafal (view profile)
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 (view profile)
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 (view profile)
Evgeny Pr (view profile)
Sieun (view profile)
thank you :D
Raymond Cheng (view profile)
Thanks for your sharing.
Sophie (view profile)
I found.. it's just a question of angle: should have done:
t =  ellipse_t.phi;
Sophie (view profile)
@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 (view profile)
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 = (xx0)*cos(t)+(yy0)*sin(t);
Y = (xx0)*sin(t)+(yy0)*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 (view profile)
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 helpme.... (sorry my english)
Love it. Works as advertised.
Exactly what I was after, perfect. Does what it says, does it well and easily.
Awesome! I just had to modify the program to plot the nonrotated ellipse instead, otherwise this was perfect!
THANKS
Incredible useful and practical mfile
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
Works well for me. Easy interface, and nice optional visualization.
Great. The plot function was a great help as well.
Execellent!!! Very well written and extremely easy to use.
Thanks!
Thank you. In presentday, i do my research about face detection in color image.
Does what it claims to. Quite helpful!
Thanks Ohad  it works well with color discrimination data. Best wishes, Peter
This method works excellent with all my noisy data. Also well written and extremely easy to use.