球到点的线性最小二乘拟合

Question

我正在寻找一种算法来找到点云和球体之间的最佳拟合。

也就是说，我想最小化

其中 C 是球体的中心，r 是球体的半径，每个 P 是我的 n 点集中的一个点。变量显然是 Cx、Cy、Cz 和 r。 就我而言，我可以事先获得已知的 r，只留下 C 的分量作为变量。

我真的不想使用任何类型的迭代最小化（例如牛顿法、Levenberg-Marquardt 等） - 我更喜欢一组线性方程或明确使用 SVD 的解决方案。

Answer 1

没有即将出现的矩阵方程。你选择E的行为很糟糕；它的偏导数甚至不是连续的，更不用说线性了。即使目标不同，这个优化问题似乎从根本上来说是非凸的；具有一个点 P 和非零半径 r，最优解集是围绕 P 的球体。

您可能应该在交流中寻求更多优化知识。

Answer 2

您可能会发现以下参考资料很有趣，但我要警告您 您需要熟悉几何代数 - 特别是共形几何代数来理解 数学。然而，该算法的实现非常简单 标准线性代数技术并且不是迭代的。

需要注意的是，至少所呈现的算法既适合中心又适合 半径，您也许可以找到一种方法来约束拟合，从而约束半径。

使用 n 维欧几里德空间中的 k 球体的总最小二乘拟合 (n+ 2)-D 等距表示。 L Dorst，《数学成像与视觉杂志》，2014 年第 1-21 页

您可以从以下位置获取副本 Leo Dorst 的研究门页面

最后一件事，我与作者没有任何联系。

Answer 3

制作矩阵方程的简短描述可以在这里找到。

我已经看到WildMagic Library使用迭代方法（至少在版本4中）

Answer 4

您可能对最佳拟合 d 维球体感兴趣，即最小化到中心的平方距离的总体方差；它有一个简单的解析解（矩阵微积分）：请参阅 Cerisier 等人的开放获取论文的附录。在《计算机杂志》中。生物。 24(11), 1134-1137 (2017), https://doi.org/10.1089/cmb.2017.0061 它在数据点加权时起作用（甚至对于连续分布也起作用；作为副产品，当 d=1 时，会检索到众所周知的不等式：峰度始终大于偏度平方加 1）。

Answer 5

如果不进行迭代，很难做到这一点。

我将按照以下步骤进行：

1. 通过平均所有点的 (X,Y,Z) 坐标来找到总体中点

2. 根据该结果，找到到中点的平均距离 Ravg，决定确定或继续..

3. 从集合中删除距离步骤 2 中找到的 Ravg 太远的点

4. 返回步骤 1（再次平均点，产生更好的中点）

当然，这需要（2）和（4）的一些条件，这取决于您的点云的质量！

Answer 6

Ian Coope 有一个有趣的算法，他使用变量的变化将问题线性化。拟合非常稳健，虽然它稍微重新定义了最优性条件，但我发现它通常在视觉上更好，尤其是对于噪声数据。

库普论文的预印本可在此处获取：https://ir.canterbury.ac.nz/bitstream/handle/10092/11104/coope_report_no69_1992.pdf.

我发现该算法非常有用，因此我在 scikit-guess 中将其实现为

skg.nsphere_fit

。假设您有一个

(m, n)

数组

p

，由 M 个维度为 N 的点组成（此处 N=3）：

r, c = skg.nsphere_fit(p)

半径

r

是标量，

c

是包含中心的

n

向量。

Answer 7

幸运的是，对于球体拟合问题确实有一个“直接”的解决方案。

展开并重新排列球体的笛卡尔方程：

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

你得到：

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

左侧部分包括球体参数（中心和半径），右侧部分是点的平方坐标之和。

设置：

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

允许将拟合问题写为线性最小二乘问题（很容易解决）A.x = b，其中x = [c_x，c_y，c_z，t]（与点的x不同！）。

在 python 中，假设你的点存储在形状 (n, 3) 的 numpy 数组中，你会有：

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

如果您有兴趣了解它是如何工作的，您可以查看这些教程。

干杯！