Linear Least Squares Fit of Sphere to Points

Question

I'm looking for an algorithm to find the best fit between a cloud of points and a sphere.

That is, I want to minimise

where C is the centre of the sphere, r its radius, and each P a point in my set of n points. The variables are obviously Cx, Cy, Cz, and r. In my case, I can obtain a known r beforehand, leaving only the components of C as variables.

I really don't want to have to use any kind of iterative minimisation (e.g. Newton's method, Levenberg-Marquardt, etc) - I'd prefer a set of linear equations or a solution explicitly using SVD.

Granule

Answer 1

There are no matrix equations forthcoming. Your choice of E is badly behaved; its partial derivatives are not even continuous, let alone linear. Even with a different objective, this optimization problem seems fundamentally non-convex; with one point P and a nonzero radius r, the set of optimal solutions is the sphere about P.

You should probably reask on an exchange with more optimization knowledge.

Parabolize

You probably want to use something like sum[i=0..n]( |P_i - C|^2 - r^2 )^2 instead, so your derivatives will behave properly. And, because your problem will be nonlinear in any case, you are probably stuck with some form of iteration. – Extensor

Answer 2

You might find the following reference interesting but I would warn you that you will need to have some familiarity with geometric algebra - particularly conformal geometric algebra to understand the mathematics. However, the algorithm is straight forward to implement with standard linear algebra techniques and is not iterative.

One caveat, the algorithm, at least as presented fits both center and radius, you may be able to work out a way to constrain the fit so the radius is constrained.

Total Least Squares Fitting of k-Spheres in n-D Euclidean Space Using an (n+ 2)-D Isometric Representation. L Dorst, Journal of Mathematical Imaging and Vision, 2014 p1-21

Your can pull in a copy from Leo Dorst's researchgate page

One last thing, I have no connection to the author.

Oden

Answer 3

Short description of making matrix equation could be found here.

I've seen that WildMagic Library uses iterative method (at least in version 4)

Waistline

Answer 4

You may be interested by the best fit d-dimensional sphere, i.e. minimizing the variance of the population of the squared distances to the center; it has a simple analytical solution (matrix calculus): see the appendix of the open access paper of Cerisier et al. in J. Comput. Biol. 24(11), 1134-1137 (2017), https://doi.org/10.1089/cmb.2017.0061 It works when the data points are weighted (it works even for continuous distributions; as a by-product, when d=1, a well-known inequality is retrieved: the kurtosis is always greater than the squared skewness plus 1).

Giverin

Answer 5

Difficult to do this without iteration.

I would proceed as follows:

1. find the overall midpoint, by averaging (X,Y,Z) coords for all points

2. with that result, find the average distance Ravg to the midpoint, decide ok or proceed..

3. remove points from your set with a distance too far from Ravg found in step 2

4. go back to step 1 (average points again, yielding a better midpoint)

Of course, this will require some conditions for (2) and (4) that depends on the quality of your points cloud !

Alpestrine

Answer 6

Ian Coope has an interesting algorithm in which he linearized the problem using a change of variable. The fit is quite robust, and although it slightly redefines the condition of optimality, I've found it to be generally visually better, especially for noisy data.

A preprint of Coope's paper is available here: https://ir.canterbury.ac.nz/bitstream/handle/10092/11104/coope_report_no69_1992.pdf.

I found the algorithm to be very useful, so I implemented it in scikit-guess as skg.nsphere_fit. Let's say you have an (m, n) array p, consisting of M points of dimension N (here N=3):

r, c = skg.nsphere_fit(p)

The radius, r, is a scalar and c is be an n-vector containing the center.

Outlast

Answer 7

Fortunately for you there is indeed a "direct" solution to the sphere fitting problem.

Expanding and rearranging the Cartesian equation of the sphere:

(x - c_x)² + (y - c_x)² + (z - c_z)² = r²

you get:

2.x.c_x + 2.y.c_y + 2.z.c_z + r² - c_x² - c_y² - c_z² = x² + y² + z²

The left part includes the sphere parameters (center and radius) and the right part is the sum of the squared coordinates of the points.

Setting:

t = r² - c_x² - c_y² - c_z²

allows to write the fitting problem as a linear least squares problem (which is quite easily solved) A.x = b with x = [c_x, c_y, c_z, t] (not the same x as for the points!).

In python, assuming that your points are stored in a numpy array of shape (n, 3), you'd have:

import numpy as np

def fit_sphere(points):

    # Coefficient matrix and values
    A = np.column_stack((2*points, np.ones(len(points))))
    b = (points**2).sum(axis=1)
    # Solve A x = b
    x, res, _, _ = np.linalg.lstsq(A, b, rcond=None)
    # Sphere parameters
    center = x[:3]
    radius = np.sqrt(x[0]**2 + x[1]**2 + x[2]**2 + x[3])
    
    return (center, radius), res

If you're interested to see how it works you can check these tutorials.

Cheers !

Stenotype