## Orthogonal least squares fitting of a sphere

Recently I presented a linear method to obtain centre and radius of a sphere given four or more points on its surface, not all beginning to the same plane.

This method is not completely satisfactory when working with more than four points if the points are affected by errors (that is always, if the points come from measurement). The least squares solution is the set of values $\alpha, \beta, \gamma, \delta$ which minimizes the (squared) 2-norm of the residual (algebraic fitting):

$\sum_{i=1}^m (\alpha x_i + \beta y_i + \gamma z_i - \delta)^2$

that is

$\sum_{i=1}^m ((x_i - x_c)^2 + (y_i - y_c)^2 + (z_i - z_c)^2 - r^2)^2$

This quantity is not very significant. Yes, it is related to the goodness of fit: it is zero if all the points lie on the sphere and increases as they depart, but its geometric meaning is awkward, as each term of the sum is proportional to the squared area of the annulus between a circumference centered in $\left[ x_c \hspace{1 mm} y_c \hspace{1 mm} z_c \right]$ and passing through $\left[ x_i \hspace{1 mm} y_i \hspace{1 mm} z_i \hspace{1 mm} \right]$, and a circumference with same centre and of radius $r$.

The squared residuals of the equations for other shapes will in general have a different geometric meaning, so that what one minimizes in this way depends on the shape.

Even more important, if the points are not evenly distributed on the sphere but are lumped in a small zone, this method will exhibit a significant bias, finding a centre too close to the lump and a radius too small. This because the area of the annulus decreases much more by moving the centre toward the lump and correspondingly decreasing the radius, than by moving the centre sideways to adjust the thickness of the annulus.

What we really should minimize is the sum of the unsigned distances from each point to the sphere; but this quantity has a discontinuous gradient in its minima, so that it is easier to minimize the sum of the squared distances from each point to the sphere (*):

$\sum_{i=1}^m (\sqrt{(x_i - x_c)^2 + (y_i - y_c)^2 + (z_i - z_c)^2} - r)^2$

This is said geometric fitting, or orthogonal fitting, or orthogonal regression. When fitting a function, the term total least squares is often found in literature, because the method takes into account not only errors on the ordinates but on the abscissas too.

Unfortunately, this problem is nonlinear, and harder to solve than algebraic fitting.

In my next post, I will describe a little variant to the method presented by Dave Eberly.

(*) Addendum 120920. I later learned that Giuseppe Calafiore of Politecnico di Torino devised a method for n-dimensional spherical and ellipsoidal fitting based on SDP (Semi-Definite Programming) which can easily deal with the $L_1$ norm too.