CODEWITHSUNDEEP
pythonpoint-cloudsleast-squaresmodel-fitting
Kowal
Kowal

Reputation: 21

Robust Fitting of a Sphere to Noisy 3D Points

I must write a code that finds the parameters of a sphere represented by a set of noisy 3D points.

I tried to estimate the sphere parameters using the linear least squares technique (see Linear Least Squares Fit of Sphere to Points) but the result I got is not satisfactory.

Are there other approaches that are more robust to noise?

import numpy as np


def fit_sphere(points):
    """Sphere fitting using Linear Least Squares"""

    A = np.column_stack((2*points, np.ones(len(points))))
    b = (points**2).sum(axis=1)
    x, res, _, _ = np.linalg.lstsq(A, b)
    center = x[:3]
    radius = np.sqrt(x[0]**2 + x[1]**2 + x[2]**2 + x[3])
   
    return (center, radius), res

Upvotes: 2

Views: 2505

Answers (1)

Tinu
Tinu

Reputation: 2523

There is a pretty simple analytic solution to your problem. Recall the sphere equation, e.g. in 3D:

(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2,

where (x_0, y_0, z_0) is the center of the sphere. Given a set of points we can derive the optimal r^2 for the sum of the Euclidean distances for each of your points denoted by i.

L = \sum_i^n (\bar{x}_i^2 + \bar{x}_i^2 + \bar{x}_i^2 + r^2)^2,

where the bar indicates we have centered coordinates, we basically subtracted the mean. The next thing to do is to set the derivative of \frac{dL}{dr^2} = 0 and solve for r^2. You can do the derivations yourself. The final result is the following:

r_{opt}^2 = \frac{1}{n} \sum_i^n \bar{x}_i^2 + \bar{x}_i^2 + \bar{x}_i^2.

Here is an example in 2D:

import numpy as np
import matplotlib.pyplot as plt

# simple 2D example with 3 points
X = np.array([[0,1],[1.5,1.5],[2,0]])
center = X.mean(0)
X_centered = X - center
r_opt = np.sum(X_centered**2)/len(X_centered)

fig, ax = plt.subplots(figsize=(4,4))
ax.scatter(X[:,0],X[:,1])
ax.scatter(center[0],center[1],color='r')
shpere = plt.Circle((center[0],center[1]),np.sqrt(r_opt),fill=False,color='r')
ax.add_artist(shpere)
ax.set_xlim((center[0]-r_opt*1.1,center[0]+r_opt*1.1))
ax.set_ylim((center[1]-r_opt*1.1,center[1]+r_opt*1.1))

enter image description here

Red dot is the center, blue dots are your noisy measurements and the red circle is the fitted sphere. The example is 2D for visualization purposes, but the results generalizes to three or more dimensions. Sorry for not being able to use any MathJax, maybe this post would fit better in math.stackexchange.com.

Upvotes: -2

Related Questions