球体与嘈杂 3D 点的稳健拟合[已关闭]

Question

我必须编写一段代码来查找由一组嘈杂的 3D 点表示的球体的参数。

我尝试使用线性最小二乘技术来估计球体参数（请参阅球体到点的线性最小二乘拟合），但我得到的结果并不令人满意。

还有其他对噪声更稳健的方法吗？

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 

Answer 1

您的问题有一个非常简单的分析解决方案。回想一下球体方程，例如3D 模式：

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

,

其中

(x_0, y_0, z_0)

是球体的中心。给定一组点，我们可以导出由

r^2

表示的每个点的欧几里得距离总和的最佳

i

。

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

,

当条形表示我们有中心坐标时，我们基本上减去了平均值。接下来要做的是设置

\frac{dL}{dr^2} = 0

的导数并求解

r^2

。您可以自己进行推导。最终结果如下：

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

.

这是一个 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))

红点是中心，蓝点是噪声测量结果，红圈是拟合球体。出于可视化目的，该示例是二维的，但结果可推广到三个或更多维度。抱歉，无法使用任何 MathJax，也许这篇文章更适合 math.stackexchange.com。