02-01-02.Mesh

Mesh

---

Open3D有一个用于3D三角形网格的数据结构，称为TriangleMesh。下面的代码展示了如何从一个ply文件中读取一个三角形网格并打印它的顶点和三角形。

print("Testing mesh in Open3D...")
mesh = o3dtut.get_knot_mesh()
print(mesh)
print('Vertices:')
print(np.asarray(mesh.vertices))
print('Triangles:')
print(np.asarray(mesh.triangles))

结果：

Testing mesh in Open3D...
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 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: True) and colors (exist: False)
A mesh with no normals and no colors does not look good.

你可以旋转和移动网格，但它被涂上了统一的灰色，看起来不像3d。原因是当前的网格没有顶点或面的法线。所以用统一颜色的阴影代替更复杂的Phong阴影。

表面法向量估计

---

让我们用表面法线绘制网格。

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的成员函数。

裁剪mesh

---

我们通过直接操作网格的triangle和triangle数据字段来移除一半的表面。这是通过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.

给mesh上色

---

pait_uniform_color用统一的颜色绘制网格。颜色在RGB空间$[0,1]$范围内。

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

结果：

Painting the mesh

mesh 属性

---

一个三角形网格有几个可以用Open3D测试的属性。一个重要的性质是流形性质，我们可以测试三角形网格是is_edge_manifold还是is_vertex_manifold。如果每条边包围一个或两个三角形，那么三角形网格是边流形。函数is_edge_manifold具有bool参数allow_boundary_edges，该参数定义是否应允许边界边缘。此外，如果三角形网格顶点的星形是边形且边连通，则三角形网格是顶点流形，例如两个或多个面仅由顶点a连接，那么三角形网格就是顶点流形。
另一个性质是自交检验。如果在一个网格中有一个三角形与另一个网格相交，函数is_self_intersecting 返回 True。水密网格可以定义为边流形、顶点流形且不自交的网格。函数is_watertight在Open3D中实现了这个检查。
我们也可以测试三角形网格，如果它是可定向的，即三角形可以以这样一种方式定向，所有法线指向外面。相应的函数在Open3D中被称为is_orientable。
下面的代码根据这些属性测试一些三角形网格，并将结果可视化。非流形边用红色表示，边界边用绿色表示，非流形顶点用绿色点表示，自相交三角形用粉色表示。

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)

check_properties('Knot', o3dtut.get_knot_mesh())
check_properties('Moebius', o3d.geometry.TriangleMesh.create_moebius(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())

结果：

Knot
  edge_manifold:          True
  edge_manifold_boundary: True
  vertex_manifold:        True
  self_intersecting:      False
  watertight:             True
  orientable:             True

Moebius
  edge_manifold:          True
  edge_manifold_boundary: False
  vertex_manifold:        True
  self_intersecting:      False
  watertight:             False
  orientable:             False

non-manifold edge
  edge_manifold:          False
  edge_manifold_boundary: False
  vertex_manifold:        True
  self_intersecting:      False
  watertight:             False
  orientable:             True

non-manifold vertex
  edge_manifold:          True
  edge_manifold_boundary: True
  vertex_manifold:        False
  self_intersecting:      False
  watertight:             False
  orientable:             True

open box
  edge_manifold:          True
  edge_manifold_boundary: False
  vertex_manifold:        True
  self_intersecting:      False
  watertight:             False
  orientable:             True

intersecting_boxes
  edge_manifold:          True
  edge_manifold_boundary: True
  vertex_manifold:        True
  self_intersecting:      True
  watertight:             False
  orientable:             True
  # visualize self-intersecting triangles

Mesh filtering

---

Open3D包含了许多过滤网格的方法。下面我们展示了实现的滤波器来平滑有噪声的三角形网格。

Average filter

最简单的过滤器是平均过滤器。给定顶点$v_i$是由相邻顶点$N$的平均值给出的。

$$v_i=\frac{v_i+{\sum }_{n\in \mathcal{N}}{v}_n}{|N|+1}$$

如下面的代码所示，这个过滤器可以用来消除网格的噪声。filter_smooth_simple函数中迭代次数的参数number_of_iterations定义了过滤器应用到网格的频率。

print('create noisy mesh')
mesh_in = o3dtut.get_knot_mesh()
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

filter with average with 5 iterations

Laplacian

另一个重要的滤镜是拉普拉斯过滤：

$$v_i=v_i\cdot \lambda \sum _{n\in N}{w}_n{v}_n-v_i,$$

其中$\lambda$为滤波器的强度，$w_n$为归一化权值，与相邻顶点的距离有关。该滤波器采用filter_smooth_laplacian实现，具有num_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])

结果：

filter with Laplacian with 10 iterations

filter with Laplacian with 50 iterations

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

结果：

filter with Taubin with 10 iterations

filter with Taubin with 100 iterations

Sampling

Open3D包含了从三角形网格中采样点云的函数。最简单的方法是sample_points_uniformly基于三角形区域从三维曲面上均匀采样点。点的数量参数number_of_points定义了从三角形表面采样的点的数量。

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

mesh = o3dtut.get_bunny_mesh()
mesh.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh])
pcd = mesh.sample_points_uniformly(number_of_points=500)
o3d.visualization.draw_geometries([pcd])

均匀采样可以在表面上产生点簇，而一种称为泊松盘采样的方法可以在表面上均匀分布点。该方法sample_points_poisson_dick实现了样本消除。它从一个采样点云开始，然后移除满足采样准则的点。该方法支持两个选项来提供初始点云：

1. 通过init_factor参数默认:该方法首先使用init_factor x num_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 = o3dtut.get_bunny_mesh()
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])

网格细分

在网格细分中，我们将每个三角形分成若干个较小的三角形。在最简单的情况下，我们计算每个三角形每边的中点，并将三角形分成四个较小的三角形。这是在细分中点函数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$连续的极限曲面，除了在特殊顶点处，它们是$C^1$连续的。这将导致更平滑的角。

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 = o3dtut.get_knot_mesh()
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=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 1440 vertices and 2880 triangles

After subdivision it has 5760 vertices and 11520 triangles

网格简化

---

有时我们想用较少的三角形和顶点来表示一个高分辨率的网格，但是低分辨率的网格仍然应该接近高分辨率的网格。为此，Open3D实现了许多网格简化方法。

顶点聚类

顶点聚类方法将所有属于给定大小体素的顶点集中到一个顶点。该方法是在simplify_vertex_clustering中实现的，作为参数的voxel_size定义了体素网格的大小和contraction定义了顶点池的大小。o3d.geometry.SimplificationContraction.Average计算一个简单的平均数。

mesh_in = o3dtut.get_bunny_mesh()
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 35947 vertices and 69451 triangles

voxel_size = 4.865594e-03
Simplified mesh has 3222 vertices and 6454 triangles

voxel_size = 9.731187e-03
Simplified mesh has 845 vertices and 1724 triangles

网格简化

另一类网格简化方法是网格抽取，以增量步骤操作。我们选择一个最小误差的三角形并删除它。这样重复，直到达到要求的三角形数量。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])

结果：

Simplified mesh has 4405 vertices and 6499 triangles

Simplified mesh has 1978 vertices and 1700 triangles。

连通区域

各种重建方法的结果。Open3D实现了一个连接组件算法cluster_connected_triangles，将每个三角形分配给一组连接三角形。它为每个三角形返回三角形簇中的簇索引，每个集群cluster_n_triangles中的三角形数量和cluster_area中集群的表面积。
例如，这在RGBD Integrarion中很有用，RGBD集成并不总是单个三角形网格，而是多个网格。 一些较小的零件是由于噪音引起的，我们很可能希望将其移除。
下面的代码显示了cluster_connected_triangles的应用以及如何使用它来消除虚假三角形。

print("Generate data")
mesh = o3dtut.get_bunny_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)

结果为：

Cluster connected triangles
[Open3D DEBUG] [ClusterConnectedTriangles] Compute triangle adjacency
[Open3D DEBUG] [ClusterConnectedTriangles] Done computing triangle adjacency
[Open3D DEBUG] [ClusterConnectedTriangles] Done clustering, #clusters=31

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