Classroom Tips and Techniques
Fitting Circles in Space to 3D Data
Robert J. Lopez
Emeritus Professor of Mathematics and Maple Fellow
© Maplesoft, a division of Waterloo Maple Inc., 2010

Introduction


An algebraic solution was given in [1] for the problem of fitting a circle to threedimensional data points not lying in the plane. In [2], this solution is crafted as a Matlab project whose steps we articulate in Table 1.
1.

Translate all points by their centroid so that Π, the plane containing the fitting circle, will pass through (or nearly through) the origin.

2.

Obtain the 3 × 3 factor U from the singularvalue decomposition of , the 3 × matrix whose columns are vectors from the origin to the translated points.

3.

The first two columns in are , an orthonormal basis for Π.

4.

The third column in is , a unit vector normal to Π.

5.

The columns in the 3 × matrix are the components of the vectors with respect to the basis vectors . Since n is orthogonal to Π, the entries in the third row of are small.

6.

A leastsquares fit to a circle is made in the plane Π by treating the first two entries in each column of as , a point in an plane.

Table 1 Essential steps in an algebraic solution to the problem of fitting a circle in space to 3D data

We will implement the algebraic solution in Maple, then provide an alternate analytic solution based on the direct minimization of the sum of squares of deviations. The analytic technique, suggested by [3], simultaneously fists a plane to the data, while finding the radius of the circle by minimizing the sum of squares of deviations of each point from a line parallel to the normal to the fitting plane.


Algebraic Solution


Each of the ten columns of the following matrix are the Cartesian coordinates of data points taken from [2].
To each data point, we add random noise of magnitude less than 1.5. The perturbations for each coordinate are in the matrix
so that the noisy data are contained in the matrix
The coordinates of the centroid of these points appear as the components of the following vector.
The coordinates of the centroid is subtracted from each point; the columns of the resulting matrix
are vectors ,, from the origin to each of the translated points.
The singular value decomposition of this matrix provides the factor and the vector S of singular values. The matrix contains orthonormal bases for the column space of and for its orthogonal complement.
The third component of
=
the vector of singular values, is a measure of how close the translated points are to the plane Π. As the magnitude of the noise in the data increases, this third singular value gets larger. Hence, the algebraic approach advocated in [1] is less robust than the analytic approach sketched in the next section.
If the translated points lie on a plane through the origin, the column space of will be twodimensional, that is, the plane Π. The orthogonal complement will be spanned by a vector orthogonal to the plane Π. Consequently, the first two columns of comprise the orthonormal basisfor the plane Π, and the third column of namely,
is a unit vector normal to Π.
That is an orthogonal matrix is verified by the calculation , implemented in Maple via
The columns of the product that is,
are the components of the along and n. Thus, the first two rows of are taken as  and coordinates of points in an plane spanned by and In this plane, a leastsquares fit of a circle to the pairs represented by these rows is effected by writing
in the form
where
(as suggested in [2]). This strategy leads to the column vectors
and to the matrix A and vector b
forming the system Ac = b.
The leastsquares solution for c is then
Consequently, in the translated and tilted plane, the fitting circle has
for the radius and center, respectively. Again, note that the center is in coordinates relative to the orthonormal basis . The coordinates of this center, in the space of the original data, are then
In the presence of noise, the third row in the matrix does not represent a zero component along the normal to the plane spanned by Thus, an essential condition for the validity of the algebraic solution is not met. The theory requires that the data lie in a plane for the algebraic method illustrated here to be applicable. Yet, as we will see, the algebraic solution is close to the analytical solution obtained below.


Analytic Solution


An analytic solution of the problem of fitting a circle to points in threedimensions begins with the following definitions.
Set equal to zero, the first expression defines the plane through the point and having
i + j + k
for its normal of length one. The second expression is the distance of the point to the line parallel to the normal, and passing through . Figure 1, using concepts from multivariate calculus, is a basis for obtaining this function.

When the point is projected onto the line whose direction is given by the unit vector , a right triangle is formed. Elementary trigonometry gives , so . Because has length one, we therefore can write
The desired result follows from

Figure 1 Distance from the point to the line with direction

The noisy data are stored in the matrix
=
To obtain an estimate of where the center of the sphere might lie, we form lists of , , and coordinates
then compute the average of the extreme values for each coordinate.
The sum of squares of the deviations of the data points from both the plane and the line is given by
and the minimum of this sum is given by
where we constrain the normal vector to have unit length. Thus, the center and radius of the sphere are
while the normal to the plane containing the fitting circle is
The equation of the plane containing the fitting circle is then
Figure 2 shows this plane, the data points, and the computed center.
Click the icon below to generate Figure 2.

Figure 2 The data points, the computed center, and the leastsquares plane containing the data points

Obtaining a graph of the fitting circle requires considerably more work.
In a spherical coordinate system centered at , and defined by
the equation of the plane containing the analyticallyfound circle becomes
from which we obtain as
To get a fourquadrant version of this result, write
which specializes to
The parametric representation of the circle in space is then
Figure 3 shows the fitted points, the fitting circle, the plane in which this circle lies, and the center of the fitting circle.
Click the icon below to generate Figure 3.

Figure 3 The leastsquares circle and the noisy data to which it has been fit.



References


[1]

Carl C. Cowen, "A Project on Circles in Space." In Resources for Teaching Linear Algebra. Edited by David Carlson, Charles R. Johnson, David Lay, A. Duane Porter, Ann Watkins, and William Watkins. Washington, D.C.: MAA, 1996.

[2]

ATLAST: Computer Exercises for Linear Algebra, Steven Leon, Eugene Herman, Richard Faulkenberry, PrenticeHall, Inc., Upper Saddle River, NJ, 1996.

[3]

Craig M. Shakarji, "LeastSquares Fitting Algorithms of the NIST Algorithm Testing System," Journal of Research of the National Institute of Standards and T2chnology, Volume 103, Number 6, NovemberDecember 1998


Legal Notice: © Maplesoft, a division of Waterloo Maple Inc. 2010. Maplesoft and Maple are trademarks of Waterloo Maple Inc. This application may contain errors and Maplesoft is not liable for any damages resulting from the use of this material. This application is intended for noncommercial, nonprofit use only. Contact Maplesoft for permission if you wish to use this application in forprofit activities.
