Python Open3D几何图形基础篇（二）网格操作

2022年 11月 16日

Python Open3D几何图形基础篇（二）网格操作

本文所参考网页：
Mesh — Open3D 0.15.1 documentation

Open3D有一个用于3D三角网格的数据结构叫TriangleMesh 。

下面的代码显示了如何从ply 文件中读取三角网格，并打印其顶点和三角形。

print("Testing mesh in Open3D...")
armadillo_mesh = o3d.data.ArmadilloMesh()
mesh = o3d.io.read_triangle_mesh(armadillo_mesh.path)

knot_mesh = o3d.data.KnotMesh()
mesh = o3d.io.read_triangle_mesh(knot_mesh.path)
print(mesh)
print('Vertices:')
print(np.asarray(mesh.vertices))
print('Triangles:')
print(np.asarray(mesh.triangles))
"""
Testing mesh in Open3D...
[Open3D INFO] Downloading https://github.com/isl-org/open3d_downloads/releases/download/20220201-data/KnotMesh.ply
[Open3D INFO] Downloaded to /home/runner/open3d_data/download/KnotMesh/KnotMesh.ply
TriangleMesh with 1440 points and 2880 triangles.
Vertices:
[[  4.51268387  28.68865967 -76.55680847]
 [  7.63622284  35.52046967 -69.78063965]
 [  6.21986008  44.22465134 -64.82303619]
 ...
 [-22.12651634  31.28466606 -87.37570953]
 [-13.91188431  25.4865818  -86.25827026]
 [ -5.27768707  23.36245346 -81.43279266]]
Triangles:
[[   0   12   13]
 [   0   13    1]
 [   1   13   14]
 ...
 [1438   11 1439]
 [1439   11    0]
 [1439    0 1428]]
"""

TriangleMesh 类有几个数据字段，例如vertices 和triangles。Open3D通过numpy提供对这些字段的直接内存访问。

3D网格可视化（Visualize a 3D mesh）

print("Try to render a mesh with normals (exist: " +
      str(mesh.has_vertex_normals()) + ") and colors (exist: " +
      str(mesh.has_vertex_colors()) + ")")
o3d.visualization.draw_geometries([mesh])
print("A mesh with no normals and no colors does not look good.")

# Try to render a mesh with normals (exist: False) and colors (exist: False)
# 下图
# A mesh with no normals and no colors does not look good.

3D网格可视化产生的图片

你可以旋转和移动这个网格，但因为它是统一用灰色绘制的，所以它看起来不是3d的。
原因是因为当前所看到的网格没有顶点或面的法线。所以，使用的是统一的颜色上色，而不是更复杂的Phong着色（Phong shading）。【~~并不是很懂这句话~~】

表面法线估计（Surface normal estimation）

使用表面法线绘制网格

print("Computing normal and rendering it.")
mesh.compute_vertex_normals()
print(np.asarray(mesh.triangle_normals))
o3d.visualization.draw_geometries([mesh])

"""
Computing normal and rendering it.
[[ 0.79164373 -0.53951444  0.28674793]
 [ 0.8319824  -0.53303008  0.15389681]
 [ 0.83488162 -0.09250101  0.54260136]
 ...
 [ 0.16269924 -0.76215917 -0.6266118 ]
 [ 0.52755226 -0.83707495 -0.14489352]
 [ 0.56778973 -0.76467734 -0.30476777]]
"""

表面法线估计的效果图

这里所使用的compute_vertex_normals 与paint_uniform_color 是mesh 的成员函数。

裁剪网格（Crop mesh）

我们通过直接操作网格的triangle 和triangle_normals数据字段来去除一半曲面。这个操作是通过numpy实现的。

print("We make a partial mesh of only the first half triangles.")
mesh1 = copy.deepcopy(mesh)
mesh1.triangles = o3d.utility.Vector3iVector(
    np.asarray(mesh1.triangles)[:len(mesh1.triangles) // 2, :])
mesh1.triangle_normals = o3d.utility.Vector3dVector(
    np.asarray(mesh1.triangle_normals)[:len(mesh1.triangle_normals) // 2, :])
print(mesh1.triangles)
o3d.visualization.draw_geometries([mesh1])

"""
We make a partial mesh of only the first half triangles.
std::vector<Eigen::Vector3i> with 1440 elements.
Use numpy.asarray() to access data.
"""

裁剪网格中的图片

网格上色（Paint mesh）

paint_uniform_color 可以为网格涂上统一的颜色。颜色为RGB空间，范围是0到1。

print("Painting the mesh")
mesh1.paint_uniform_color([1, 0.706, 0])
o3d.visualization.draw_geometries([mesh1])

网格上色生成的图片

网格属性（Mesh properties）

三角网格有一些可以用Open3D检测的属性。

一个重要的属性是流形属性，我们可以测试三角网格是否是边流形，是否是边流形（is_edge_manifold ），是否是顶点流形（is_vertex_manifold）。如果每条边都是一个或两个三角形的边界，那么这个三角网格是边流形的。

函数is_edge_manifold 具有bool 值的参数allow_boundary_edges ，用于定义是否应允许边界边（if boundary edges should be allowed.）。

此外，如果顶点的星形（the star of the vertex）是边流形的且是边连接的（edge-connected），则三角网格是顶点流形的。例如，两个或多个面仅由一个顶点而非一条边连接。

另一个属性是测试自相交。如果网格中存在与另一个网格相交的三角形，则函数is_self_intersecting 返回True 。

不漏水网格（watertight mesh）可以被定义为是边流形的、顶点流形的且不自交的网格。在Open3D中可以通过函数is_watertight 来实现这个检测。

如果网格是可定向的，我们也可以测试三角网格的这个属性。可定向的，也就是说，三角形的定向方式可以使所有法线指向外边（the triangles can be oriented in such a way that all normals point towards
the outside. ）。Open3D中对应的函数叫is_orientable。

下面的代码根据这些属性测试了一些三角网格，并将结果可视化。非流形边显示为红色，边界边显示为绿色，非流形顶点显示为绿色点，自交三角显示为粉色。

示例代码：

# o3dtut 是官方代码open3d_tutorial.py 
# 也就是 import open3d_tutorial as o3dtut 
# 
def check_properties(name, mesh):
    mesh.compute_vertex_normals()

    edge_manifold = mesh.is_edge_manifold(allow_boundary_edges=True)
    edge_manifold_boundary = mesh.is_edge_manifold(allow_boundary_edges=False)
    vertex_manifold = mesh.is_vertex_manifold()
    self_intersecting = mesh.is_self_intersecting()
    watertight = mesh.is_watertight()
    orientable = mesh.is_orientable()

    print(name)
    print(f"  edge_manifold:          {edge_manifold}")
    print(f"  edge_manifold_boundary: {edge_manifold_boundary}")
    print(f"  vertex_manifold:        {vertex_manifold}")
    print(f"  self_intersecting:      {self_intersecting}")
    print(f"  watertight:             {watertight}")
    print(f"  orientable:             {orientable}")

    geoms = [mesh]
    if not edge_manifold:
        edges = mesh.get_non_manifold_edges(allow_boundary_edges=True)
        geoms.append(o3dtut.edges_to_lineset(mesh, edges, (1, 0, 0)))
    if not edge_manifold_boundary:
        edges = mesh.get_non_manifold_edges(allow_boundary_edges=False)
        geoms.append(o3dtut.edges_to_lineset(mesh, edges, (0, 1, 0)))
    if not vertex_manifold:
        verts = np.asarray(mesh.get_non_manifold_vertices())
        pcl = o3d.geometry.PointCloud(
            points=o3d.utility.Vector3dVector(np.asarray(mesh.vertices)[verts]))
        pcl.paint_uniform_color((0, 0, 1))
        geoms.append(pcl)
    if self_intersecting:
        intersecting_triangles = np.asarray(
            mesh.get_self_intersecting_triangles())
        intersecting_triangles = intersecting_triangles[0:1]
        intersecting_triangles = np.unique(intersecting_triangles)
        print("  # visualize self-intersecting triangles")
        triangles = np.asarray(mesh.triangles)[intersecting_triangles]
        edges = [
            np.vstack((triangles[:, i], triangles[:, j]))
            for i, j in [(0, 1), (1, 2), (2, 0)]
        ]
        edges = np.hstack(edges).T
        edges = o3d.utility.Vector2iVector(edges)
        geoms.append(o3dtut.edges_to_lineset(mesh, edges, (1, 0, 1)))
    o3d.visualization.draw_geometries(geoms, mesh_show_back_face=True)


knot_mesh_data = o3d.data.KnotMesh()
knot_mesh = o3d.io.read_triangle_mesh(knot_mesh_data.path)
check_properties('KnotMesh', knot_mesh)
check_properties('Mobius', o3d.geometry.TriangleMesh.create_mobius(twists=1))
check_properties("non-manifold edge", o3dtut.get_non_manifold_edge_mesh())
check_properties("non-manifold vertex", o3dtut.get_non_manifold_vertex_mesh())
check_properties("open box", o3dtut.get_open_box_mesh())
check_properties("intersecting_boxes", o3dtut.get_intersecting_boxes_mesh())

这里只显示2个例子的图具体可以上官方手册看
网格属性效果图1
网格属性效果图2

网格滤波器（Mesh filtering）

Open3D中有一些网格滤波器的方法，下面将展示实现了的滤波器用它来平滑有噪声的三角形网格。

平均滤波器（Average filter）

最简单的滤波器是平均滤波器。一个给定的顶点

v_i

$v_{i}$ 由相邻顶点

\mathcal{N}

$N$ 的平均值给出。

∑

∈

∣

v_i = \frac{v_i + { \sum_{n\in\mathcal{N}}^{}}v_n}{ \left | N \right | + 1 }

$v_{i} = \frac{v _{i} + \sum _{n \in N} v _{n}}{∣ N ∣ + 1}$

如下面的代码所示，这个滤波器可用于对网格进行去噪。在函数filter_smooth_simple 中的参数number_of_iterations 定义的是对网格使用滤波器的次数。

print('create noisy mesh')
knot_mesh = o3d.data.KnotMesh()
mesh_in = o3d.io.read_triangle_mesh(knot_mesh.path)
vertices = np.asarray(mesh_in.vertices)
noise = 5
vertices += np.random.uniform(0, noise, size=vertices.shape)
mesh_in.vertices = o3d.utility.Vector3dVector(vertices)
mesh_in.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_in])

print('filter with average with 1 iteration')
mesh_out = mesh_in.filter_smooth_simple(number_of_iterations=1)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])

print('filter with average with 5 iterations')
mesh_out = mesh_in.filter_smooth_simple(number_of_iterations=5)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])

create noisy mesh
生成有噪声的网格
filter with average with 1 iteration
迭代了1次的效果图
filter with average with 5 iterations
递归了3次的效果图

拉普拉斯滤波器（Laplacian）

另一个重要的网格滤波器叫拉普拉斯滤波器，它是这样定义的：

v

i

=

v

i

⋅

λ

∑

n

∈

N

w

n

v

n

−

v

i

v_i = v_i \cdot \lambda \sum_{n \in N }^{} w_nv_n-v_i

$v_{i} = v_{i} \cdot λ n \in N \sum w_{n} v_{n} - v_{i}$

其中，

λ

\lambda

$λ$ 是滤波器的强度

w

n

w_n

$w_{n}$ 是与相邻顶点的距离相关的归一化权重。这个滤波器由函数filter_smooth_laplacian 实现，有参数number_of_iterations 和lambda 。

示例代码如下：
```
print('filter with Laplacian with 10 iterations')
mesh_out = mesh_in.filter_smooth_laplacian(number_of_iterations=10)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])

print('filter with Laplacian with 50 iterations')
mesh_out = mesh_in.filter_smooth_laplacian(number_of_iterations=50)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])
```
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
filter with Laplacian with 10 iterations

filter with Laplacian with 50 iterations
Taubin 滤波器（Taubin filter）

平均滤波器和拉普拉斯滤波器存在的问题是它们都会导致三角网格收缩。Taubin1995表明，使用两个不同

λ

\lambda

$λ$ 参数的拉普拉斯滤波器可以防止网格收缩。filter_smooth_taubin这个函数实现了这个滤波器。

示例代码如下：
```
print('filter with Taubin with 10 iterations')
mesh_out = mesh_in.filter_smooth_taubin(number_of_iterations=10)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])

print('filter with Taubin with 100 iterations')
mesh_out = mesh_in.filter_smooth_taubin(number_of_iterations=100)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])
```
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
filter with Taubin with 10 iterations

filter with Taubin with 100 iterations

采样（Sampling）

Open3D中有从三角网格采样点云的函数。最简单的方法是sample_points_uniformly ，它基于三角形区域从三维曲面均匀采样点。

代码如下：

mesh = o3d.geometry.TriangleMesh.create_sphere()
mesh.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh])
pcd = mesh.sample_points_uniformly(number_of_points=500)
o3d.visualization.draw_geometries([pcd])

采样前
采样后

另外还有一个采样兔子的样例效果差不多就不赘述了。

均匀采样能在曲面上生成点簇（yield clusters of points），然而泊松圆盘采样可以在曲面上均匀分布（evenly distribute）点。

方法 sample_points_poisson_disk 实现了样本消除（sample elimination），它从采样点云开始，移除点以满足采样标准。

这个方法支持两种选项来提供初始点云：
1. 默认使用参数init_factor：这个方法首先通过 init_factor x number_of_points 的方式从网格均匀采样点云，然后将其用于消除；
2. 可以提供一个点云，并将其传递给 sample_points_poisson_disk方法。然后，使用该点云进行消除。

示例代码如下：

mesh = o3d.geometry.TriangleMesh.create_sphere()
pcd = mesh.sample_points_poisson_disk(number_of_points=500, init_factor=5)
o3d.visualization.draw_geometries([pcd])

pcd = mesh.sample_points_uniformly(number_of_points=2500)
pcd = mesh.sample_points_poisson_disk(number_of_points=500, pcl=pcd)
o3d.visualization.draw_geometries([pcd])

均匀采样效果图
对上图用泊松圆盘采样的效果图

这里也是还有一个采样兔子的样例效果差不多就不赘述了。（~~这个方法我没看出什么效果~~ ）

网格再分割（Mesh subdivision）

在网格再分割中，我们将每个三角形划分为多个较小的三角形。

在最简单的情况下，我们计算每个三角形每条边的中点，并将三角形分成四个较小的三角形。这个方法在函数subdivide_midpoint 中实现了。

3D曲面和面积保持不变，但其顶点和三角形的数量增加了。参数number_of_iterations 定义了这个过程重复了多少次。

mesh = o3d.geometry.TriangleMesh.create_box()
mesh.compute_vertex_normals()
print(
    f'The mesh has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles'
)
o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True)
mesh = mesh.subdivide_midpoint(number_of_iterations=1)
print(
    f'After subdivision it has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles'
)
o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True)

The mesh has 8 vertices and 12 triangles
使用前

After subdivision it has 26 vertices and 48 triangles
使用后

Open3D基于[Loop1987]实现了一种额外的再分割方法。这个方法基于四次长方体样条曲线，该样条曲线在任何地方生成 $C^2 C2连续极限曲面(continuous limit surfaces)，但在 C 1 C^1 C1连续的特殊顶点除外。这将导致更平滑的拐角。$

示例代码如下：

mesh = o3d.geometry.TriangleMesh.create_sphere()
mesh.compute_vertex_normals()
print(
    f'The mesh has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles'
)
o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True)
mesh = mesh.subdivide_loop(number_of_iterations=2)
print(
    f'After subdivision it has {len(mesh.vertices)} vertices and {len(mesh.triangles)} triangles'
)
o3d.visualization.draw_geometries([mesh], zoom=0.8, mesh_show_wireframe=True)

The mesh has 762 vertices and 1520 triangles
使用前
After subdivision it has 12162 vertices and 24320 triangles
使用后
另外还有一个采样结的样例效果差不多就不赘述了。

网格简化（Mesh simplification）

有时，我们希望用较少的三角形和顶点表示高分辨率网格，但低分辨率网格仍应接近高分辨率网格。为此，Open3D实现了一些网格简化方法。

顶点聚类（Vertex clustering）

顶点聚类方法将落入给定大小的体素的所有顶点汇集到单个顶点。

这个方法是通过simplify_vertex_clustering 实现的，并且有两个参数voxel_size 与contraction ，voxel_size 定义了定义体素栅格（voxel grid）的大小，contraction 定义了顶点是如何汇集的。o3d.geometry.SimplificationContraction.Average 用于计算简单的平均数。

示例代码如下

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

print(
    f'Input mesh has {len(mesh_in.vertices)} vertices and {len(mesh_in.triangles)} triangles'
)
o3d.visualization.draw_geometries([mesh_in])

voxel_size = max(mesh_in.get_max_bound() - mesh_in.get_min_bound()) / 32
print(f'voxel_size = {voxel_size:e}')
mesh_smp = mesh_in.simplify_vertex_clustering(
    voxel_size=voxel_size,
    contraction=o3d.geometry.SimplificationContraction.Average)
print(
    f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles'
)
o3d.visualization.draw_geometries([mesh_smp])

voxel_size = max(mesh_in.get_max_bound() - mesh_in.get_min_bound()) / 16
print(f'voxel_size = {voxel_size:e}')
mesh_smp = mesh_in.simplify_vertex_clustering(
    voxel_size=voxel_size,
    contraction=o3d.geometry.SimplificationContraction.Average)
print(
    f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles'
)
o3d.visualization.draw_geometries([mesh_smp])

Input mesh has 1440 vertices and 2880 triangles

voxel_size = 5.650959e+00
Simplified mesh has 1355 vertices and 2720 triangles

voxel_size = 1.130192e+01
Simplified mesh has 860 vertices and 1773 triangles

网格抽取（Mesh decimation）
- 另一种网格简化的方法是以增加步数（incremental steps）的方式来进行的网格抽取。这个算法将会选择一个三角形，使误差度量（error metric）最小化，并将其删除。重复此操作，直到达到所需数量的三角形。
- Open3D 通过 simplify_quadric_decimation 实现了这个算法，也就最小化误差二次曲面（到相邻平面的距离）。参数target_number_of_triangles 定义了抽取算法的停止标准。
```
mesh_smp = mesh_in.simplify_quadric_decimation(target_number_of_triangles=6500)
print(
    f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles'
)
o3d.visualization.draw_geometries([mesh_smp])

mesh_smp = mesh_in.simplify_quadric_decimation(target_number_of_triangles=1700)
print(
    f'Simplified mesh has {len(mesh_smp.vertices)} vertices and {len(mesh_smp.triangles)} triangles'
)
o3d.visualization.draw_geometries([mesh_smp])
```
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
Simplified mesh has 1440 vertices and 2880 triangles

Simplified mesh has 850 vertices and 1700 triangles

连接组件（Connected components）

各种许多重建方法的结果（~~不是很懂这句话的意思~~ The result of various reconstruction methods. ）。

Open3D 通过cluster_connected_triangles 函数实现了连接组件的算法，这个算法将每个三角形分配给一簇（a cluster）连接了的三角形。这个函数将会为每个三角形返回三个值triangle_clusters 、cluster_n_triangles 、cluster_area。

triangle_clusters ：簇的索引值

cluster_n_triangles ：每个簇中三角形的数量

cluster_area：簇的曲面表面积

比如，这在RGBD Integration中很有用，RGBD Integration并不常是一个三角网格，而是多个网格。一些较小的部分是由于噪音造成的，我们很可能希望将其移除。

以下的代码展示了cluster_connected_triangles 的使用与以及如何用它去除错误的三角形。

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

mesh = mesh.subdivide_midpoint(number_of_iterations=2)
vert = np.asarray(mesh.vertices)
min_vert, max_vert = vert.min(axis=0), vert.max(axis=0)
for _ in range(30):
    cube = o3d.geometry.TriangleMesh.create_box()
    cube.scale(0.005, center=cube.get_center())
    cube.translate(
        (
            np.random.uniform(min_vert[0], max_vert[0]),
            np.random.uniform(min_vert[1], max_vert[1]),
            np.random.uniform(min_vert[2], max_vert[2]),
        ),
        relative=False,
    )
    mesh += cube
mesh.compute_vertex_normals()
print("Show input mesh")
o3d.visualization.draw_geometries([mesh])

# Generate data
# Show input mesh

print("Cluster connected triangles")
with o3d.utility.VerbosityContextManager(
        o3d.utility.VerbosityLevel.Debug) as cm:
    triangle_clusters, cluster_n_triangles, cluster_area = (
        mesh.cluster_connected_triangles())
triangle_clusters = np.asarray(triangle_clusters)
cluster_n_triangles = np.asarray(cluster_n_triangles)
cluster_area = np.asarray(cluster_area)

print("Show mesh with small clusters removed")
mesh_0 = copy.deepcopy(mesh)
triangles_to_remove = cluster_n_triangles[triangle_clusters] < 100
mesh_0.remove_triangles_by_mask(triangles_to_remove)
o3d.visualization.draw_geometries([mesh_0])

Show mesh with small clusters removed

print("Show largest cluster")
mesh_1 = copy.deepcopy(mesh)
largest_cluster_idx = cluster_n_triangles.argmax()
triangles_to_remove = triangle_clusters != largest_cluster_idx
mesh_1.remove_triangles_by_mask(triangles_to_remove)
o3d.visualization.draw_geometries([mesh_1])

Show largest cluster