将曲线样条拟合到 3D 点云

Question

目标

我有一个 3D 刻面模型（例如 .off 文件），它可以看起来像一个管道/管子（参见示例图片）。目标是导出一个近似样条曲线（线和样条曲线的最佳组合），它使用 python 表示该管的 3D 骨架。

最先进的

同一领域的 Stackoverflow 帖子：

• how-to-fit-a-line-through-a-3d-pointcloud

一般：

• Fitting-Spline-Curves-through-Set-of-Unorganized Point Cloud
• Skeletons from point cloud
• Powercrust / NN-Crust？我读到了这些，但是我找不到python实现，也许我太笨了，无法自己实现。据我了解，基础是我已经知道的 delaunay/voronoi。但是我不知道任何进一步的步骤

我的方法（到目前为止）

从示例 facet 模型（图 1）开始，我使用 python 包将 3d 模型转换为点云（图 2）。该点云可用于体素化表示（图 3）。因此，这三种类型的数据是我的出发点。

基本上，这个问题对我来说似乎并不太复杂，但是我缺少一个起始逻辑。大多数研究论文对于各种更深远的任务来说过于复杂了。 一个想法是做一个 PCA 来导出组件的主轴，并沿着这些轴扫描。但是，这似乎不会以高效的方式产生良好的结果。 另一个想法是使用体素化网格并检测由于体素邻接导致的路径。 另一个想法是使用 KD-Tree 来评估最近点以检测正确的平面，以便通过其平面法线定义样条线方向。

我尝试的一种方法是从点云中选择 N 个随机点并搜索半径内的所有邻居 (cKDTree.query_ball_point)。我计算了所有相邻点的中心。这导致了图 4 中的结果。第一种方法的结果似乎很好，但它或多或少是对半径参数的调整。

图片 1：

图片 2：

图 3：

图 4：

Answer 1

Delaunay/Voronoi 方法可用于解决此问题，因为中轴是 Voronoi 图的子图 （例如，参见 Attali、Boissonnat 和 Edelsbrunner 的 this paper）。

在下面，我将演示从小半径 10 和大半径 100 的四分之一圆环表面采样的点示例的方法（中间路径/骨架从点 (100, 0, 0) 开始并在 ( 0, 100, 0))。

Voronoi 图是 3D Delaunay 四面体化的对偶图（从现在起我将使用三角剖分这个术语）。 可以使用 scipy 的 scipy.spatial.Delaunay 包计算 Delaunay 三角剖分。

下图是样本点（本例中为 200 个）及其完整的 Delaunay 三角剖分图 （三角剖分是用函数from here 绘制的）。

对应于 Delaunay 四面体的 Voronoi 顶点是四面体的外接球体的中心。 以下是计算这些 Delaunay 中心的函数，它是我之前的回答 here 对 2D 函数的扩展。

def compute_delaunay_tetra_circumcenters(dt):
"""
Compute the centers of the circumscribing circle of each tetrahedron in the Delaunay triangulation.
:param dt: the Delaunay triangulation
:return: array of xyz points
"""
simp_pts = dt.points[dt.simplices]
# (n, 4, 3) array of tetrahedra points where simp_pts[i, j, :] holds the j'th 3D point (of four) of the i'th tetrahedron
assert simp_pts.shape[1] == 4 and simp_pts.shape[2] == 3

# finding the circumcenter (x, y, z) of a simplex defined by four points:
# (x-x0)**2 + (y-y0)**2 + (z-z0)**2 = (x-x1)**2 + (y-y1)**2 + (z-z1)**2
# (x-x0)**2 + (y-y0)**2 + (z-z0)**2 = (x-x2)**2 + (y-y2)**2 + (z-z2)**2
# (x-x0)**2 + (y-y0)**2 + (z-z0)**2 = (x-x3)**2 + (y-y3)**2 + (z-z3)**2
# becomes three linear equations (squares are canceled):
# 2(x1-x0)*x + 2(y1-y0)*y + 2(z1-z0)*y = (x1**2 + y1**2 + z1**2) - (x0**2 + y0**2 + z0**2)
# 2(x2-x0)*x + 2(y2-y0)*y + 2(z2-z0)*y = (x2**2 + y2**2 + z2**2) - (x0**2 + y0**2 + z0**2)
# 2(x3-x0)*x + 2(y3-y0)*y + 2(z3-z0)*y = (x3**2 + y3**2 + z3**2) - (x0**2 + y0**2 + z0**2)

# building the 3x3 matrix of the linear equations
a = 2 * (simp_pts[:, 1, 0] - simp_pts[:, 0, 0])
b = 2 * (simp_pts[:, 1, 1] - simp_pts[:, 0, 1])
c = 2 * (simp_pts[:, 1, 2] - simp_pts[:, 0, 2])
d = 2 * (simp_pts[:, 2, 0] - simp_pts[:, 0, 0])
e = 2 * (simp_pts[:, 2, 1] - simp_pts[:, 0, 1])
f = 2 * (simp_pts[:, 2, 2] - simp_pts[:, 0, 2])
g = 2 * (simp_pts[:, 3, 0] - simp_pts[:, 0, 0])
h = 2 * (simp_pts[:, 3, 1] - simp_pts[:, 0, 1])
i = 2 * (simp_pts[:, 3, 2] - simp_pts[:, 0, 2])

v1 = (simp_pts[:, 1, 0] ** 2 + simp_pts[:, 1, 1] ** 2 + simp_pts[:, 1, 2] ** 2) - (simp_pts[:, 0, 0] ** 2 + simp_pts[:, 0, 1] ** 2 + simp_pts[:, 0, 2] ** 2)
v2 = (simp_pts[:, 2, 0] ** 2 + simp_pts[:, 2, 1] ** 2 + simp_pts[:, 2, 2] ** 2) - (simp_pts[:, 0, 0] ** 2 + simp_pts[:, 0, 1] ** 2 + simp_pts[:, 0, 2] ** 2)
v3 = (simp_pts[:, 3, 0] ** 2 + simp_pts[:, 3, 1] ** 2 + simp_pts[:, 3, 2] ** 2) - (simp_pts[:, 0, 0] ** 2 + simp_pts[:, 0, 1] ** 2 + simp_pts[:, 0, 2] ** 2)

# solve a 3x3 system by inversion (see https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_3_%C3%97_3_matrices)
A = e*i - f*h
B = -(d*i - f*g)
C = d*h - e*g
D = -(b*i - c*h)
E = a*i - c*g
F = -(a*h - b*g)
G = b*f - c*e
H = -(a*f - c*d)
I = a*e - b*d

det = a*A + b*B + c*C

# multiplying inv*[v1, v2, v3] to get solution point (x, y, z)
x = (A*v1 + D*v2 + G*v3) / det
y = (B*v1 + E*v2 + H*v3) / det
z = (C*v1 + F*v2 + I*v3) / det

return (np.vstack((x, y, z))).T

我们想过滤掉原始曲面之外的四面体（参见上图中的长四面体）。 这可以通过在原始表面上测试四面体来完成。 然而，一个更简单的方法，它非常适合管/管表面的输入是 过滤掉具有大外接半径的四面体。 这就是alpha-shape 算法的作用。 这在我们的上下文中很容易做到，因为半径只是中心和任何四面体点之间的距离。

下图为过滤掉半径大于20的四面体后的Delaunay三角剖分。

我们现在可以使用这些构建块来构建满足半径条件的四面体的 Voronoi 子图。 下面的函数使用 Delaunay 三角剖分中的连通性信息来构建 Voronoi 子图，表示为边列表。

def compute_voronoi_vertices_and_edges(points, r_thresh=np.inf):
"""
Compute (finite) Voronoi edges and vertices of a set of points.
:param points: input points.
:param r_thresh: radius value for filtering out vertices corresponding to
Delaunay tetrahedrons with large radii of circumscribing sphere (alpha-shape condition).
:return: array of xyz Voronoi vertex points and an edge list.
"""
dt = Delaunay(points)
xyz_centers = compute_delaunay_tetra_circumcenters(dt)

# filtering tetrahedrons that have radius > thresh
simp_pts_0 = dt.points[dt.simplices[:, 0]]
radii = np.linalg.norm(xyz_centers - simp_pts_0, axis=1)
is_in = radii < r_thresh

# build an edge list from (filtered) tetrahedrons neighbor relations
edge_lst = []
for i in range(len(dt.neighbors)):
    if not is_in[i]:
        continue  # i is an outside tetra
    for j in dt.neighbors[i]:
        if j != -1 and is_in[j]:
            edge_lst.append((i, j))

return xyz_centers, edge_lst

结果仍然不够充分，如下图所示，其中子图边缘为黑色线段。 原因是 3D Delaunay 三角剖分因四面体薄而臭名昭著 （Shewchuk 在this paper 中所谓的 slivers, needles 和 caps）， 这会导致结果中出现外部“尖峰”。

虽然有一些通用方法可以消除这些不需要的尖峰 （例如，参见Amenta and Bern）， 在管表面的情况下，有一个简单的解决方案。 我们正在寻找的路径可以计算为图中的最短欧几里得路径，从距离管起点最近的点开始，到距离终点最近的点结束。 以下代码使用 networkx 图表执行此操作，权重设置为边的长度。

# get closest vertex to start and end points
xyz_centers, edge_lst = compute_voronoi_vertices_and_edges(pts, r_thresh=20.)
kdt = cKDTree(xyz_centers)
dist0, idx0 = kdt.query(np.array([100., 0, 0]))
dist1, idx1 = kdt.query(np.array([0, 100., 0]))

# compute shortest weighted path
edge_lengths = [np.linalg.norm(xyz_centers[e[0], :] - xyz_centers[e[1], :]) for e in edge_lst]
g = nx.Graph((i, j, {'weight': dist}) for (i, j), dist in zip(edge_lst, edge_lengths))
path_s = nx.shortest_path(g,source=idx0,target=idx1, weight='weight')

下图显示了原始200分的结果。

这是 1000 个点的密集样本的结果。

现在您可以通过路径点传递近似样条 - 插值或最小二乘拟合。 您可以按照链接中的建议使用scipy.interpolate.UnivariateSpline 或 scipy.interpolate.splrep 完成 here 或任何其他标准样条实现。