表面重建

在许多场景下，我们希望生成密集的三维几何，即三角网格（triangle mesh）。然而，从多视图立体（multi-view stereo）方法或深度传感器中，我们只能获得非结构化的点云（point cloud）。要从这种非结构化的输入得到三角网格，就需要进行表面重建（surface reconstruction）。文献中已有多种方法，Open3D 目前实现了以下几种：

Alpha shapes（Alpha 形状）[Edelsbrunner1983]
Ball pivoting（滚球法）[Bernardini1999]
Poisson surface reconstruction（泊松表面重建）[Kazhdan2006]

Alpha shapes

Alpha shape（Alpha 形状）[Edelsbrunner1983] 是凸包（convex hull）的推广。正如这里所描述的，可以直观地把 alpha shape 想象成下面这样：假设有一大块冰淇淋，里面嵌着以点集 $S$ 表示的硬巧克力碎块。我们用一把球形的冰淇淋勺，挖掉所有能够触及、又不会碰到巧克力碎块的冰淇淋部分，甚至可以在内部挖出洞（也就是仅从外部移动勺子无法触及的部分）。最终我们会得到一个由球冠、圆弧和点围成的（不一定是凸的）物体。如果我们把所有弯曲的表面都拉直成三角形和线段，就得到了点集 $S$ 的 alpha shape 的直观描述。

Open3D 实现了 create_from_point_cloud_alpha_shape 方法，其中包含一个权衡参数 alpha。

bunny = o3d.data.BunnyMesh()
mesh = o3d.io.read_triangle_mesh(bunny.path)
mesh.compute_vertex_normals()

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 shape，可以先只计算一次凸包，再把它传给 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)

Ball pivoting

滚球算法（Ball Pivoting Algorithm, BPA）[Bernardini1999] 是一种与 alpha shapes 相关的表面重建方法。直观地说，想象一个给定半径的三维球落到点云上：如果它碰到任意 3 个点（并且不会从这 3 个点之间漏下去），就会生成一个三角形。随后，算法从已有三角形的边开始转动（pivot）；每当它碰到 3 个点且球不会漏下去时，就再生成一个三角形。

Open3D 在 create_from_point_cloud_ball_pivoting 中实现了该方法。该方法接受一个半径列表作为参数，对应于在点云上滚动的各个球的半径。

bunny = o3d.data.BunnyMesh()
gt_mesh = o3d.io.read_triangle_mesh(bunny.path)
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])

泊松表面重建

泊松表面重建（Poisson surface reconstruction）方法 [Kazhdan2006] 通过求解一个正则化的优化问题来获得光滑的表面。正因如此，泊松表面重建往往优于前面提到的那些方法——后者直接把点云中的点当作结果网格的顶点而不做任何修改，因此得到的表面并不光滑。

Open3D 实现了 create_from_point_cloud_poisson 方法，它本质上是 Kazhdan 代码的封装。该函数的一个重要参数是 depth，它定义了用于表面重建的八叉树（octree）深度，从而决定了结果三角网格的分辨率。depth 值越高，网格的细节就越丰富。

eagle = o3d.data.EaglePointCloud()
pcd = o3d.io.read_point_cloud(eagle.path)

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])

输出：

[Open3D INFO] Downloading https://github.com/isl-org/open3d_downloads/releases/download/20220201-data/EaglePointCloud.ply
[Open3D INFO] Downloaded to /home/runner/open3d_data/download/EaglePointCloud/EaglePointCloud.ply
PointCloud with 796825 points.

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])

在 Debug 日志级别下，可以看到泊松求解器的进展：

run Poisson surface reconstruction
[Open3D DEBUG] Input Points / Samples: 796825 / 368254
[Open3D DEBUG] #   Got kernel density: 0.08908891677856445 (s), 392.31640625 (MB) / 392.31640625 (MB) / 476 (MB)
[Open3D DEBUG] #     Got normal field: 0.6136579513549805 (s), 490.06640625 (MB) / 490.06640625 (MB) / 490 (MB)
[Open3D DEBUG] Point weight / Estimated Area: 2.623551e-06 / 2.090511e+00
[Open3D DEBUG] #       Finalized tree: 0.48689985275268555 (s), 615.74609375 (MB) / 615.74609375 (MB) / 615 (MB)
[Open3D DEBUG] #  Set FEM constraints: 1.4530129432678223 (s), 576.54296875 (MB) / 615.74609375 (MB) / 615 (MB)
[Open3D DEBUG] #Set point constraints: 0.23029398918151855 (s), 576.54296875 (MB) / 615.74609375 (MB) / 615 (MB)
[Open3D DEBUG] Leaf Nodes / Active Nodes / Ghost Nodes: 2945433 / 3365000 / 1209
[Open3D DEBUG] Memory Usage: 576.543 MB
[Open3D DEBUG] # Linear system solved: 1.8742828369140625 (s), 615.16796875 (MB) / 615.74609375 (MB) / 615 (MB)
[Open3D DEBUG] Got average: 0.06455492973327637 (s), 615.16796875 (MB) / 615.74609375 (MB) / 615 (MB)
[Open3D DEBUG] Iso-Value: 5.028478e-01 = 4.006817e+05 / 7.968250e+05
[Open3D DEBUG] #          Total Solve:       9.2 (s),     759.9 (MB)
TriangleMesh with 563112 points and 1126072 triangles.

泊松表面重建也会在点密度较低的区域生成三角形，甚至会向某些区域外推（参见上面鹰模型输出的底部）。create_from_point_cloud_poisson 函数还有第二个返回值 densities，它为每个顶点指示其密度。密度值越低，说明该顶点仅由输入点云中较少的点所支撑。

下面的代码中，我们用伪彩色（pseudo color）在三维中可视化密度。紫色表示低密度，黄色表示高密度。

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$ 分位数（quantile）的顶点（及其相连的三角形）。

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])

输出：

remove low density vertices
TriangleMesh with 557480 points and 1113212 triangles.

法向量估计

在上面的例子中，我们假设点云已经带有朝外的法向量。然而，并非所有点云都自带法向量。Open3D 可以用 estimate_normals 来估计点云法向量，该方法为每个三维点局部拟合一个平面以推导法向量。但估计出的法向量朝向可能并不一致。orient_normals_consistent_tangent_plane 会利用最小生成树来传播法向量的朝向。

bunny = o3d.data.BunnyMesh()
gt_mesh = o3d.io.read_triangle_mesh(bunny.path)

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)