Open3d学习计划（4)网格

点云PCL博主

发布于 2020-06-28 14:41:23

3.3K1

发布于 2020-06-28 14:41:23

文章被收录于专栏：点云PCL点云PCL

Open3D是一个开源库，支持快速开发和处理3D数据。Open3D在c++和Python中公开了一组精心选择的数据结构和算法。后端是高度优化的，并且是为并行化而设置的。

本系列学习计划有Blue同学作为发起人，主要以Open3D官方网站的教程为主进行翻译与实践的学习计划。点云PCL公众号作为免费的3D视觉，点云交流社区，期待有使用Open3D或者感兴趣的小伙伴能够加入我们的翻译计划，贡献免费交流社区，为使用Open3D提供中文的使用教程。

三角网格

open3d有一种被称为TriangleMesh的3d三角网格的数据结构。下面的代码展示了如何从一个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 ...
>>>>geometry::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直接对这些字段进行访问。

可视化网格

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)

这个网格没有法线和颜色看起来不是那么美观。你可以旋转和移动这个网格，但是由于它是纯灰色的所以看起来不是那么像一个’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这两个网格计算函数。

裁剪网格

通过直接操作网格的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.

网格上色

网格的上色和点云的上色是一致的。

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


>>>>Painting the mesh

网格属性

三角网格有几个可以用open3d测试的属性。一个重要的属性是流形性质(manifold property)，可以使用is_edge_manifold去测试网格是不是边缘流形（edge manifold）和is_vertex_manifold去测试是否所有顶点为流形。如果一个三角网格是边缘流形，并且每个边缘包括一个或两个三角形。is_edge_manifold函数有一个bool型的参数allow_boundary_edges用来指定是否允许边界的边缘(这里是API中关于这个参数的解释：If true, than non-manifold edges are defined as edges with more than two adjacent triangles, otherwise each edge that is not adjacent to two triangles is defined as non-manifold.感觉这不好翻译所以在这里放出原文)。此外，如果定点的星形边（star of the vertex）是边缘流形（edge-manifold）和边缘连接（edge-connected）的话，三角形网格是顶点流形。比如两个或者更多的面可能只有一个顶点连接而不是通过边。另一个属性是自交测试。如果在一个网格中存在以另一个网格相交的三角形，is_self_intersecting这个函数就会返回true。一个水密网格(watertight mesh)能够被定义成一个边缘流形edge-manifold），顶点流形（vertex manifold）和不自交（not self intersecting）。在open3d中通过is_watertight接口实现这种检测。我们也可以测试一个网格是否是可定向的(orientable)，也就是说三角形可以通过所有法线指向外面的方式定向。这个通过is_orientable实现。下面的代码测试了这些属性并且可视化。非流形边缘用红色表示，边界边缘（boundary edge）用绿色标识，非流形顶点用绿色点，自交的三角形用粉色显示。

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

>>>>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:             Tru

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

网格滤波

Open3d包含许多网格滤波的算法。接下来我们将会展示相关的一些滤波接口。

均值滤波

最简单的是均值滤波。一个顶点 vi\ v_{i} vi 的值是通过相邻顶点的平均给出的。公式如下：

如下面代码所示，该滤波器能用以网格去噪。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

五次迭代后

拉普拉斯算子（Laplacian）

另一个重要的网格过滤器是拉普拉斯算子，其定义如下：

这里的 λ\lambdaλ 是滤波器强度， wn\ w_{n} wn 是与邻点距离相关的归一化权重。这个滤波器的接口是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])


>>>>filter with Laplacian with 10 iterations

Taubin滤波

均值滤波和Laplacian滤波有一个问题是他们会使三角网格收缩。[Taubin1995] 展示了使用两种不同λ\lambdaλ 参数的Laplacian滤波器来防止网格收缩。这个滤波器的实现接口是：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])

采样

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_disk实现了该功能。因此该算法从一个采样后的点云开始，移除点以满足采样标准。这个算法支持两个初始点云的选择方法。

默认通过参数init_factor：首先通过init_factor x number_of_points来从网格中均匀采样点云，之后进行消除。
可以直接提供一个点云数据给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])~

网格细分

网格细分就是把每个三角形划分为更小的三角形。最简单的方式就是，计算三角形每个边的中点，将其划分为四个较小的三角形。这个通过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)

Open3d实现了基于[Loop1987]的附加细分方法。该方法基于四次箱样条，该四边形样条线会在除 C2连续的非凡顶点之外的所有地方生成 C1连续的极限曲面。（The method is based on a quartic box spline, which generate C2continuous limit surfaces everywhere except at extraordinary vertices where they are C1 continuous.）这样可以得到更加平滑的拐角。

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)

网格简化

有时候我们想用较少的三角形来表示一个高分辨率的网格，但是低分辨率的网格仍然应该接近高分辨率的网格。为此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])~

网格抽取（mesh decimation）

网格细分的另一种方式是逐步执行的网格抽取。我们选取一个最小化误差指标的三角形移除。重复此过程直到满足指定的三角形数量时停止。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])

连通分量（connected components）

各种重建算法的结果（比如RGBD Integration 并不是只有一个三角网格而是有多个网格。一些较小的部分（相比如主物体较小）是由于噪声引起的，我们会想要移除它。Open3d实现了一个连通分量算法cluster_connected_triangles，他将每个三角形分配给一个连接的三角集群。他会从集群中返回每一个三角形的索引triangle_cluters，和每一个集群中三角形的数量cluter_n_triangles还有集群的表面积cluster_area。下面的代码展示cluster_connected_triangles的应用和如何使用它来删除假（spurious）三角形。

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