MATLAB Answers

How do I fit an ellipse to my data in MATLAB?

141 views (last 30 days)
I have some measured data stored in vectors "x" and "y". I want to find the ellipse that best fits this data (in the least squares sense).

Accepted Answer

MathWorks Support Team
MathWorks Support Team on 2 May 2012
This question can be viewed as both a matrix problem and as a nonlinear least squares question.
The ellipse equation is:
(x - c(1))^2/r(1)^2 + (y - c(2))^2/r(2)^2 = 1
LINEAR LEAST SQUARES
--------------------
First, view the problem as a linear least squares problem in terms of the parametric equations:
x = a(1) + a(2)*cos(t);
y = a(3) + a(4)*sin(t) ;
Here, you are trying to find "a" to determine the best fit of x and y (given t) to these equations in the least-squares sense. (Assume you do not know where the ellipse is centered. If it is centered at the origin, then a(1) and a(3) are zero and can be left out of the equations.)
Assume the data (x, y, and t) is generated by
t = (0:pi/10:2*pi)';
x = 5 + 4*cos(t) + rand(size(t));
y = 8 + 2*sin(t) + rand(size(t));
Then you want to perform the operation:
minimize || W*a - [x;y] ||_2
a
where W is:
i = size(t); % t is a column vector
W=[ones(i) cos(t) zeros(i) zeros(i) ;
zeros(i) zeros(i) ones(i) sin(t) ] ;
Determine "a" by solving the over determined system:
W*af = [x;y]
using backslash
af = W \ [x;y];
Measure the goodness of the fit by looking at the 2-norm of the residual vector:
norm(W*a - [x;y])
You can also visually check the results:
plot(x,y,'*'); hold on
xnew = af(1) + af(2)*cos(t);
ynew = af(3) + af(4)*sin(t);
plot(xnew, ynew ,'r')
NONLINEAR LEAST SQUARES
------------------------
If "t" is not known (or even if it is), you can also determine the best fit, in the least squares sense, of x and y to an ellipse. That is, you want to perform the operation:
minimize || f(a,x,y) ||_2
a
where f = (( x-a(1) ).^2)/a(2).^2 + (( y-a(3) ).^2)/a(4).^2 - 1 . In other words, you are minimizing the sum of the squares of the residuals:
min (sum ( f .* f ) )
a
You can use the nonlinear least squares function (LSQNONLIN) in the Optimization Toolbox to solve this problem. Using the x and y data from the previous problem, and the initial guess [10 10 10 10], we can determine the ellipse with:
a0 = [10 10 10 10];
f = @(a) ((x-a(1)).^2)/a(2).^2 + ((y-a(3)).^2)/a(4).^2 -1;
options = optimset('Display','iter');
af = lsqnonlin(f, a0, [], [], options);
Anonymous functions were introduced in MATLAB 7.0 (R14). If you are using a previous version, you will need to create an inline function using:
f = inline('((x-a(1)).^2)/a(2).^2 + ((y-a(3)).^2)/a(4).^2 -1', 'a','x','y');
and pass x and y as additional parameters to LSQNONLIN:
af = lsqnonlin(f, a0, [], [], options, x, y);
Since you have variables in the denominator and numerator, they can shrink and grow together and you could obtain very large or very small coefficients. To avoid this, fix a(1) and a(2) (providing the center of the ellipse) to be the average of the data x and y. Then, vary the other two parameters:
a0 = [10 10];
options = optimset('Display','iter');
c = [mean(x) mean(y)];
f = @(a) ((x-c(1)).^2)/a(1).^2 + ((y-c(2)).^2)/a(2).^2 -1;
af = lsqnonlin(f, a0, [], [], options);
plot(x,y,'*'), hold on
plot(c(1), c(2), 'r*')
t=0:pi/10:2*pi;
plot(c(1) + af(1)*cos(t), c(2) + af(2)*sin(t), 'r')
If you are using a version of MATLAB prior to MATLAB 7.0 (R14), you will need to create an inline function
f = inline('((x-c(1)).^2)/a(1).^2 + ((y-c(2)).^2)/a(2).^2 -1', 'a', 'x', 'y', 'c');
and change the LSQNONLIN command to:
af = lsqnonlin(f, a0, [], [], options, x, y, c);
If you do not have the Optimization Toolbox, you can also use the FMINSEARCH function to solve this problem. FMINSEARCH is not specialized like LSQNONLIN to handle the least squares problem, but the objective function can be modified to account for this. For FMINSEARCH, the sum of squares evaluation must be performed within the objective funtion.
f = @(a) norm(((x-c(1)).^2)/a(1).^2 + ((y-c(2)).^2)/a(2).^2 -1);
af = fminsearch(f, a0, options);
If you are using a version of MATLAB prior to MATLAB 7.0 (R14), you will need to create an inline function:
f = inline('norm(((x-c(1)).^2)/a(1).^2 + ((y-c(2)).^2)/a(2).^2 -1)', 'a', 'x', 'y', 'c');
af = fminsearch(f, a0, options, x, y, c)

  4 Comments

Show 1 older comment
Gert Kruger
Gert Kruger on 27 Oct 2016
Thanks for the answer, I would just like to point out that the center of the ellipse is not necessarily the "average of the data x and y" if it was not uniformly sampled.
Jessica Hiscocks
Jessica Hiscocks on 8 Mar 2017
This answer was extremely useful to me. However the use of rand(size(t)) will only generate positive values, meaning that the returned center will always be shifted from the input values. This will give the residual vector a higher, non-representative value. Also this line should probably be adjusted to read norm(W*af - [x;y])
Kumar Arumugam
Kumar Arumugam on 21 Apr 2019
can someone point me to the source of the linear least square program explained above?

Sign in to comment.

More Answers (2)

Shahryar Ahmad
Shahryar Ahmad on 14 Jan 2018
Edited: Shahryar Ahmad on 14 Jan 2018
You can also use Ellipse Direct Fit

  0 Comments

Sign in to comment.


SAI MANOJ KONDAPALLI
SAI MANOJ KONDAPALLI on 19 Nov 2018
Can someone help me with this error.
I am trying to run the code for Fitting a conic to a given set of points, using data from coordinates of an ellipse, detected using Canny edge detection algorithm
First it is asking me to use (.*)
timesDot.JPG
And then, when I change to (.*), It is showing the following error
timesStar.JPG

  1 Comment

Sign in to comment.

Sign in to answer this question.