Fitting Circles and Ellipses to Data
Using the Least-Squares Method

Marc Renault

Abstract

       We address the problem of fitting circles and ellipses to given points in the plane.  There are two approaches we can use to find a circle or ellipse which best fits our data: algebraic fitting, and geometric fitting.  We note that even within these two approaches, there are different ways to solve the problem.  The best algebraic fit depends on which constraint we use, and the best geometric fit depends on which geometric quantity we decide to minimize.

       The main content of this paper is derived from [2] , but [3] was also extremely helpful (particularly with algebraic fitting), and [1] gives a nice alternative approach to the geometric fitting for circles.

1. General Concepts

       Let us first make clear what the difference is between an ellipse which best fits data algebraically, and one which best fits data geometrically.   An ellipse which is given by the algebraic equation F(x) = 0 is an ``algebraic fit'' if its parameters (e.g. the coefficients) are determined in a least squares sense.  An ellipse is a ``geometric fit'' if the sum of the squares of the distances to the given points is minimal.  The geometric fit, much more so than the algebraic fit, gives an ellipse which visually approximates the data points.

       As we will see, finding a best algebraic fit requires solving a linear least-squares problem, while finding a best geometric fit requires solving a nonlinear least-squares problem.  The latter problem is this:  Let u = (u₁,  u₂,¼,  u_n)^T be a vector of unknowns and let f(u) = 0 be a system of m nonlinear equations in those unknowns.  Then we wish to find the vector [(u)\tilde] which satisfies

m å i = 1  fi( ~ u   )2 = min u  m å i = 1  fi(u)2.

Another way of saying this is that we want the u which will minimize the square of the 2-norm of ( f₁(u),  f₂(u),¼,  f_m(u)) ^T.  Throughout this paper ||  ^.  || will denote the 2-norm.

       All the papers I researched use the ``Gauss-Newton'' algorithm for solving the nonlinear least-squares problem.  This is an iterative method which solves a series of linear least-squares problems, updating its best guess with every iteration.

2. Fitting a Circle

2.1 Finding an Algebraic Fit

       Let x,  b Î R² and let a ¹ 0.  Then the expression  

F(x) = axTx+bTx+c = 0

determines a circle in the plane.  Notice, however, that multiplying a,  b, and c by a constant determines the same circle.  (And if a = 0 then F(x) is a line.)

       Let m points in the plane, (x₁₁,x₁₂),  (x₂₁,x₂₂),¼,(x_m1,x_m2), be given.  Let u = ( a,b₁,b₂,c) ^T and let

B = é ê ê ê ë x11+x12 x11 x12 1 : : : : xm1+xm2 xm1 xm2 1 ù ú ú ú û .

Then consider the (linear) system of m equations in 4 unknowns

Bu = 0.

Ideally, if all m points actually lie on a circle, then this system has a solution.  If all m points do not lie on a circle then there is no solution.  Our only hope is to consider Bu = r and find a u which minimizes ||r||.  Furthermore, we need to put some constraints on u since, as noted previously, there are many different choices for u which determine the same circle.  It is common to require u₁ = a = 1, but GGS [2] use the constraint ||u|| = 1.  Then the matter of finding a least-squares solution to the linear system is routine and once a good u is found, the center and radius of the circle are easily determined.  See the dashed circle in figure 2.1.

Remark 1        Although the process of finding the best algebraic fit for a set of points is fairly straightforward, the geometric meaning of the resulting circle is unclear.  One often ends up with a circle which, upon visual inspection, does not appear closely related to the set of points we are approximating.  However, GGS [2] comment that ``the algebraic solution is often useful as a starting vector for methods minimizing the geometric distance.''  A starting vector being necessary, of course, for the Gauss-Newton method.

Remark 2        We will see that in the case of fitting an ellipse algebraically, the choice of constraint really makes a difference in the resulting figure.  This is a major topic in [3].  It is unclear whether choosing a different constraint in the case of the circle, for example choosing a = 1, will make a difference in the resulting circle.  This is not addressed in any of the papers.

2.2 Finding a Geometric Fit

       Consider the m data points x₁,...,x_m in R².  Given a circle with center z and radius r, let

di(z1,z2,r) = |   ||z-xi||-r  |

the distance from the data point x_i to the circle.  Letting u = (z₁,z₂,r), we wish to find [(u)\tilde] satisfying the expression

m å i = 1  di( ~ u   )2 = min u  m å i = 1  di(u)2

thus determining the circle which minimizes the sum of the squares of the distances from the data points to the circle.  Now we have a nonlinear least squares problem which can be solved using the Gauss-Newton method. The solid circle in figure 2.1 below is a geometric fit to the data points.

Remark 1        Finding the geometric fit seems much more intuitive.  The disadvantage to this approach is that we must solve a nonlinear least-squares problem.  Of course the method for doing this is well known, and an initial guess can be easily found by first finding an algebraic fit.

Remark 2        In [2], GGS also develop the method for finding a geometric fit based on the parametric form of the circle which, they claim, is commonly used.

3. Fitting an Ellipse

3.1 Finding an Algebraic Fit

       This is where things really start to get interesting.  Any conic section in the plane can be expressed by the equation

xTAx+bTx+c = 0.

where A is a 2 by 2 real symmetric matrix.  When A is positive definite and c < 0 we are assured that the above equation gives us an ellipse.  Expanding the above equation gives us

a11x12+2a12x1x2+a22x22+b1x1+b2x2+c = 0.

Let our m data points be (x₁₁,x₁₂),  (x₂₁,x₂₂),¼,(x_m1,x_m2) as before, and let

B = é ê ê ê ë x112 x11x12 x122 x11 x12 1 : : : : : : xm12 xm1xm2 xm22 xm1 xm2 1 ù ú ú ú û .

Let u = (a₁₁,  2a₁₂,  a₂₂,  b₁,  b₂,  c)^T.  Then as before, we would like to find [(u)\tilde] such that

||B ~ u   || = min u  ||Bu||.

Furthermore, we must again put some constraints on our choice of u. Varah [3] considers [3] possible constraints:

1. SVD Constraint:

   Minimize ||Bu|| subject to the constraint that ||u|| = 1.  We can see clearly that the solution to this problem is the vector (of norm 1) which corresponds to the smallest singular value of B.

2. Linear Constraint:

   Minimize ||Bu|| subject to the constraint that v^Tu = 1 where v is some fixed vector chosen in advance.  In other words, we require that some linear combination of the elements of u always equal 1.  For example, suppose as in the circle case that u = ( a,b₁,b₂,c) ^T.  Then, setting v = (1,0,0,0)^T and requiring v^Tu = 1 is equivalent to requiring that a = 1.  Varah proves that the solution to this problem is the unique vector

   u = (BTB)-1v vT(BTB)-1v .

3. Quadratic Constraint:

   Minimize ||Bu|| subject to the constraint that ||Su|| = 1 for some general matrix S where the multiplication Su makes sense.  Varah gives a solution for u under this constraint, which is beyond the scope of this paper.

       Hence, no matter what our constraint is, we can find a u which determines an ellipse.  However, as Varah remarks, ``when the data are not well-fitted, the solutions obtained using the three constraints can be very different indeed.''  See figure 3.1 below. In particular, note the three ellipses approximating the lower set of data points.

       Furthermore, suppose we have a set of data points, P, and an ellipse E which fits our data algebraically subject to some constraint.  Let T be any Euclidean transformation (a combination of rotation and translation).  Ideally, we would like the ellipse which algebraically fits T(P) to be T(E).  However, it turns out that constraint (a) fails with respect to this criterion.  In figure 3.1, we see that the ellipse which is generated using the SVD constraint changes considerably when the data points are translated.  However, the ellipses generated using the other constraints are translated along with the data points.  In fact, constraints (b) and (c) are invariant under all Euclidean transformations.

       Finally we remark that although the algebraic fit has no immediately apparent relation to the geometric fit, the algebraic fit for some constraints appears to visually approximate the data points geometrically.  Compare the two graphs in figure 3.2 below.

3.2 Finding a Geometric Fit

       GGS approaches this problem strictly from the parametric point of view.  We will see what the nonlinear least-squares problem looks like in this case.  The set of points x satisfying

x = é ê ë acosf bsinf ù ú û

for some 0 £ f < 2p is an ellipse centered at the origin with an axis of length a running along the x-axis and an axis of length b running along the y-axis.  Thus an ellipse centered at z with an axis of length a at an angle of a to the x-axis is given by the set of points x satisfying

x = z+Q(a) é ê ë acosf bsinf ù ú û ,  where   Q(a) = é ê ë cosa -sina sina cosa ù ú û

for some 0 £ f < 2p.

       Let our usual m data points be given.  Let u = (z₁,z₂,a,b,a).  Then the distance between point x_i and the ellipse is given by

gi(u) = ê ê ê   ||xi-z||-Q(a) é ê ë acosfi bsinfi ù ú û    ê ê ê

where f_i is the angle between the x-axis and the line joining x_i and z.  So, finally, our nonlinear least-squares problem is to find the vector [(u)\tilde] which satisfies

m å i = 1  gi( ~ u   )2 = min u  m å i = 1  gi(u)2.

4. An Alternative Characterization of the Geometric Fitting Problem

       We conclude this paper with an observation from [1] concerning the geometric fit for circles.  Previously we sought

min u  m å i = 1  di(u)2

where d_i(z₁,z₂,r) = |   ||z-x_i||-r  | .  This makes sense and clearly we are minimizing the sum of the squares of the distances from the data points to the circle.  Coope [1] calls this the ``Total Least-Squares'' problem (TLS).  A disadvantage to this approach is that apparently outliers are given undue weight.  Consider figure 4.1 where an outlier noticeably shifts the circle.

       Coope proposes using the following ``Linear Least-Squares'' (LLS) problem instead: find

min u  m å i = 1  di(u)2

(1)

where d_i(u) = ||z-x_i||²-r².  At first this appears to be another nonlinear least-squares problem, but through a clever change of variables, Coope shows that this is actually a linear problem.  First, note that

di(u) = zTz-2zTxi+xiTxi-r2.

 Let

y: = é ê ê ê ë 2z1 2z2 r2-zTz ù ú ú ú û ,  and   bi: = é ê ê ê ë xi1 xi2 1 ù ú ú ú û .

We see that b_i^Ty = 2z^Tx_i-z^Tz+r², and so expression (1) becomes

min y  m å i = 1  (biTy-xiTxi)2

which we can write more compactly as

min y  | | By-w| |

where B = [ b₁  |  b₂  |  ¼   |  b_m] and w is the vector with components w_i = || x_i| | ².  Hence, (1) is indeed a linear least-squares problem, and as such is easily solved.  Changing variables back again allows us to easily recover the center and radius of the circle.

       The advantage of this approach is that the circle which is created is less tolerant of outliers.  Consider figure 4.2 which uses the same data as the previous figure.

Remark        It is interesting to note that Coope's paper [1] preceeds the other two, yet [2] and [3] do not mention this approach at all and do not use [1] as a reference.

       Additionally, it would be nice to see this approach used on a wider selection of data.  If the data points look like those in [2] and [3], what will the LLS circle look like?  How will it compare to the TLS circles?

References

1. I. D. Coope, Circle Fitting by Linear and Nonlinear Least Squares, J Optim Theory Appl, 76 (1993) pp. 381-388.

2. W. Gander, G. H. Golub, and R. Strebel, Least-Squares Fitting of Circles and Ellipses, BIT, 34 (1994) pp. 558-578.

3. J. M. Varah, Least Squares Data Fitting with Implicit Functions, BIT, 36 (1996), pp. 842-854.

---

File translated from T_EX by T_TH, version 2.25.
On 09 Feb 2000, 17:55.