【小白深度教程 1.20】手把手教你使用 Open3D（2）点云聚类、分割和重建（含 Python 代码）
1. 表面重建概念
在许多场景中，我们需要生成密集的 3D 几何体，例如三角网格。然而，从多视图立体方法或深度传感器中，我们只能获得非结构化的点云。为了从这些非结构化输入中生成三角网格，我们需要执行表面重建。在文献中有多种方法，而 Open3D 目前实现了以下几种：
- Alpha 形状 [Edelsbrunner1983]
- 球形旋转 [Bernardini1999]
- 泊松表面重建 [Kazhdan2006]
2. Alpha 形状
Alpha 形状 [Edelsbrunner1983] 是凸包的泛化。正如这里所描述的，可以直观地将 Alpha 形状理解为：想象一大块包含点集 S 的冰激凌，这些点是硬质巧克力块。使用一个球形的冰激凌勺，我们挖出所有能够在不碰到巧克力块的情况下从冰激凌块中取出的部分，甚至挖出内部的孔洞（例如，不能从外部简单移动勺子就能到达的部分）。最终得到的是一个由盖帽、弧和点界定的对象（不一定是凸的）。如果我们将所有的圆形面拉直为三角形和线段，那么我们就有了点集 S 的 Alpha 形状的直观描述。
Open3D 实现了 create_from_point_cloud_alpha_shape 方法，该方法涉及折中参数 alpha。
mesh = o3dtut.get_bunny_mesh()
pcd = mesh.sample_points_poisson_disk(750)
o3d.visualization.draw_geometries([pcd])
alpha = 0.03
print(f"alpha={alpha:.3f}")
mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(pcd, alpha)
mesh.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh], mesh_show_back_face=True)
点云数据：
调整 alpha = 0.03：
该实现基于点云的凸包。如果我们希望从给定的点云计算多个 Alpha 形状，那么只需计算一次凸包并将其传递给 create_from_point_cloud_alpha_shape，即可节省一些计算量。
tetra_mesh, pt_map = o3d.geometry.TetraMesh.create_from_point_cloud(pcd)
for alpha in np.logspace(np.log10(0.5), np.log10(0.01), num=4):
print(f"alpha={alpha:.3f}")
mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(
pcd, alpha, tetra_mesh, pt_map)
mesh.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh], mesh_show_back_face=True)
调整 alpha = 0.5：
调整 alpha = 0.136：
调整 alpha = 0.037：
调整 alpha = 0.01：
3. 球形旋转
球形旋转算法（BPA）[Bernardini1999] 是一种与 Alpha 形状相关的表面重建方法。直观地讲，可以想象一个给定半径的 3D 球体，我们将其放到点云上。如果球体击中任意 3 个点（且不会从这 3 个点之间掉落），它就会创建一个三角形。然后，算法从已有三角形的边缘开始旋转，每当球体击中 3 个点且不会掉落时，就会创建另一个三角形。
Open3D 在 create_from_point_cloud_ball_pivoting 中实现了此方法。该方法接受一个半径列表作为参数，这些半径对应于在点云上旋转的各个球体的半径。
注意：
该算法假设 PointCloud 具有法线。
gt_mesh = o3dtut.get_bunny_mesh()
gt_mesh.compute_vertex_normals()
pcd = gt_mesh.sample_points_poisson_disk(3000)
o3d.visualization.draw_geometries([pcd])
radii = [0.005, 0.01, 0.02, 0.04]
rec_mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting(
pcd, o3d.utility.DoubleVector(radii))
o3d.visualization.draw_geometries([pcd, rec_mesh])
4. 泊松表面重建
泊松表面重建方法 [Kazhdan2006] 通过解决正则化优化问题来获得光滑的表面。因此，相比于前面提到的方法，泊松表面重建更具优势，因为前述方法会产生非光滑的结果，点云中的点直接成为结果三角网格的顶点而未进行任何修改。
Open3D 实现了 create_from_point_cloud_poisson 方法，该方法基本上是 Kazhdan 代码的封装。函数的一个重要参数是 depth，它定义了用于表面重建的八叉树的深度，从而影响最终三角网格的分辨率。较高的深度值意味着生成的网格具有更多细节。
注意：
该算法假设 PointCloud 具有法线。
pcd = o3dtut.get_eagle_pcd()
print(pcd)
o3d.visualization.draw_geometries([pcd],
zoom=0.664,
front=[-0.4761, -0.4698, -0.7434],
lookat=[1.8900, 3.2596, 0.9284],
up=[0.2304, -0.8825, 0.4101])
print('run Poisson surface reconstruction')
with o3d.utility.VerbosityContextManager(
o3d.utility.VerbosityLevel.Debug) as cm:
mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(
pcd, depth=9)
print(mesh)
o3d.visualization.draw_geometries([mesh],
zoom=0.664,
front=[-0.4761, -0.4698, -0.7434],
lookat=[1.8900, 3.2596, 0.9284],
up=[0.2304, -0.8825, 0.4101])
输出：
run Poisson surface reconstruction
[Open3D DEBUG] Input Points / Samples: 796825 / 368254
[Open3D DEBUG] # Got kernel density: 0.170430 (s), 295.309 (MB) / 295.309 (MB) / 328 (MB)
[Open3D DEBUG] # Got normal field: 1.0846 (s), 391.949 (MB) / 391.949 (MB) / 391 (MB)
[Open3D DEBUG] Point weight / Estimated Area: 2.623552e-06 / 2.090512e+00
[Open3D DEBUG] # Finalized tree: 1.1184 (s), 516.758 (MB) / 516.758 (MB) / 516 (MB)
[Open3D DEBUG] # Set FEM constraints: 3.67671 (s), 470.125 (MB) / 516.758 (MB) / 516 (MB)
[Open3D DEBUG] #Set point constraints: 0.517555 (s), 470.125 (MB) / 516.758 (MB) / 516 (MB)
[Open3D DEBUG] Leaf Nodes / Active Nodes / Ghost Nodes: 2945433 / 3365000 / 1209
[Open3D DEBUG] Memory Usage: 470.125 MB
[Open3D DEBUG] # Linear system solved: 4.45958 (s), 504.602 (MB) / 516.758 (MB) / 516 (MB)
[Open3D DEBUG] Got average: 0.223773 (s), 504.602 (MB) / 516.758 (MB) / 516 (MB)
[Open3D DEBUG] Iso-Value: 5.028479e-01 = 4.006817e+05 / 7.968250e+05
[Open3D DEBUG] # Total Solve: 21.2 (s), 656.6 (MB)
TriangleMesh with 563112 points and 1126072 triangles.
泊松表面重建还会在点密度较低的区域生成三角形，甚至会在某些区域进行外推（参见上方鹰形输出的底部）。create_from_point_cloud_poisson 函数有一个第二个返回值 densities，用于指示每个顶点的密度。低密度值表示该顶点仅由输入点云中的少量点支持。
在下面的代码中，我们使用伪彩色在 3D 中可视化密度。紫色表示低密度，而黄色表示高密度。
print('visualize densities')
densities = np.asarray(densities)
density_colors = plt.get_cmap('plasma')(
(densities - densities.min()) / (densities.max() - densities.min()))
density_colors = density_colors[:, :3]
density_mesh = o3d.geometry.TriangleMesh()
density_mesh.vertices = mesh.vertices
density_mesh.triangles = mesh.triangles
density_mesh.triangle_normals = mesh.triangle_normals
density_mesh.vertex_colors = o3d.utility.Vector3dVector(density_colors)
o3d.visualization.draw_geometries([density_mesh],
zoom=0.664,
front=[-0.4761, -0.4698, -0.7434],
lookat=[1.8900, 3.2596, 0.9284],
up=[0.2304, -0.8825, 0.4101])
可视化密度：
我们可以进一步利用密度值去除支撑较少的顶点和三角形。在下面的代码中，我们移除了所有密度值低于所有密度值的 0.01 分位数的顶点（以及相连的三角形）。
print('remove low density vertices')
vertices_to_remove = densities < np.quantile(densities, 0.01)
mesh.remove_vertices_by_mask(vertices_to_remove)
print(mesh)
o3d.visualization.draw_geometries([mesh],
zoom=0.664,
front=[-0.4761, -0.4698, -0.7434],
lookat=[1.8900, 3.2596, 0.9284],
up=[0.2304, -0.8825, 0.4101])
5. 法线估计
在上述示例中，我们假设点云具有指向外部的法线。然而，并非所有点云都自带关联的法线。Open3D 可以使用 estimate_normals 来估计点云法线，该方法通过对每个 3D 点局部拟合一个平面来推导法线。然而，估计的法线方向可能并不一致。orient_normals_consistent_tangent_plane 使用最小生成树传播法线方向，以确保一致性。
gt_mesh = o3dtut.get_bunny_mesh()
pcd = gt_mesh.sample_points_poisson_disk(5000)
pcd.normals = o3d.utility.Vector3dVector(np.zeros(
(1, 3))) # invalidate existing normals
pcd.estimate_normals()
o3d.visualization.draw_geometries([pcd], point_show_normal=True)
pcd.orient_normals_consistent_tangent_plane(100)
o3d.visualization.draw_geometries([pcd], point_show_normal=True)
扫一扫