1、主要参考

（1）官网的介绍

Surface reconstruction — Open3D 0.16.0 documentation

（2）大佬的blog

三维点云重建 — open3d python_Coding的叶子的博客-CSDN博客_三维点云重建

（3）视频

Surface Reconstruction of Point Clouds - Point Cloud Processing in Open3D with_哔哩哔哩_bilibili

2、表面重建（Surface reconstruction）

在许多情况下，我们想要生成一个密集的3D几何图形，例如，一个三角形网格（triangle mesh）。然而，从多视点立体方法，或深度传感器，我们只能得到一个非结构化的点云。为了从这个非结构化输入中得到一个三角形网格，我们需要执行曲面重构。在文献中有一些方法，Open3D目前实现如下:

• Alpha shapes [Edelsbrunner1983]

• Ball pivoting [Bernardini1999]

• Poisson surface reconstruction [Kazhdan2006]

 2.1 alpha形状（Alpha shapes）

（1）原理

alpha shapes[Edelsbrunner1983]是凸包的一般化。正如这里所描述的，人们可以直观地认为alpha shapes如下:想象一个巨大的冰淇淋质量包含硬巧克力块点S。用这种球形的冰淇淋勺子，我们可以把冰淇淋块上所有我们够得到而不碰到巧克力块的部分都刻出来，甚至在里面也刻出了洞(例如，简单地从外面移动勺子就够不到的部分)。我们最终会得到一个(不一定是凸的)由帽、弧和点组成的对象。如果我们现在把所有的圆面拉直成三角形和线段，我们就有了一个直观的描述，叫做S的alpha shapes。

（2）

Alpha shapes 是一种散点外轮廓的提取方法。open3d中对应的函数为create_from_point_cloud_alpha_shape，其关键参数为alpha。alpha是该方法在搜索外轮廓时的半径大小。alpha值越小，网格的细节就越多，分辨率越高。  

 （3）主要函数

mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(pcd, alpha=2)

（4）测试代码

没错，还是那支兔子

import open3d as o3d
import numpy as npbunny = 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)

1）泊松磁盘采样点750的结果如下

2） alpha = 0.03时的重建结果

 该算法的实现基于点云的凸包。如果我们想从一个给定的点云计算多个alpha形状，那么我们可以通过只计算一次凸包并将其传递给create_from_point_cloud_alpha_shape来节省一些计算量。

 注意：下面不同alpha的测试代码

import open3d as o3d
import numpy as npbunny = 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)#以下方法节省一些计算量
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)

3）alpha = 0.5时的重建结果

 4）alpha = 0.136时的重建结果

5）alpha = 0.037时的重建结果

6）alpha = 0.01时的重建结果

  2.2球旋转（ Ball pivoting）

（1）原理

球旋转算法(BPA) [Bernardini1999]是一种与alpha形状相关的曲面重建方法。直观地，想象一个半径给定的3D球，我们把它扔在点云上。如果它击中任何3个点(并且它没有从这3个点掉下去)，它就会创建一个三角形。然后，算法开始从现有三角形的边缘旋转，每次它碰到3个球没有掉下去的点，我们就创建另一个三角形。

（2）函数

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

注意：该算法假设点云数据有法线。

 （3）输入还是兔子，泊松磁盘采样3000个点

import open3d as o3d
import numpy as npbunny = 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

radii = [0.005, 0.01, 0.02, 0.04]

（4）测试代码

import open3d as o3d
import numpy as npbunny = 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])

（5）测试结果

   2.3泊松曲面重建（ Poisson surface reconstruction）

2.3.1 泊松重建

（1）原理

泊松曲面重构方法[Kazhdan2006]解决正则化优化问题，获得光滑曲面。由于这个原因，泊松曲面重构可能比上面提到的方法更可取，因为它们会产生非光滑的结果，因为点云的点也是没有任何修改的三角形网格的顶点。

（2）函数

Open3D实现了create_from_point_cloud_poisson方法，它基本上是Kazhdan代码的包装器。函数的一个重要参数是depth，它定义了用于曲面重构的八叉树的深度，因此隐含了生成的三角形网格的分辨率。一个更高的深度值意味着一个网格有更多的细节。 

注意：该算法假设点云数据有法线。

 （3）测试的数据如下，比较威猛的老鹰：）

1）文件大小20多M

2）代码 

import open3d as o3d
import numpy as npeagle = 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])

 3）实物图

 （4）重建测试代码

import open3d as o3d
import numpy as npeagle = 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])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])

（5）重建后果然很不错

2.3.2 低密度值

（1）描述 

 泊松曲面重构也会在低密度的区域创建三角形，甚至外推到一些区域(见上图老鹰输出的底部)。create_from_point_cloud_poisson函数有第二个密度返回值，表示每个顶点的密度。低密度值意味着顶点仅由来自输入点云的少量点支持。

（2） 在下面的代码中，我们使用伪颜色在3D中可视化密度。紫色表示低密度，黄色表示高密度。

import open3d as o3d
import numpy as np
from matplotlib import pyplot as plt##(1)input 
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])##(2)Poisson surface reconstruction
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])##(3) Violet indicates low density and yellow indicates a high density
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])

（3）测试结果

 （4）我们可以进一步使用密度值来删除支持度低的顶点和三角形。在下面的代码中，我们删除密度值低于所有密度值的0.01分位数的所有顶点(和连接的三角形)。

import open3d as o3d
import numpy as np
from matplotlib import pyplot as plt##(1)input 
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])##(2)Poisson surface reconstruction
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])##(3) Violet indicates low density and yellow indicates a high density
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])##(4) remove low density vertices
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）移除低密度值后的结果

 2.4 法向估计（Normal estimation）

（1）描述

在上面的例子中，我们假设点云有指向外的法线。然而，并不是所有的点云都有相关的法线。Open3D可以使用estimate_normals来估计点云法线，它在每个3D点的局部拟合一个平面来推导法线。然而，估计的法线可能不是一致的方向。Orient_normals_consistent_tangent_plane使用最小生成树传播法线方向。

PS注意：Orient_normals_consistent_tangent_plane要在estimate_normals来估计点云法线的基础上使用。

（2）estimate_normals获得法向的测试代码

1）代码

import open3d as o3d
import numpy as np
from matplotlib import pyplot as plt##(1)estimate_normals 
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 normalspcd.estimate_normals()
o3d.visualization.draw_geometries([pcd], point_show_normal=True)

2）测试结果

 （3）使用orient_normals_consistent_tangent_plane获得法向

1）代码

import open3d as o3d
import numpy as np
from matplotlib import pyplot as plt##(1)estimate_normals 
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 normalspcd.estimate_normals()
o3d.visualization.draw_geometries([pcd], point_show_normal=True)##(2)orient_normals_consistent_tangent_plane 
pcd.orient_normals_consistent_tangent_plane(100)
o3d.visualization.draw_geometries([pcd], point_show_normal=True)

2）结果