Open3d之网格(Mesh)操作
网格
open3d有一种被称为TriangleMesh的3D三角网格的数据结构。下面的代码展示了如何从一个ply文件读取三角网格数据并打印它的顶点和三角形。
# -*-coding:utf-8 -*-
import copy
import numpy as np
import open3d as o3d
# 加载三角网格
mesh = o3d.io.read_triangle_mesh("/home/ancy/PycharmProjects/learn/Open3D/examples/test_data/knot.ply")
# 打印网格信息
print(mesh)
# 打印网格定点信息
print('Vertices:')
print(np.asarray(mesh.vertices))
# 打印网格的三角形
print('Triangles:')
print(np.asarray(mesh.triangles))
TriangleMesh类是类有一些数据字段,例如vertices和triangles。Open3D通过numpy提供对这些字段的直接内存访问。
可视化3D网格
print("尝试使用法线 (存在: " + str(mesh.has_vertex_normals()) +
") 和颜色 (exist: " + str(mesh.has_vertex_colors()) + ")去渲染网格")
# 因为网格不存在法线和颜色故看起来不美观
o3d.visualization.draw_geometries([mesh])
你可以旋转和移动这个网格,但是由于它是纯灰色的所以看起来不是并不像一个3D图形。这是因为当前网格没有顶点或面的法线,所以使用统一颜色着色而不是用更复杂的Phong着色(冯氏着色)。
表面法线估计
使用曲面法线来绘制网格。
# 计算顶点的法线
mesh.compute_vertex_normals()
# 打印法线信息
print(np.asarray(mesh.triangle_normals))
# 可视化网格
o3d.visualization.draw_geometries([mesh])
裁剪网格
通过numpy直接操作网格的triangle和triangle_normals字段去除一半的曲面。
# 只取前半个三角形网格
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])
网格上色
网格的上色和点云的上色是一致的, 使用paint_uniform_color给网格涂上统一的颜色,颜色在[0,1]范围的RGB空间中。
# 给网格上色
mesh1.paint_uniform_color([1, 0.706, 0])
# 可视化网格
o3d.visualization.draw_geometries([mesh1])
网格属性
三角网格有几个可以用Open3D测试的属性。一个重要的属性是流形性质(manifold property),可使用is_edge_manifold去测试网格是否为边缘流形(edge manifold),使用is_vertex_manifold去测试所有顶点是否为流形。如果每个边缘包括一个或两个三角形,则这个三角网格是边缘流形。is_edge_manifold函数有一个bool型的参数allow_boundary_edges用来指定是否允许边界的边缘。此外,如果顶点的星形边是边缘流形和边缘连接的话,则三角形网格是顶点流形。比如两个或者更多的面可能只有一个顶点连接而不是通过边。
另一个属性是自相交测试。如果在一个网格中存在与另一个网格相交的三角形,is_self_intersecting这个函数就会返回true。一个水密网格能够被定义成一个边缘流形,顶点流形和不自交。在Open3D中通过is_watertight接口实现这种检测。
我们也可测试一个网格是否为可定向的,即三角形可以以所有法线指向外部的方式定向。这个通过is_orientable实现。
下面的代码测试了这些属性并且可视化。非流形边缘用红色表示,边界边缘用绿色标识,非流形顶点用绿色点,自交的三角形用粉色显示。
def edges_to_lineset(mesh, edges, color):
"""将边缘转为线集"""
ls = o3d.geometry.LineSet()
ls.points = mesh.vertices
ls.lines = edges
colors = np.empty((np.asarray(edges).shape[0], 3))
colors[:] = color
ls.colors = o3d.utility.Vector3dVector(colors)
return ls
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(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(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(edges_to_lineset(mesh, edges, (1, 0, 1)))
o3d.visualization.draw_geometries(geoms, mesh_show_back_face=True)
def edges_to_lineset(mesh, edges, color):
"""将边缘转为线集"""
ls = o3d.geometry.LineSet()
ls.points = mesh.vertices
ls.lines = edges
colors = np.empty((np.asarray(edges).shape[0], 3))
colors[:] = color
ls.colors = o3d.utility.Vector3dVector(colors)
return ls
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(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(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(edges_to_lineset(mesh, edges, (1, 0, 1)))
o3d.visualization.draw_geometries(geoms, mesh_show_back_face=True)
def get_non_manifold_edge_mesh():
verts = np.array(
[[-1, 0, 0], [0, 1, 0], [1, 0, 0], [0, -1, 0], [0, 0, 1]],
dtype=np.float64,
)
triangles = np.array([[0, 1, 3], [1, 2, 3], [1, 3, 4]])
mesh = o3d.geometry.TriangleMesh()
mesh.vertices = o3d.utility.Vector3dVector(verts)
mesh.triangles = o3d.utility.Vector3iVector(triangles)
mesh.compute_vertex_normals()
mesh.rotate(
mesh.get_rotation_matrix_from_xyz((np.pi / 4, 0, np.pi / 4)),
center=mesh.get_center(),
)
return mesh
def get_non_manifold_vertex_mesh():
verts = np.array(
[
[-1, 0, -1],
[1, 0, -1],
[0, 1, -1],
[0, 0, 0],
[-1, 0, 1],
[1, 0, 1],
[0, 1, 1],
],
dtype=np.float64,
)
triangles = np.array([
[0, 1, 2],
[0, 1, 3],
[1, 2, 3],
[2, 0, 3],
[4, 5, 6],
[4, 5, 3],
[5, 6, 3],
[4, 6, 3],
])
mesh = o3d.geometry.TriangleMesh()
mesh.vertices = o3d.utility.Vector3dVector(verts)
mesh.triangles = o3d.utility.Vector3iVector(triangles)
mesh.compute_vertex_normals()
mesh.rotate(
mesh.get_rotation_matrix_from_xyz((np.pi / 4, 0, np.pi / 4)),
center=mesh.get_center(),
)
return mesh
def get_open_box_mesh():
mesh = o3d.geometry.TriangleMesh.create_box()
mesh.triangles = o3d.utility.Vector3iVector(np.asarray(mesh.triangles)[:-2])
mesh.compute_vertex_normals()
mesh.rotate(
mesh.get_rotation_matrix_from_xyz((0.8 * np.pi, 0, 0.66 * np.pi)),
center=mesh.get_center(),
)
return mesh
def get_intersecting_boxes_mesh():
mesh0 = o3d.geometry.TriangleMesh.create_box()
T = np.eye(4)
T[:, 3] += (0.5, 0.5, 0.5, 0)
mesh1 = o3d.geometry.TriangleMesh.create_box()
mesh1.transform(T)
mesh = mesh0 + mesh1
mesh.compute_vertex_normals()
mesh.rotate(
mesh.get_rotation_matrix_from_xyz((0.7 * np.pi, 0, 0.6 * np.pi)),
center=mesh.get_center(),
)
return mesh
check_properties('Knot', mesh)
check_properties('Moebius', o3d.geometry.TriangleMesh.create_moebius(twists=1))
check_properties("non-manifold edge", get_non_manifold_edge_mesh())
check_properties("non-manifold vertex", get_non_manifold_vertex_mesh())
check_properties("open box", get_open_box_mesh())
check_properties("intersecting_boxes", get_intersecting_boxes_mesh())
网格滤波
Open3d包含许多网格滤波的算法。接下来我们将会展示相关的一些滤波接口。
均值滤波
最简单的是均值滤波。给定的顶点由相邻顶点的平均值N给出的。公式如下:
可以使用此过滤器对网格进行降噪,如下面的代码所。filter_smooth_simple函数的参数number_of_iterations用来定义应用于网格的滤波器的频率。
# 创建具有噪声的网格
mesh_in = copy.deepcopy(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])
# 使用均值滤波迭代1次
mesh_out = mesh_in.filter_smooth_simple(number_of_iterations=1)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])
# 使用均值滤波迭代5次
mesh_out = mesh_in.filter_smooth_simple(number_of_iterations=5)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])
一次滤波迭代之后
五次迭代后
拉普拉斯算子(Laplacian)
另一个重要的网格过滤器是拉普拉斯算子,其定义如下:
其中λ是滤波器的强度,是与相邻顶点的距离相关的归一化权重. 这个滤波器的接口是filter_smooth_laplacian,有两个参数: number_of_iterations 和 lambda。
# 使用拉普拉斯滤波迭代10次
mesh_out = mesh_in.filter_smooth_laplacian(number_of_iterations=10)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])
# 使用拉普拉斯滤波迭代50次
mesh_out = mesh_in.filter_smooth_laplacian(number_of_iterations=50)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])
使用拉普拉斯滤波迭代10次
使用拉普拉斯滤波迭代50次
Taubin滤波
均值滤波和Laplacian滤波有一个问题是他们会使三角网格收缩。[Taubin1995] 展示了使用两种不同λ参数的Laplacian滤波器来防止网格收缩。这个滤波器的实现接口是:filter_smooth_taubin。
# 使用taubin滤波迭代10次
mesh_out = mesh_in.filter_smooth_taubin(number_of_iterations=10)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])
# 使用taubin滤波迭代100次
mesh_out = mesh_in.filter_smooth_taubin(number_of_iterations=100)
mesh_out.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh_out])
使用taubin滤波迭代10次
使用taubin滤波迭代100次
采样
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])
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])~
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]的附加细分方法。该方法基于四次样条曲线,该样条曲线除了在特殊顶点处生成连续的极限曲面外,其他地方都生成连续的极限曲面。这样可以得到更加平滑的拐角。
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 4409 vertices and 6500 triangles
Simplified mesh has 1982 vertices and 1700 triangles
连通分量
各种重建方法的结果。 Open3D实现了一个连接组件算法cluster_connected_triangles,该算法将每个三角形分配给一组连接的三角形。 它为每个三角形返回cluster_index中的簇索引,每个簇返回cluster_n_triangles中三角形的数目以及cluster_area中簇的表面积。
下面的代码展示cluster_connected_triangles的应用和如何使用它来删除假(spurious)三角形。
# 生成数据
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()
# 显示输入网格
o3d.visualization.draw_geometries([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)
# 删除小的连通分量
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])
# 显示最大的连通分量
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])
更多推荐
所有评论(0)