**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]];** |

**Plot of Data Points**

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

**Least Squares Computation**

The general conic:

**> ** |
**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; |

Sum of squares of residuals:

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

Minimize sum of squares of residuals subject to constraint . (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});** |

The values of the parameters are given by

Substitute into F:

Modest rearrangement:

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

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]); |

**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); |

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

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

Graph of unrotated ellipse in standard form:

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

Calculation of semimajor and semiminor axes:

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

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

* *

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.