Open3d点云对象详解

PointCloud是open3d中用于点云处理的类，封装了包括几何变换、数据滤波、聚类分割等一系列实用算法。如无特别说明，本例中所有例程均基于斯坦福兔子的点云模型，下载地址：斯坦福标准模型

# 此行代码后面不再重复引入
import open3d as o3d
# 载入斯坦福兔子 rabbit.pcd文件需在当前python工作的文件夹中
pcd = o3d.io.read_point_cloud("rabbit.pcd")	

读取和清除点云

一般点云数据的读取方法属于open3d.io的内容，但点云类也提供了一些生成点云的方法，最简单的是创建一个点云对象后，为其点集进行赋值，下例中xyz为Nx3的矩阵

import open3d as o3d
pcd = o3d.geometry.PointCloud()
pcd.points = o3d.utility.Vector3dVector(xyz)

如果想让pcd恢复到没有点云的状态，可以通过pcd.clear()来清除。

此外，点云类还提供了两个静态方法，create_from_depth_image和create_from_rgbd_image，使得用户可以通过深度图像或rgbd图像来生成点云。

由于手头没有深度图像，所以暂且用常见的RGB图像来演示一下这个功能，用下面这张图像生成深度图像，并绘制点云

import open3d as o3d

raw = o3d.io.read_image("test.jpg")
rgbd = o3d.geometry.RGBDImage.create_from_color_and_depth(raw, raw, convert_rgb_to_intensity=False)

inter = o3d.camera.PinholeCameraIntrinsic()
# 此为相机内部参数，是我乱写的 1920x1080为像素尺寸；(600,600)为x,y方向的焦距；(640,360)为相机位置
inter.set_intrinsics(1920, 1080, 600, 600, 640, 360)
pcd = o3d.geometry.PointCloud().create_from_rgbd_image(rgbd, inter)
o3d.visualization.draw_geometries([pcd])

绘图结果为

感觉很有冲击力。

点云的颜色可以通过paint_uniform_color来指定，后面会经常用到这个方法。

点云属性

点云类中封装了一些方法，用以提取点云的特征

对象属性

>>> pcd.dimension()
3           # 点云维度为3
>>> pcd.has_colors()
False       # 点云无颜色
>>> pcd.has_covariances()
False       # 对象中不含协方差矩阵
>>> pcd.has_normals()
False       # 不含法向量
>>> pcd.has_points()
True        # 点云中有点
>>> pcd.is_empty()
False       # 点云不为空
>>> pcd.get_geometry_type()
<Type.PointCloud: 1>    #几何体的类型是点云

点云特征

通过下面三个函数，可以获取点云对象$x,y,z$三个坐标的最大值、最小值以及质心。

>>> pcd.get_max_bound()
array([8.7769  , 9.210494, 4.985259])
>>> pcd.get_min_bound()
array([-6.793   , -6.222866, -7.082071])
>>> pcd.get_center()
array([-1.00036164e-06, -1.81378116e-08,  4.44904992e-07])

点云框线

点云边框

open3d为点云对象提供两种生成边框的方法，分别是生成定向边框的get_oriented_bounding_box和生成轴对齐边框的get_axis_aligned_bounding_box。

obb = pcd.get_oriented_bounding_box()
obb.color = [1,0,0]
aabb = pcd.get_axis_aligned_bounding_box()
aabb.color = [0,1,0]

效果为

其中红色为定向边框，相当于最高效地将斯坦福兔子圈在了框中；绿色边框则横平竖直，每条框线均与坐标轴平行。

凸包

凸包的本质也是一种点云框线，但其采取的框选方案要更加紧密，即选取点云中最外侧的点连接成一个凸多面体。具体情况看下图就能明白

hull, _ = pcd.compute_convex_hull()
hull_ls = o3d.geometry.LineSet.create_from_triangle_mesh(hull)
o3d.visualization.draw_geometries([pcd, hull_ls])

几何变换

点云类提供了平移、缩放以及旋转的空间变换方法，分别用到函数translate、scale以及rotate，示例如下：

import copy

# 此为原始点云
pcd0 = o3d.io.read_point_cloud("rabbit.pcd")

# 1号点云，沿x轴移动20米
pcd1 = copy.deepcopy(pcd0).translate((20, 0, 0))
pcd1.paint_uniform_color([1, 0, 0])  #设为红色，100 表示rgb

# 2号点云，演示缩放
pcd2 = copy.deepcopy(pcd0).translate((40, 0, 0))
# 前一个参数为缩放比例；后一个参数为缩放后的位置，为必填项
pcd2.scale(0.5, center=pcd2.get_center())
pcd2.paint_uniform_color([0, 1, 1])

# 3号点云，通过欧拉角旋转
pcd3 = copy.deepcopy(pcd0).translate((0, -20, 0))
# 采用Euler角的方法生成旋转角，表示绕y轴旋转90°
thEuler = pcd3.get_rotation_matrix_from_xyz((0,np.pi/2,0))
pcd3.rotate(thEuler)
pcd3.paint_uniform_color([0, 1, 0])

# 4号点云，通过轴角法旋转
pcd4 = copy.deepcopy(pcd0).translate((20, -20, 0))
# 通过轴角表示法生成旋转角，表示绕y轴旋转60°
th = np.array([0, np.pi/3, 0]).T
thAxis = pcd4.get_rotation_matrix_from_axis_angle(th)
pcd4.rotate(thAxis)
pcd4.paint_uniform_color([0, 0, 1])

# 5号点云，通过四元数旋转
pcd5 = copy.deepcopy(pcd0).translate((40, -20, 0))
# 通过四元数法生成转角，表示绕x轴旋转180°
quart = np.array([0,0,0,1]).T
thQuart = pcd5.get_rotation_matrix_from_quaternion(quart)
pcd5.rotate(thQuart)
pcd5.paint_uniform_color([1, 0, 1])

# 5绘图函数可以输入点云列表
pcds = [pcd0, pcd1, pcd2, pcd3, pcd4, pcd5]
o3d.visualization.draw_geometries(pcds)

效果如下

其中平移和缩放没什么好说的，对于旋转，需要知道欧拉角有三个分量，分别是横滚、俯仰以及航向，代表绕x、y、z轴旋转。

考虑一架正在飞行的飞机，以某一时刻前后为x轴、左右为y轴、上下为z轴。则航向角就是飞机前进方向偏离的角度；俯仰角就是飞机头尾姿态的俯仰；横滚角描述的就是飞机翅膀的摆动。

欧拉角用于描述静态角度是没问题的，但用于表示旋转时会导致万向节死锁，简单来说就是飞机的航向角变化90°之后，其横滚轴变成了俯仰轴，而四元数则没有这个顾虑。

设欧拉角为$(\psi ,\theta ,\varphi )$，则四元数可表示为

$$\mathbf{Q}=\left[\begin{array}{c}\cos (\psi /2)\\ 0\\ 0\\ \sin (\psi /2)\end{array}\right]\left[\begin{array}{cccc}\cos (\theta /2)& 0& \sin (\theta /2)& 0\end{array}\right]\left[\begin{array}{c}\cos (\varphi /2)\\ \sin (\varphi /2)\\ 0\\ 0\end{array}\right]$$

点云类中的生成旋转矩阵的方法均为静态方法，可在不建立对象的情况下调用。

在演示变换的过程中调用了三个，其前缀均为get_rotation_matrix_from_，结尾是axis_angle表示通过欧拉角生成旋转矩阵；quaternion通过四元数；from_xyz则通过旋转向量。

open3d支持通过不同顺序的xyz数组创建旋转矩阵，由于三个坐标总计有6个组合，故而提供了6个静态方法。例如，针对向量$(x,y,z)$，其创建旋转矩阵的方法为get_rotation_matrix_from_xyz；对于$(y,x,z)$，只需将xyz换为yxz即可，非常方便记忆。

法线

对于不包含法线的点云，通过estimate_normals方法可以生成法线。

>>> pcd.estimate_normals()
>>> pcd.has_normals()
True    # 可见pcd已经生成了法线
>>> o3d.visualization.draw_geometries([pcd], point_show_normal=True)

结果表明，把法线画出来实在是太丑了，可谓san值狂掉，属于精神污染了。

estimate_normals中有两个输入参数，其中search_param为k-d树的索引参数，用以生成法线时索引邻近点；fast_normal_computation为布尔型的参数，顾名思义，为True时将开启加速，但可能会导致数值不稳定。

点云类中对法线的操作主要是两种，一是归一化；二是定向。

normalize_normals可对法线进行归一化，并返回归一化法线后的点云。

定向则有三个方法：

• orient_normals_consistent_tangent_plane：通过切平面定向
• orient_normals_to_align_with_direction：通过对其坐标轴定向
• orient_normals_towards_camera_location：通过相机位置定向

特征计算

协方差矩阵

对于点云中的某点$r_i$，其协方差矩阵为

$$C_i=\frac{1}{N}\sum _{i=1}^N(r_i-\bar{r})\cdot (r_i-\bar{r}{)}^T$$

点云对象通过estimate_covariances函数计算每个点的协方差矩阵。此外，compute_mean_and_covariance可计算整个点云及其协方差的均值。

>>> M_and_C = pcd.compute_mean_and_covariance()
>>> pprint(M_and_C)
(array([-1.00036164e-06, -1.81378116e-08,  4.44904992e-07]),
 array([[16.80007645, -5.77246083,  0.46407276],
       [-5.77246083, 17.24825378, -2.5936246 ],
       [ 0.46407276, -2.5936246 ,  7.9322634 ]]))

最邻近距离和马氏距离

点云类为单点提供了两种距离属性，分别是最邻近距离和Mahalanobis距离，对应的函数分别是compute_nearest_neighbor_distance和compute_mahalanobis_distance。

最邻近距离即当前点和离它最近的其他点的距离；Mahalanobis距离的定义则为

$$D_M(r_i)=\sqrt{(r_i-\bar{r}{)}^T{C}_i(r_i-\bar{r})}$$

其中$C_i$为该点对应的协方差矩阵。

索引、采样和滤波

点云索引

select_by_index(ids, invert)方法，可以通过点号对点云进行选取。其中ids为点号列表，invert参数默认为False，表示返回ids中的点；若将invert置为True，将返回ids之外的点。

import open3d as o3d
import numpy as np
pcd = o3d.io.read_point_cloud("rabbit.pcd")
idx = np.arange(10000)
# 索引对应的点
pIn = pcd.select_by_index(idx)
pIn.paint_uniform_color([1, 0, 0])
# 索引外的点云
pOut = pcd.select_by_index(idx, invert=True)
pOut.paint_uniform_color([0, 1, 0])
o3d.visualization.draw_geometries([pIn, pOut])

效果为

无效点剔除

remove_non_finite_points方法用于剔除NaN或者Inf之类的无穷值。其两个bool型的输入参数remove_nan和remove_infinite分别用于指明这两者。

统计滤波和邻域滤波

这两中滤波方法都是先得到符合要求的点索引，然后通过索引滤波，将这些点挑选出来，输出输出为滤波后的点云和点的索引号。

# 上接索引滤波的内容

pcd1 = copy.deepcopy(pcd).translate((20, 0, 0))
pcd2 = copy.deepcopy(pcd).translate((40, 0, 0))

# 统计滤波，参数分别表示K邻域点的个数和标准差乘数
sPcd, sInd = pcd1.remove_statistical_outlier(6, 2.0)

# 半径滤波，输入参数为邻域球内最少点数和邻域半径
rPcd, rInd = pcd2.remove_radius_outlier(9, 0.   )

o3d.visualization.draw_geometries([sPcd, rPcd])

效果如下

这两种算法的逻辑是一样的，对于某点$x$，选取距离$x$最近的一些点，如果这些点的标准差小于设定值，则符合统计滤波的标准；如果均小于邻域半径，则符合半径滤波的标准。

下采样

点云对象共有三种下采样方案

	API	输入参数
随机下采样	random_down_sample	采样率
等序下采样	uniform_down_sample	采样间隔
体素下采样	voxel_down_sample	体素尺寸

其中随机下采样即根据采样率，随机选择一些点；等序则根据输入的采样间隔，则取等间隔的点的序号。

体素采样稍显复杂，会构建三维体素格网，然后输出格网内的点云质心，而非原始数据；若存在法线或颜色，则通通取均值。

此外，还有voxel_down_sample_and_trace函数，在体素采样的基础上，可以规定体素的边界范围。

downpcd = pcd.voxel_down_sample(20)

聚类算法

DBSCAN聚类

DBSCAN，即Density-Based Spatial Clustering of Applications with Noise，基于密度的噪声应用空间聚类。

在DBSCAN算法中，将数据点分为三类：

1. 核心点：若样本$x_i$的$\epsilon$邻域内至少包含了$M$个点，则为核心点
2. 边界点：若样本$x_i$的$\epsilon$邻域内包含的点数小于$M$，但在其他核心点的$\epsilon$邻域内，则为边界点
3. 噪声：既非核心点也非边界点则为噪声

可见，DBSCAN算法需要两个参数，分别是邻域半径$\epsilon$和点数$M$。

在open3d中，提供了cluster_dbscan接口，示例如下

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

pcd = o3d.io.read_point_cloud("rabbit.pcd")

eps = 0.5       # 同一聚类中最大点间距
M = 50          # 有效聚类的最小点数
Labels = np.array(pcd.cluster_dbscan(eps, M))
print(np.max(Labels))                   # 得到结果为3
cs = plt.get_cmap("jet")(Labels/3)      # 伪彩映射
cs[labels < 0] = 0  # labels = -1 的簇为噪声，以黑色显示
pcd.colors = o3d.utility.Vector3dVector(cs[:, :3])
o3d.visualization.draw_geometries([pcd])

结果如图所示

由于斯坦福兔子的扫描点比较连续，所以分割效果不太还，为了更好地演示DBSCAN算法，可以用手骨模型，下载地址：斯坦福标准模型

RANSAC 平面分割

RANSAC，即RANdom SAmple Consensus，随机抽样一致算法。

以平面上的点集举例，假设点集中有一条直线$L$，$L$外的点很少，均为噪声。

那么第一步，随机选取两个点连成一条直线$\hat{L}$，那么这条直线有可能就是$L$，也有可能是噪声连出来的莫名其妙的一条线。

接下来，随机抽取点集中的一些点，如果随机抽取的大部分点都落在$L$附近，那么就说明$\hat{L}$有很大的概率就是$L$；否则说明不太像是$L$。随着抽取出的直线越来越多，最后可以得到最接近$L$的直线，从而完成了对点集的分割。

在Open3d中，提供了基于RANSAC算法的平面分割接口segment_plane

pcd = o3d.io.read_point_cloud("rabbit.pcd")

d = 0.2             # 内点到平面模型的最大距离
n = 5               # 用于拟合平面的采样点数
nIter = 50          # 最大迭代次数

# 返回模型系数plane和内点索引ids，并赋值
plane, ids = pcd.segment_plane(d, n, nIter)

# 平面方程
[a, b, c, d] = plane

# 平面内点点云
iCloud = pcd.select_by_index(ids)
iCloud.paint_uniform_color([0, 0, 1.0])

# 平面外点点云
oCloud = pcd.select_by_index(ids, invert=True)
oCloud.paint_uniform_color([1.0, 0, 0])

# 可视化平面分割结果
o3d.visualization.draw_geometries([iCloud, oCloud])

最后得到的结果为

本来以为平面会出现在兔子的底座上，没想到会把兔子一分为二。