Application Center - Maplesoft
  Submit your work
Application Center Applications Classroom Tips and Techniques: Fitting an Ellipse to Noisy Data Preview

App Preview:

Classroom Tips and Techniques: Fitting an Ellipse to Noisy Data

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application

 

[Inserted Image]

Fitting an Ellipse to Data

Maplesoft, a division of Waterloo Maple Inc., 2004

Introduction

This worksheet uses Maple to generate an ellipse which is fitted to a given set of data points.

The following Maple techniques are highlighted:

Use of the Optimization package to find an optimal ellipse

Plotting the ellipse

Translating the ellipse to standard form

Initialization

> restart;
interface( warnlevel = 0 ):
with( LinearAlgebra ):
with( Optimization ):
with( plots ):
interface( rtablesize=15 ):

Data

A slightly modified set of data points - needed to see if a new algorithm would work on points not lying on an ellipse.

 

> Q := [[10,1.5],[15.87785,4], [19.51057,3.934443], [19.51057,3], [15.87785,1.5], [10,-0.59], [4.122147,-1.8], [0.489435,-2], [0.489435,-1.00739], [4.122147,0.5]];

Q := [[10, 1.5], [15.87785, 4], [19.51057, 3.934443], [19.51057, 3], [15.87785, 1.5], [10, -.59], [4.122147, -1.8], [.489435, -2], [.489435, -1.00739], [4.122147, .5]]

Plot of Data Points

> p1 := plot(Q, style=point, symbol=circle, color=black):
p1;

[Plot]

Least Squares Computation

The general conic:

 

> F := a*x^2+b*x*y+c*y^2+d*x+e*y+f;

F := a*x^2+b*x*y+c*y^2+d*x+e*y+f

The 10 data points generate the following 10 equations

 

> for k from 1 to 10 do
        eq||k := eval(F, {x=Q[k][1], y=Q[k][2]});
od;

eq1 := 100*a+15.0*b+2.25*c+10*d+1.5*e+f

eq2 := 252.1061206*a+63.51140*b+16*c+15.87785*d+4*e+f

eq3 := 380.6623417*a+76.76322556*b+15.47984172*c+19.51057*d+3.934443*e+f

eq4 := 380.6623417*a+58.53171*b+9*c+19.51057*d+3*e+f

eq5 := 252.1061206*a+23.816775*b+2.25*c+15.87785*d+1.5*e+f

eq6 := 100*a-5.90*b+.3481*c+10*d-.59*e+f

eq7 := 16.99209589*a-7.4198646*b+3.24*c+4.122147*d-1.8*e+f

eq8 := .2395466192*a-.978870*b+4*c+.489435*d-2*e+f

eq9 := .2395466192*a-.4930519246*b+1.014834612*c+.489435*d-1.00739*e+f

eq10 := 16.99209589*a+2.0610735*b+.25*c+4.122147*d+.5*e+f

Sum of squares of residuals:

 

> E := add((eq||k)^2,k=1..10):

Minimize sum of squares of residuals subject to constraint 4*a*c-b^2 = 1 .  (See the paper Direct Least Square Fitting of Ellipses by Fitzgibbon, Pilu, and Fisher, pulled of the web by Google this morning.)

 

> V := Minimize(E, {4*a*c-b^2=1});

V := [44.7391881223692850, [e = 12.6215148366343328, c = 3.40240907787633073, a = .340856178828092715, b = -1.90759760651448107, f = 12.2687559064070886, d = -5.02190104894557177]]

The values of the parameters are given by

 

> W := V[2];

W := [e = 12.6215148366343328, c = 3.40240907787633073, a = .340856178828092715, b = -1.90759760651448107, f = 12.2687559064070886, d = -5.02190104894557177]

Substitute into F:

 

> sol := eval(F, W);

sol := .340856178828092715*x^2-1.90759760651448107*x*y+3.40240907787633073*y^2-5.02190104894557177*x+12.6215148366343328*y+12.2687559064070886

Modest rearrangement:

 

> Sol := sol/coeff(sol,x,2);

Sol := .9999999999*x^2-5.596488271*x*y+9.981949248*y^2-14.73319646*x+37.02885740*y+35.99393724

Graph of least-squares ellipse and data points:

 

> p2 := implicitplot(Sol, x=0..20, y=-2.5..4.5, grid=[80,80], color=red):
display([p1,p2]);

[Plot]

Translation and Rotation to Standard Form

Formulas for determining center, and rotation angle:

 

> H := (2*c*d-b*e)/(b^2-4*a*c);
K := (2*a*e-b*d)/(b^2-4*a*c);
Theta := arctan(b/(a-c))/2;
P := Matrix([[2*a,b,d],[b,2*c,e],[d,e,2*f]]);
Trans := a*u^2+b*u*v+c*v^2-delta/2/(b^2-4*a*c);

H := (2*c*d-b*e)/(b^2-4*a*c)

K := (2*a*e-b*d)/(b^2-4*a*c)

Theta := 1/2*arctan(b/(a-c))

P := Matrix([[2*a, b, d], [b, 2*c, e], [d, e, 2*f]])

Trans := a*u^2+b*u*v+c*v^2-1/2*delta/(b^2-4*a*c)

Evaluation of the center (h,k) and rotation angle:

 

> h := eval(H,W);
k := eval(K,W);
theta := eval(Theta,W);
delta := Determinant(eval(P,W));
trans := eval(Trans,W);

h := 10.09635195

k := .9755237880

theta := .2786093397

delta := -13.85278071

trans := .340856178828092715*u^2-1.90759760651448107*u*v+3.40240907787633073*v^2-6.926390369

Equation in standard form in unrotated system:

 

> SOL := fnormal(simplify(eval(trans,{u=cos(theta)*x-sin(theta)*y, v=sin(theta)*x+cos(theta)*y})));
XY := select(has,SOL,{x,y}):
C := SOL-XY:
final := (XY = -C)/(-C);

SOL := 0.6802272166e-1*x^2-6.926390369+3.675242535*y^2

final := 0.9820803914e-2*x^2+.5306144093*y^2 = .9999999998

Graph of unrotated ellipse in standard form:

 

> implicitplot(final,x=-12..12,y=-2..2, scaling=constrained, tickmarks=[5,3]);

[Plot]

Calculation of semimajor and semiminor axes:

 

> semimajor := 1/sqrt(coeff(lhs(final),x,2));
semiminor := 1/sqrt(coeff(lhs(final),y,2));

semimajor := 10.09082048

semiminor := 1.372810146

Graph of smallest inscribed circle and largest circumscribed circle:

 

> p3 := plot([h+semimajor*cos(t),k+semimajor*sin(t),t=0..2*Pi],color=black):
p4 := plot([h+semiminor*cos(t),k+semiminor*sin(t),t=0..2*Pi],color=green):
display([p||(1..4)],scaling=constrained);

[Plot]

 


Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material.. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.

[Inserted Image]