MATLAB: How to fit an ellipse to the data in MATLAB

circlecurvedataellipsefitfittingleastOptimization Toolboxregressionsquares

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

Best Answer

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)
Related Question