algorithm - 最小封闭圆柱体

Question

是否有一种算法可以为 3D 点云找到半径最小的封闭圆柱体？我知 Prop 有最小外接圆的 2D 情况已解决(例如此线程 Smallest enclosing circle in Python, error in the code )，但是否有任何适用于 3D 的工作方法？

---

编辑 1:OBB。以下是弧形点云的示例。这个工具找到了最小的外接圆https://www.nayuki.io/page/smallest-enclosing-circle

圆是由三个点定义的，其中两个点几乎位于一条直径上，因此很容易估计中心轴在哪里。点的“装箱”将产生明显偏离真实中心的方框中心。

我的结论是，OBB 方法并不通用。

---

EDIT2:PCA。下面是紧点云 vs 的 PCA 分析示例。具有异常值的点云。对于紧密的点云，PCA 可以令人满意地预测圆柱方向。但是，如果与主云相比有少量异常值，那么 PCA 基本上会忽略它们，从而产生与封闭圆柱体的真实轴相距很远的向量。在下面的示例中，封闭圆柱体的真实几何轴以黑色显示。

我得出结论，PCA 方法并不通用。

---

EDIT3:OBB 与 PCA 和 OLS。一个主要区别 - OBB 仅依赖于几何形状，而 PCA 和 OLS 依赖于点的总数，包括那些不影响形状的中间点。为了使它们更有效率，可以包括数据准备步骤。首先，找到凸包。第二，排除所有内部点。然后，船体上的点可能分布不均匀。我建议删除所有这些，只留下多边形船体，并用网格覆盖它，其中节点将是新点。将 PCA 或 OLS 应用于这种新的点云应该可以提供更准确的圆柱轴估计。

如果 OBB 提供了一个尽可能平行于封闭圆柱体轴的轴，那么所有这些都是不必要的。

---

EDIT4:已发布的方法。@meowgoesthedog:Michel Petitjean 的论文(“关于最小封闭圆柱体问题的代数解”)可能有所帮助，但我没有足够的资格将其转换为工作程序。作者自己做了(这里的模块 CYL http://petitjeanmichel.free.fr/itoweb.petitjean.freeware.html)。但在论文的结论中，他说:“和目前的软件，名为 CYL，可在 http://petitjeanmichel.free.fr/itoweb.petitjean.freeware.html 免费下载，既没有声称提供了这些方法的最佳实现，也没有声称比其他方法更好地工作圆柱体计算软件。”论文中的其他短语也给人留下了印象，这是一种实验方法，尚未得到彻底验证。无论如何，我会尝试使用它。

@Ripi2:Timothy M. Chan 的这篇论文对我来说也有点太复杂了。我不是那种水平的数学专家，能够转换为工具。

@Helium_1s2:这可能是一个很好的建议，但是，与上面的两篇论文相比，它的细节要少得多。此外，未经验证。

---

EDIT5:回复 user1717828。两个最远的点与圆柱轴。一个反例 - 立方体形状中的 8 个点，适合一个圆柱体。两点之间的最大距离 - 绿色对角线。明显不平行于圆柱轴。

Ripi2 的“中间点”方法:它仅适用于 2D。在 3D 情况下，圆柱轴可能不会与任何两点之间的单个线段相交。

Answer 1

概念性答案

1. 找到距离h 最大的两个点。它们位于圆柱体的表面上，连接它们的线将平行于圆柱体的轴线。

2. 将所有点投影到垂直于该轴的平面上。

3. 在该平面上找到它们之间距离最大 d 的两个点。它们定义了直径 d 为圆柱体直径的圆。

4. 包含所有点的具有最小体积*的圆柱体有

.

* 这假设只有一对点，它们之间的距离最大，定义了圆柱体的轴。如果有可能两对点共享最大值，则对每对点重复步骤 2-4 并选择直径最小的圆柱体。

Python 答案

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook

from numpy.linalg import norm
from scipy.spatial.distance import pdist, squareform

如果您还没有积分，请生成积分:

np.random.seed(0)
N = 30
M = np.random.randint(-3,3,(N,3))
print(M)

[[ 1  2 -3]
 [ 0  0  0]
 [-2  0  2]
 [-1  1 -3]
 [-3  1 -1]
 ...
 [ 1 -3  1]
 [ 0 -1  2]]

计算每对可能的点之间的距离，并选择距离最大的一对。

max_dist_pair = list(pd.DataFrame(squareform(pdist(M))).stack().idxmax())
p1 = M[max_dist_pair[0]]
p2 = M[max_dist_pair[1]]
print(f"Points defining cylinder faces: {p1}, {p2}")
print(f"Length of cylinder: {norm(p1-p2)}")

Points defining cylinder faces: [-1 -3 -3], [1 2 2]
Length of cylinder: 7.3484692283495345

用蓝色显示所有点，用红色显示最大分离的点。

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(*M.T, c = ['red' if i in max_dist_pair else 'blue' for i in range(N)])

ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")

plt.show()

这是旋转的同一个图，所以我们沿着两个红点之间的轴看。

上面的 View 与投影在垂直于圆柱体轴的平面上的点是一回事。找到包含此平面中的点的最小圆。为此，我们找到每个点相对于轴的位移，然后找到两点之间的最大距离。

perp_disp = (np.cross(p2-p1, M-p1))/norm(p2-p1) # Get perpendicular displacement vectors.
print(perp_disp)

[[-3.40206909  1.36082763  0.        ]
 [ 0.         -0.13608276  0.13608276]
 [ 1.36082763 -2.04124145  1.4969104 ]
 [-2.72165527  0.          1.08866211]
 [-1.36082763 -1.90515869  2.44948974]
 [ 0.68041382 -0.95257934  0.68041382]
 [ 2.72165527  0.68041382 -1.76907593]
 ...
 [ 0.          0.27216553 -0.27216553]
 [ 0.         -0.40824829  0.40824829]
 [ 2.72165527  0.27216553 -1.36082763]
 [ 2.04124145 -0.68041382 -0.13608276]]

最大距离是通过执行上面使用的相同 pdist 技巧找到的。

max_perp_disp_pair = list(pd.DataFrame(squareform(pdist(perp_disp))).stack().idxmax())
perp_p1 = M[max_perp_disp_pair[0]]
perp_p2 = M[max_perp_disp_pair[1]]
print(perp_p1, perp_p2)

[ 1  2 -3] [-3 -2  1]

最后，我们得到了圆柱体的直径。

print(norm(perp_p1 - perp_p2))
6.92820323028

可以包含这些点的圆柱体的最小体积是

注意事项

• 使用 Numpy 的成对距离函数 pdist 找到点之间的最大距离。然后使用 squareform 对其进行格式化，将其放入 Pandas DataFrame 中，这样就可以使用 idxmax 轻松找到两个点的索引。在没有 Pandas 的情况下，可能有更好的方法来做到这一点。

• 如果 np.cross 部分让您摸不着头脑，那么这只是找出点和线之间的最小距离。如果您有兴趣，我可以跟进更多细节，但是如果您绘制两条线的叉积，您会得到一个平行四边形，其中两条非平行边由线给出。这个平行四边形的面积与矩形的面积相同，矩形的长度等于其中一条线，宽度等于点到线的距离。