Fitting a Circle to a Number of Points

by Thomas Schramm, RZ TUHH, Germany

NOTE: This worksheet demonstrates how Maple can be used to fit a circle to a set of Datapoints

Introduction

At first glance some fitting problems may appear to be non-linear, whereas with a bit of thought they can be solved using linear fitting techniques.

For example the naive approach to the problem of fitting data to a circle would be to fit the parameter (x0, y0, r) of an equation of the following type to the data:

> eq := (x-x0)^2 + (y-y0)^2 = r^2;

Unfortunately the parameters occur in a quadratic form:

> expand(eq);

But all is not lost, if a set of datapoints are known to lie on a circle, the average quadratic distance of the points from a circle should be minimized. This problem can be written so that it is linear in the parameters to be determined. The general equation for a circle is:

where the center of the circle is [a, b], and the radius is given by .

The coordinates of the datapoints are inserted in the equations to give a system of linear equations. This over-determined system can be directly solved using the "leastsqrs" function in Maple's "linalg" package.

> restart:

The general equation for a circle is:

> Kreis := (x,y)->x^2+y^2-2*a*x-2*b*y-c=0;

A set of [x,y] coordinates for which we might want to fit a circle are:

> Px := [6, 3, 9, 2, 5]:   Py := [2, 2, 5, 3, 2]:
    Punkte := zip((x,y)->[x,y], Px, Py);

Punkte := [[6, 2], [3, 2], [9, 5], [2, 3], [5, 2]]

Using this data we construct the corresponding system of equations:

> gls := { op( zip( Kreis, Px, Py) ) };

gls := { 40 - 12 a - 4 b - c = 0, 13 - 6 a - 4 b - c = 0, 106 - 18 a - 10 b - c = 0, 13 - 4 a - 6 b - c = 0, 29 - 10 a - 4 b - c = 0 }

and solve this over-determined system using:

> lsg := linalg[leastsqrs](gls, {a,b,c});

The radius of the circle is given by:

> r := evalf( subs( lsg, sqrt(c + a^2 + b^2) ) );

r := 4.321897536

Substituting the values for a, b and c into our circle equation, gives the parameteric description of our circle-fit as:

> KF := subs(lsg, Kreis(x,y));

The solution curve (circle) of the implicit equation can be directly plotted using:

> g1 := plots[implicitplot]( KF, x=0..10, y=-11..11, thickness=2):

And the datapoints using:

> g2 := plot( Punkte, style=point, symbol=cross, color=blue):

We display both the data and the circle together.

> plots[display]( {g1,g2}, scaling=constrained );

Disclaimer: While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.