Point cloud（点云）

译者注：原教程包含大量渲染的可视化效果图，本译文未包含这些图片。请运行文中的代码以查看相应的可视化结果。

本教程演示点云（point cloud）的基本用法。文中的示例代码假设已导入以下模块：

import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt

import open3d_tutorial as o3dtut  # Open3D 官方教程辅助工具

可视化点云

本教程的第一部分读取一个点云并将其可视化。

print("Load a ply point cloud, print it, and render it")
ply_point_cloud = o3d.data.PLYPointCloud()
pcd = o3d.io.read_point_cloud(ply_point_cloud.path)
print(pcd)
print(np.asarray(pcd.points))
o3d.visualization.draw_geometries([pcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])

Load a ply point cloud, print it, and render it
PointCloud with 196133 points.
[[0.65234375 0.84686458 2.37890625]
 [0.65234375 0.83984375 2.38430572]
 [0.66737998 0.83984375 2.37890625]
 ...
 [2.00839925 2.39453125 1.88671875]
 [2.00390625 2.39488506 1.88671875]
 [2.00390625 2.39453125 1.88793314]]
[Open3D WARNING] GLFW Error: Failed to detect any supported platform
[Open3D WARNING] GLFW initialized for headless rendering.

read_point_cloud 从文件中读取一个点云。它会根据文件扩展名尝试解码文件。支持的文件类型列表请参见 File IO。

draw_geometries 用于可视化点云。使用鼠标/触控板可以从不同视角查看几何体。

它看起来像是一个密集的表面，但实际上是以 surfel（面元）形式渲染的点云。该 GUI 支持多种键盘功能。例如，- 键可以减小点（surfel）的大小。

Note

按 H 键可以打印出该 GUI 的完整键盘指令列表。有关可视化 GUI 的更多信息，请参见 Visualization 与 Customized visualization。

Note

在 macOS 上，GUI 窗口可能无法接收键盘事件。这种情况下，请尝试使用 pythonw 而不是 python 来启动 Python。

体素降采样

体素降采样（voxel downsampling）使用一个规则的体素网格（voxel grid），从输入点云中生成一个均匀降采样的点云。它通常作为许多点云处理任务的预处理步骤。该算法分两步执行：

1. 将点归入各个体素（voxel）。
2. 每个被占据的体素通过对其内部所有点取平均，生成且仅生成一个点。

print("Downsample the point cloud with a voxel of 0.05")
downpcd = pcd.voxel_down_sample(voxel_size=0.05)
o3d.visualization.draw_geometries([downpcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])

Downsample the point cloud with a voxel of 0.05
[Open3D WARNING] GLFW initialized for headless rendering.

顶点法向量估计

点云的另一个基本操作是点法向量估计（point normal estimation）。按 N 键可以查看点法向量。- 和 + 键可用于控制法向量的长度。

print("Recompute the normal of the downsampled point cloud")
downpcd.estimate_normals(
    search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30))
o3d.visualization.draw_geometries([downpcd],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024],
                                  point_show_normal=True)

Recompute the normal of the downsampled point cloud
[Open3D WARNING] GLFW initialized for headless rendering.

estimate_normals 为每个点计算法向量。该函数查找相邻点，并使用协方差分析（covariance analysis）计算相邻点的主轴。

该函数接收一个 KDTreeSearchParamHybrid 类的实例作为参数。两个关键参数 radius = 0.1 和 max_nn = 30 分别指定搜索半径和最大最近邻数量。它有 10cm 的搜索半径，并且最多只考虑 30 个邻居以节省计算时间。

Note

协方差分析算法会产生两个相反的方向作为法向量的候选。在不知道几何体全局结构的情况下，这两个方向都可能正确。这就是所谓的法向量朝向问题（normal orientation problem）。如果存在原始法向量，Open3D 会尝试使法向量与原始法向量对齐；否则，Open3D 会进行随机猜测。如果朝向是一个需要关注的问题，则需要进一步调用 orient_normals_to_align_with_direction 和 orient_normals_towards_camera_location 等朝向函数。

访问估计的顶点法向量

估计得到的法向量可以从 downpcd 的 normals 变量中获取。

print("Print a normal vector of the 0th point")
print(downpcd.normals[0])

Print a normal vector of the 0th point
[-0.27566603 -0.89197839 -0.35830543]

要查看其他变量，请使用 help(downpcd)。法向量可以使用 np.asarray 转换为 NumPy 数组。

print("Print the normal vectors of the first 10 points")
print(np.asarray(downpcd.normals)[:10, :])

Print the normal vectors of the first 10 points
[[-0.27566603 -0.89197839 -0.35830543]
 [-0.04230441 -0.99410664 -0.09981149]
 [-0.00399871 -0.99965423 -0.02598917]
 [-0.93768261 -0.07378998  0.3395679 ]
 [-0.43476205 -0.62438493 -0.64894177]
 [-0.09739809 -0.9928602  -0.06886388]
 [-0.27498453 -0.67317361 -0.68645524]
 [-0.11728718 -0.95516445 -0.27185399]
 [-0.00816546 -0.99965616 -0.02491762]
 [-0.11067463 -0.99205156 -0.05987351]]

有关 NumPy 数组的更多示例，请参见 Working with NumPy。

裁剪点云

print("Load a polygon volume and use it to crop the original point cloud")
demo_crop_data = o3d.data.DemoCropPointCloud()
pcd = o3d.io.read_point_cloud(demo_crop_data.point_cloud_path)
vol = o3d.visualization.read_selection_polygon_volume(demo_crop_data.cropped_json_path)
chair = vol.crop_point_cloud(pcd)
o3d.visualization.draw_geometries([chair],
                                  zoom=0.7,
                                  front=[0.5439, -0.2333, -0.8060],
                                  lookat=[2.4615, 2.1331, 1.338],
                                  up=[-0.1781, -0.9708, 0.1608])

Load a polygon volume and use it to crop the original point cloud
[Open3D WARNING] GLFW initialized for headless rendering.

read_selection_polygon_volume 读取一个指定多边形选择区域的 json 文件。vol.crop_point_cloud(pcd) 过滤掉部分点，只保留椅子。

为点云着色

print("Paint chair")
chair.paint_uniform_color([1, 0.706, 0])
o3d.visualization.draw_geometries([chair],
                                  zoom=0.7,
                                  front=[0.5439, -0.2333, -0.8060],
                                  lookat=[2.4615, 2.1331, 1.338],
                                  up=[-0.1781, -0.9708, 0.1608])

Paint chair
[Open3D WARNING] GLFW initialized for headless rendering.

paint_uniform_color 将所有点绘制为统一的颜色。颜色在 RGB 空间中，取值范围为 [0, 1]。

点云距离

Open3D 提供了 compute_point_cloud_distance 方法，用于计算从源点云到目标点云的距离。也就是说，它为源点云中的每个点计算到目标点云中最近点的距离。

在下面的示例中，我们使用该函数计算两个点云之间的差异。需要注意的是，该方法也可用于计算两个点云之间的 Chamfer 距离（倒角距离）。

# Load data
demo_crop_data = o3d.data.DemoCropPointCloud()
pcd = o3d.io.read_point_cloud(demo_crop_data.point_cloud_path)
vol = o3d.visualization.read_selection_polygon_volume(demo_crop_data.cropped_json_path)
chair = vol.crop_point_cloud(pcd)

dists = pcd.compute_point_cloud_distance(chair)
dists = np.asarray(dists)
ind = np.where(dists > 0.01)[0]
pcd_without_chair = pcd.select_by_index(ind)
o3d.visualization.draw_geometries([pcd_without_chair],
                                  zoom=0.3412,
                                  front=[0.4257, -0.2125, -0.8795],
                                  lookat=[2.6172, 2.0475, 1.532],
                                  up=[-0.0694, -0.9768, 0.2024])

[Open3D WARNING] GLFW initialized for headless rendering.

包围盒

PointCloud 几何类型与 Open3D 中的所有其他几何类型一样，具有包围盒（bounding volume）。目前，Open3D 实现了 AxisAlignedBoundingBox（轴对齐包围盒）和 OrientedBoundingBox（有向包围盒），它们也可用于裁剪几何体。

aabb = chair.get_axis_aligned_bounding_box()
aabb.color = (1, 0, 0)
obb = chair.get_oriented_bounding_box()
obb.color = (0, 1, 0)
o3d.visualization.draw_geometries([chair, aabb, obb],
                                  zoom=0.7,
                                  front=[0.5439, -0.2333, -0.8060],
                                  lookat=[2.4615, 2.1331, 1.338],
                                  up=[-0.1781, -0.9708, 0.1608])

[Open3D WARNING] GLFW initialized for headless rendering.

凸包

点云的凸包（convex hull）是包含所有点的最小凸集。Open3D 包含 compute_convex_hull 方法用于计算点云的凸包。该实现基于 Qhull。

在下面的示例代码中，我们首先从一个网格（mesh）采样得到点云，并计算以三角网格形式返回的凸包。然后，我们将凸包可视化为红色的 LineSet（线集）。

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

pcl = mesh.sample_points_poisson_disk(number_of_points=2000)
hull, _ = pcl.compute_convex_hull()
hull_ls = o3d.geometry.LineSet.create_from_triangle_mesh(hull)
hull_ls.paint_uniform_color((1, 0, 0))
o3d.visualization.draw_geometries([pcl, hull_ls])

[Open3D WARNING] GLFW initialized for headless rendering.

DBSCAN 聚类

给定一个来自（例如）深度传感器的点云，我们希望将局部的点云聚簇归并到一起。为此，我们可以使用聚类算法（clustering algorithm）。Open3D 实现了 DBSCAN [Ester1996]，这是一种基于密度的聚类算法。该算法在 cluster_dbscan 中实现，需要两个参数：eps 定义一个聚簇中到邻居的距离，min_points 定义构成一个聚簇所需的最少点数。该函数返回标签（labels），其中标签 -1 表示噪声。

ply_point_cloud = o3d.data.PLYPointCloud()
pcd = o3d.io.read_point_cloud(ply_point_cloud.path)

with o3d.utility.VerbosityContextManager(
        o3d.utility.VerbosityLevel.Debug) as cm:
    labels = np.array(
        pcd.cluster_dbscan(eps=0.02, min_points=10, print_progress=True))

max_label = labels.max()
print(f"point cloud has {max_label+1} clusters")
colors = plt.get_cmap("tab20")(labels / (max_label if max_label > 0 else 1))
colors[labels < 0] = 0
pcd.colors = o3d.utility.Vector3dVector(colors[:, :3])
o3d.visualization.draw_geometries([pcd],
                                  zoom=0.455,
                                  front=[-0.4999, -0.1659, -0.8499],
                                  lookat=[2.1813, 2.0619, 2.0999],
                                  up=[0.1204, -0.9852, 0.1215])

[Open3D DEBUG] Precompute neighbors.
[Open3D DEBUG] Done Precompute neighbors.
[Open3D DEBUG] Compute Clusters
Precompute neighbors.[========================================] 100%
[Open3D DEBUG] Done Compute Clusters: 10
point cloud has 10 clusters
[Open3D WARNING] GLFW initialized for headless rendering.

Note

该算法会预先计算所有点在 epsilon 半径内的邻居。如果所选的 epsilon 过大，这可能需要消耗大量内存。

平面分割

Open3D 还支持使用 RANSAC 从点云中分割几何基元（geometric primitives）。要找到点云中支持度最高的平面，我们可以使用 segment_plane。该方法有三个参数：distance_threshold 定义一个点到估计平面的最大距离（在此距离内仍被视为内点，inlier），ransac_n 定义为估计平面而随机采样的点数，num_iterations 定义随机平面被采样和验证的次数。随后，该函数以 ((a,b,c,d)) 的形式返回平面，使得对于平面上的每个点 ((x,y,z))，都有 (ax + by + cz + d = 0)。该函数还会返回内点索引的列表。

pcd_point_cloud = o3d.data.PCDPointCloud()
pcd = o3d.io.read_point_cloud(pcd_point_cloud.path)

plane_model, inliers = pcd.segment_plane(distance_threshold=0.01,
                                         ransac_n=3,
                                         num_iterations=1000)
[a, b, c, d] = plane_model
print(f"Plane equation: {a:.2f}x + {b:.2f}y + {c:.2f}z + {d:.2f} = 0")

inlier_cloud = pcd.select_by_index(inliers)
inlier_cloud.paint_uniform_color([1.0, 0, 0])
outlier_cloud = pcd.select_by_index(inliers, invert=True)
o3d.visualization.draw_geometries([inlier_cloud, outlier_cloud],
                                  zoom=0.8,
                                  front=[-0.4999, -0.1659, -0.8499],
                                  lookat=[2.1813, 2.0619, 2.0999],
                                  up=[0.1204, -0.9852, 0.1215])

Plane equation: -0.06x + -0.10y + 0.99z + -1.06 = 0
[Open3D WARNING] GLFW initialized for headless rendering.

平面块检测

除了查找支持度最高的单个平面之外，Open3D 还包含一种基于鲁棒统计方法的平面块检测（planar patch detection）算法 [ArujoAndOliveira2020]。该算法首先将点云划分为更小的块（使用八叉树，octree），然后尝试对每个块拟合一个平面。如果平面通过了鲁棒平面性检验，则被接受。通过对点子集取中位数点位置和中位数点法向量并估计平面 (ax + by + cz + d = 0)，即可拟合出一个平面。鲁棒平面性检验由两个主要部分组成。首先，求出每个点法向量与拟合平面法向量之间夹角的分布。如果该分布的离散程度过大（即所有关联点法向量之间的方差过大），则拒绝该平面。其次，计算从拟合平面到每个点的距离分布。如果该分布的离散程度过大（使用共面性度量，coplanarity metric，参见 [ArujoAndOliveira2020] 的图 4），则拒绝该平面。在找到初始的一组平面之后，会使用一个迭代过程来生长并合并平面，最终得到一个更小、更稳定的平面集合。随后，可以使用这些平面关联点集的二维凸包对它们进行边界框定，从而提取出平面块（planar patch）。

要查找点云中的平面块列表，我们可以使用 detect_planar_patches。该方法可接受六个参数：normal_variance_threshold_deg 控制点法向量之间允许的方差大小，默认值为 (60^\circ)。较小的值通常会产生更少但质量更高的平面。coplanarity_deg 控制点到平面距离的允许分布，默认值为 (75^\circ)。较大的值会促使点更紧密地分布在拟合平面周围。outlier_ratio 设定在拒绝之前，拟合平面关联点集中允许的最大离群点比例，默认值为 0.75。min_plane_edge_length 用于拒绝误报——一个平面块的最大边必须大于该值才会被视为真正的平面块。如果保持为 0，算法默认使用点云最大维度的 1%。min_num_points 决定关联八叉树的深度以及尝试拟合平面时必须存在的点数。如果保持为 0，算法默认使用点云中点数的 0.1%。search_param 是 geometry::KDTreeSearchParam 的一个实例，默认为 geometry::KDTreeSearchParamKNN。在生长和合并平面时会用到每个点的 k 个最近邻。较大的 k 值通常会产生质量更高的平面块，但计算开销也更大。该函数随后返回检测到的平面块列表，以 geometry::OrientedBoundingBox 对象的形式表示，其中 R 的第三列（即 (z)）指示平面块的法向量。(z) 方向上的范围非零，以便 OrientedBoundingBox 能够包含参与平面检测的点。可以使用 geometry::TriangleMesh::CreateFromOrientedBoundingBox 工厂函数对平面块进行可视化，并使用 scale 参数沿法向量”压平”包围盒，即 CreateFromOrientedBoundingBox(obox, scale=[1, 1, 0.0001])。

dataset = o3d.data.PCDPointCloud()
pcd = o3d.io.read_point_cloud(dataset.path)
assert (pcd.has_normals())

# using all defaults
oboxes = pcd.detect_planar_patches(
    normal_variance_threshold_deg=60,
    coplanarity_deg=75,
    outlier_ratio=0.75,
    min_plane_edge_length=0,
    min_num_points=0,
    search_param=o3d.geometry.KDTreeSearchParamKNN(knn=30))

print("Detected {} patches".format(len(oboxes)))

geometries = []
for obox in oboxes:
    mesh = o3d.geometry.TriangleMesh.create_from_oriented_bounding_box(obox, scale=[1, 1, 0.0001])
    mesh.paint_uniform_color(obox.color)
    geometries.append(mesh)
    geometries.append(obox)
geometries.append(pcd)

o3d.visualization.draw_geometries(geometries,
                                  zoom=0.62,
                                  front=[0.4361, -0.2632, -0.8605],
                                  lookat=[2.4947, 1.7728, 1.5541],
                                  up=[-0.1726, -0.9630, 0.2071])

Detected 10 patches
[Open3D WARNING] GLFW initialized for headless rendering.

隐点去除

假设你想从给定视角渲染一个点云，但背景中的点因为没有被其他点遮挡而泄露到前景中。为此，我们可以应用隐点去除（hidden point removal）算法。Open3D 实现了 [Katz2007] 的方法，它可以在不进行表面重建或法向量估计的情况下，近似地求出从给定视角看点云的可见性。

print("Convert mesh to a point cloud and estimate dimensions")
armadillo = o3d.data.ArmadilloMesh()
mesh = o3d.io.read_triangle_mesh(armadillo.path)
mesh.compute_vertex_normals()

pcd = mesh.sample_points_poisson_disk(5000)
diameter = np.linalg.norm(
    np.asarray(pcd.get_max_bound()) - np.asarray(pcd.get_min_bound()))
o3d.visualization.draw_geometries([pcd])

Convert mesh to a point cloud and estimate dimensions
[Open3D WARNING] GLFW initialized for headless rendering.

print("Define parameters used for hidden_point_removal")
camera = [0, 0, diameter]
radius = diameter * 100

print("Get all points that are visible from given view point")
_, pt_map = pcd.hidden_point_removal(camera, radius)

print("Visualize result")
pcd = pcd.select_by_index(pt_map)
o3d.visualization.draw_geometries([pcd])

Define parameters used for hidden_point_removal
Get all points that are visible from given view point
Visualize result
[Open3D WARNING] GLFW initialized for headless rendering.