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

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.

Load packages and define initial data by clicking on the icon below.

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 two-dimensional, 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 least-squares 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 least-squares 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.

An analytic solution of the problem of fitting a circle to points in three-dimensions 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.

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.

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 analytically-found circle becomes

from which we obtain  as

To get a four-quadrant 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.