爆肝5万字 Open3D 点云数据处理基础（Python版））-轻识

#爆肝5万字❤️Open3D 点云数据处理基础（Python版）

历时一个月，Open3D 点云数据处理基础（Python版） 终于初具雏形！目前主要是一些基础内容，后续会继续完善，希望互相学习！

1 概述

Open3D是一个开源库，支持快速开发处理3D数据的软件。Open3D后端是用C++实现的，经过高度优化并通过Python的前端接口公开。Open3D提供了三种数据结构：点云（point cloud）、网格（mesh）和RGB-D图像。对于每个表示，open3D都实现了一整套基本处理算法，如I/O、采样、可视化和数据转换。此外，还包括一些常用的算法，如法线估计、ICP配准等。

open3D包含了九个模块，如下表所示：

模块	功能
Geometry 几何模块	数据结构和基本处理算法
Camera 相机模块	相机模型和相机轨迹
Odometry 里程计模块	RGB-D图像的跟踪与对齐
Registration 配准模块	全局和局部配准
Integration 积分模块	体积积分
I/O 输入输出模块	读写3维数据
Visualization 可视化模块	使用OpenGL呈现3D数据的可自定义GUI
Utility 辅助功能模块	辅助功能，如控制台输出、文件系统和特征包装器

本文旨在快速入门Python点云数据处理，原理性知识会在后续专题中给出，这里不做过多描述。

2 安装

2.4 Open3D 0.13.0 安装

在成功安装 Anaconda3 之后，在所有应用中找到 Anaconda3，选择 Anaconda Powershell Prompt

打开之后如下所示

手动输入

pip install open3d

开始安装 open3d-0.13.0，需要漫长的等待…

安装成功！

2.5 新建一个 Python 项目

文件 > 新建项目，输入（若路径不存在，则自动创建）或选择位置路径，选择已经由Conda配置好的解释器

新建Python文件

输入新建的Python文件名，按回车完成新建

测试代码

import open3d as o3d
import numpy as np

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("test.pcd")
print(pcd)

print("->正在保存点云")
o3d.io.write_point_cloud("write.pcd", pcd, True) # 默认false，保存为Binarty；True 保存为ASICC形式
print(pcd)

3 点云读写

代码：

import open3d as o3d
import numpy as np

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("test.pcd")
print(pcd)

print("->正在保存点云")
o3d.io.write_point_cloud("write.pcd", pcd, True) # 默认false，保存为Binarty；True 保存为ASICC形式
print(pcd)

输出结果：

->正在加载点云... 
PointCloud with 356478 points.
->正在保存点云
PointCloud with 356478 points.

默认情况下，Open3D尝试通过文件扩展名推断文件类型。支持以下点云文件类型：

格式	描述
`xyz`	每一行包含 `[x, y, z]`, 其中的`x`, `y`, `z` 是三维坐标
`xyzn`	每一行包含 `[x, y, z, nx, ny, nz]`, 其中的 `nx`, `ny`, `nz` 是法线向量
`xyzrgb`	每一行包含 `[x, y, z, r, g, b]`, 其中的 `r`, `g`, `b` 是范围在 `[0, 1]`的float类型
`pts`	第一行是表示点数的整数。每个后续行包含`[x, y, z, i, r, g, b]`, 其中的 `r`, `g`, `b` 是`uint8`类型
`ply`	查看 , ply文件能够同时包含点云和网格数据

pcd = o3d.io.read_point_cloud("../../test_data/my_points.txt", format='xyz')

4 点云可视化

2.1 可视化单个点云

代码：

import open3d as o3d
import numpy as np

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("test.pcd")
print(pcd)

print("->正在可视化点云")
o3d.visualization.draw_geometries([pcd])

输出结果：

->正在加载点云... 
PointCloud with 356478 points.
->正在可视化点云

可视化结果展示：

可以通过键盘上的 +、-，实时调整点的大小。

2.2 同一窗口可视化多个点云

函数描述：

o3d.visualization.draw_geometries([pcd1, pcd2, ... ,pcdn])

代码：

import open3d as o3d
import numpy as np

print("->正在加载点云1... ")
pcd1 = o3d.io.read_point_cloud("bunny.pcd")
print(pcd1)

print("->正在加载点云2...")
pcd2 = o3d.io.read_point_cloud("bunny0.pcd")
print(pcd2)

print("->正在同时可视化两个点云...")
o3d.visualization.draw_geometries([pcd1, pcd2])

输出结果：

->正在加载点云1... 
PointCloud with 35947 points.
->正在加载点云2...
PointCloud with 35947 points.
->正在同时可视化两个点云...

可视化结果：

2.3 可视化的属性设置

函数原型1：

draw_geometries(geometry_list, window_name='Open3D', width=1920, height=1080, left=50, top=50, point_show_normal=False, mesh_show_wireframe=False, mesh_show_back_face=False)

参数说明

参数名	描述
`geometry_list`	要可视化的几何体列表
`window_name`	`(str, optional, default='Open3D')`，可视化窗口的显示标题
`width`	`(int, optional, default=1920)`，可视化窗口的宽度
`height`	`(int, optional, default=1080)`，可视化窗口的高度
`left`	`(int, optional, default=50)`，可视化窗口的左边距
`top`	`(int, optional, default=50)`，可视化窗口的顶部边距
`point_show_normal`	`(bool, optional, default=False)`，如果设置为true，则可视化点法线，需要事先计算点云法线
`mesh_show_wireframe`	`(bool, optional, default=False)`，如果设置为true，则可视化网格线框
`mesh_show_back_face`	`(bool, optional, default=False)`，可视化网格三角形的背面

函数原型2：

draw_geometries(geometry_list, window_name='Open3D', width=1920, height=1080, left=50, top=50, point_show_normal=False, mesh_show_wireframe=False, mesh_show_back_face=False, lookat, up, front, zoom)

参数说明

参数名	描述
`geometry_list`	要可视化的几何体列表
`window_name`	`(str, optional, default='Open3D')`，可视化窗口的显示标题
`width`	`(int, optional, default=1920)`，可视化窗口的宽度
`height`	`(int, optional, default=1080)`，可视化窗口的高度
`left`	`(int, optional, default=50)`，可视化窗口的左边距
`top`	`(int, optional, default=50)`，可视化窗口的顶部边距
`point_show_normal`	`(bool, optional, default=False)`，如果设置为true，则可视化点法线
`mesh_show_wireframe`	`(bool, optional, default=False)`，如果设置为true，则可视化网格线框
`mesh_show_back_face`	`(bool, optional, default=False)`，可视化网格三角形的背面
`lookat`	`(numpy.ndarray[float64[3, 1]])`，相机的注视向量
`up`	`(numpy.ndarray[float64[3, 1]])`，相机的上方向向量
`front`	`(numpy.ndarray[float64[3, 1]])`，相机的前矢量
`zoom`	`(float)`，相机的缩放倍数

代码：

import open3d as o3d

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print(pcd)

# 法线估计
radius = 0.01   # 搜索半径
max_nn = 30     # 邻域内用于估算法线的最大点数
pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius, max_nn))     # 执行法线估计

# 可视化
o3d.visualization.draw_geometries([pcd], 
                                  window_name = "可视化参数设置",
                                  width = 600,
                                  height = 450,
                                  left = 30,
                                  top = 30,
                                  point_show_normal = True)

可视化结果：

5 k-d tree 与 Octree

5.1 k-d tree

o3d.geometry.KDTreeFlann(pcd)：建立 KDTree- search_knn_vector_3d(search_Pt, k)：K 近邻搜索- search_radius_vector_3d(search_Pt, radius)：半径 R 近邻搜索- search_hybrid_vector_3d(search_pt, radius, max_nn)：混合邻域搜索，返回半径 radius 内不超过 max_nn 个近邻点代码：

import open3d as o3d
import numpy as np


print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print(pcd)

# 将点云设置为灰色
pcd.paint_uniform_color([0.5, 0.5, 0.5])

# 建立KDTree
pcd_tree = o3d.geometry.KDTreeFlann(pcd)

# 将第1500个点设置为紫色
pcd.colors[1500] = [0.5, 0, 0.5]

# 使用K近邻，将第1500个点最近的5000个点设置为蓝色
print("使用K近邻，将第1500个点最近的5000个点设置为蓝色")
k = 5000    # 设置K的大小
[num_k, idx_k, _] = pcd_tree.search_knn_vector_3d(pcd.points[1500], k)    # 返回邻域点的个数和索引
np.asarray(pcd.colors)[idx_k[1:], :] = [0, 0, 1]  # 跳过最近邻点（查询点本身）进行赋色
print("k邻域内的点数为：", num_k)

# 使用半径R近邻，将第1500个点半径（0.02）范围内的点设置为红色
print("使用半径R近邻，将第1500个点半径（0.02）范围内的点设置为红色")
radius = 0.02   # 设置半径大小
[num_radius, idx_radius, _] = pcd_tree.search_radius_vector_3d(pcd.points[1500], radius)   # 返回邻域点的个数和索引
np.asarray(pcd.colors)[idx_radius[1:], :] = [1, 0, 0]  # 跳过最近邻点（查询点本身）进行赋色
print("半径r邻域内的点数为：", num_radius)

# 使用混合邻域，将半径R邻域内不超过max_num个点设置为绿色
print("使用混合邻域，将第1500个点半径R邻域内不超过max_num个点设置为绿色")
max_nn = 200   # 半径R邻域内最大点数
[num_hybrid, idx_hybrid, _] = pcd_tree.search_hybrid_vector_3d(pcd.points[1500], radius, max_nn)
np.asarray(pcd.colors)[idx_hybrid[1:], :] = [0, 1, 0]  # 跳过最近邻点（查询点本身）进行赋色
print("混合邻域内的点数为：", num_hybrid)

print("->正在可视化点云...")
o3d.visualization.draw_geometries([pcd])

输出结果：

->正在加载点云... 
PointCloud with 35947 points.
使用K近邻，将第1500个点最近的5000个点设置为蓝色
k邻域内的点数为：5000
使用半径R近邻，将第1500个点半径（0.02）范围内的点设置为红色
半径r邻域内的点数为：751
使用混合邻域，将第1500个点半径R邻域内不超过max_num个点设置为绿色
混合邻域内的点数为：200
->正在可视化点云...

结果展示：

5.2 Octree

八叉树（Octree）是一种树数据结构，其中每个内部节点有八个子节点。八叉树通常用于三维点云的空间划分。八叉树的非空叶节点包含属于同一空间细分的一个或多个点。八叉树是三维空间的有用描述，可用于快速查找附近的点。Open3D具有几何体类型的八叉树，可用于创建、搜索和遍历具有用户指定的最大树深度的八叉树max_depth。

5.2.1 从点云中构建Octree

可以使用convert_from_point_cloud方法从点云构造八叉树。通过遵循从根节点到“最大深度”（depth max_depth）处的相应叶节点的路径，将每个点插入到树中。随着树深度的增加，内部（最终是叶子）节点表示三维空间的较小分区。

如果点云具有颜色，则对应的叶节点将采用上次插入点的颜色。size_expand参数会增加根八叉树节点的大小，使其略大于原始点云边界以容纳所有点。

代码：

import open3d as o3d
import numpy as np

# --------------------------- 加载点云 ---------------------------
print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("tree.pcd")
print("原始点云：", pcd)
# ==============================================================

# ------------------------- 构建Octree --------------------------
print('octree 分割')
octree = o3d.geometry.Octree(max_depth=4)
octree.convert_from_point_cloud(pcd, size_expand=0.01)
print("->正在可视化Octree...")
o3d.visualization.draw_geometries([octree])
# ==============================================================

输出结果：

->正在加载点云... 
原始点云：PointCloud with 5746 points.
octree 分割
->正在可视化Octree...

结果展示：

5.2.2 从体素栅格中构建Octree

还可以使用create_from_voxel_grid的方法，从Open3D体素网格VoxelGrid几何体构建八叉树。每个体素被视为三维空间中的一个点，坐标对应于体素的原点。每个叶节点采用其相应体素的颜色。

代码：

import open3d as o3d
import numpy as np

# --------------------------- 加载点云 ---------------------------
print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("tree.pcd")
print("原始点云：", pcd)
# ==============================================================

# ------------------------- 构建Octree --------------------------
print('体素化')
voxel_grid = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd, voxel_size=0.2)
print("体素：", voxel_grid)
print("正在可视化体素...")
o3d.visualization.draw_geometries([voxel_grid])

print('Octree 分割')
octree = o3d.geometry.Octree(max_depth=4)
octree.create_from_voxel_grid(voxel_grid)
print("Octree：", octree)
print("正在可视化Octree...")
o3d.visualization.draw_geometries([octree])
# ==============================================================

输出结果：

->正在加载点云... 
原始点云：PointCloud with 5746 points.
体素化
体素：VoxelGrid with 861 voxels.
正在可视化体素...
Octree 分割
Octree：Octree with origin: [64.647, -54.2659, -20.2166], size: 7, max_depth: 4
正在可视化Octree...

可视化结果：||

6 点云滤波

6.1 体素下采样

体素下采样使用规则体素栅格从输入点云创建均匀下采样点云。该算法分两步操作：

将点云进行进行体素划分- 对所有非空体素，取体素内点云的质心作为该体素的点的位置。代码：

import open3d as o3d
import numpy as np

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("test.pcd")
print(pcd)

print("->正在可视化原始点云")
o3d.visualization.draw_geometries([pcd])

print("->正在体素下采样...")
voxel_size = 0.5
downpcd = pcd.voxel_down_sample(voxel_size)
print(downpcd)

print("->正在可视化下采样点云")
o3d.visualization.draw_geometries([downpcd])

输出结果：

->正在加载点云... 
PointCloud with 356478 points.
->正在可视化原始点云
->正在体素下采样...
PointCloud with 11141 points.
->正在可视化下采样点云

可视化展示：||

6.2 统计滤波

statistical_outlier_removal 会移除距离相邻点更远的点。它需要两个输入参数：

num_neighbors，指定为了计算给定点的平均距离，需要考虑多少个邻居。即K邻域的点数- std_ratio，允许根据点云平均距离的标准偏差设置阈值水平。该数值越低，滤除的点数就越多代码：

import open3d as o3d
import numpy as np

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("desk.pcd")
print("原始点云：", pcd)

# ------------------------- 统计滤波 --------------------------
print("->正在进行统计滤波...")
num_neighbors = 20 # K邻域点的个数
std_ratio = 2.0 # 标准差乘数
# 执行统计滤波，返回滤波后的点云sor_pcd和对应的索引ind
sor_pcd, ind = pcd.remove_statistical_outlier(num_neighbors, std_ratio)
sor_pcd.paint_uniform_color([0, 0, 1])
print("统计滤波后的点云：", sor_pcd)
sor_pcd.paint_uniform_color([0, 0, 1])
# 提取噪声点云
sor_noise_pcd = pcd.select_by_index(ind,invert = True)
print("噪声点云：", sor_noise_pcd)
sor_noise_pcd.paint_uniform_color([1, 0, 0])
# ===========================================================

# 可视化统计滤波后的点云和噪声点云
o3d.visualization.draw_geometries([sor_pcd, sor_noise_pcd])

输出结果：

->正在加载点云... 
原始点云：PointCloud with 41049 points.
->正在进行统计滤波...
统计滤波后的点云：PointCloud with 40061 points.
噪声点云：PointCloud with 988 points.

可视化结果：

6.3 半径滤波

radius_outlier_removal 移除给定球体中几乎没有邻居的点。需要两个参数：

num_points，邻域球内的最少点数，低于该值的点为噪声点

radius ，邻域半径的大小

代码：

import open3d as o3d
import numpy as np


print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("desk.pcd")
print("原始点云：", pcd)

# ------------------------- 半径滤波 --------------------------
print("->正在进行半径滤波...")
num_points = 20  # 邻域球内的最少点数，低于该值的点为噪声点
radius = 0.05    # 邻域半径大小
# 执行半径滤波，返回滤波后的点云sor_pcd和对应的索引ind
sor_pcd, ind = pcd.remove_radius_outlier(num_points, radius)
sor_pcd.paint_uniform_color([0, 0, 1])
print("半径滤波后的点云：", sor_pcd)
sor_pcd.paint_uniform_color([0, 0, 1])
# 提取噪声点云
sor_noise_pcd = pcd.select_by_index(ind,invert = True)
print("噪声点云：", sor_noise_pcd)
sor_noise_pcd.paint_uniform_color([1, 0, 0])
# ===========================================================

# 可视化半径滤波后的点云和噪声点云
o3d.visualization.draw_geometries([sor_pcd, sor_noise_pcd])

输出结果：

->正在加载点云... 
原始点云：PointCloud with 41049 points.
->正在进行半径滤波...
半径滤波后的点云：PointCloud with 40146 points.
噪声点云：PointCloud with 903 points.

可视化结果：

7 点云特征提取

7.1 法线估计

代码：

import open3d as o3d
import numpy as np

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print(pcd)

print("->正在估计法线并可视化...")
radius = 0.01   # 搜索半径
max_nn = 30     # 邻域内用于估算法线的最大点数
pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius, max_nn))     # 执行法线估计
o3d.visualization.draw_geometries([pcd], point_show_normal=True)

print("->正在打印前10个点的法向量...")
print(np.asarray(pcd.normals)[:10, :])

输出结果：

->正在加载点云... 
PointCloud with 35947 points.
->正在估计法线并可视化...
->正在打印前10个点的法向量...
[[-0.22344398 -0.96962557  0.09949394]
 [-0.30282456 -0.91827757  0.25507564]
 [-0.0930339  -0.77633579 -0.62341594]
 [ 0.06452443 -0.96881599 -0.23923249]
 [ 0.24771039 -0.96349484  0.10157387]
 [ 0.1890532  -0.97541781  0.11322096]
 [-0.26920394 -0.95010988  0.15754506]
 [ 0.72941317  0.51298568  0.45255067]
 [ 0.83949302  0.5402317   0.05831961]
 [-0.32325253  0.62920765  0.7068278 ]]

可视化结果：

8 点云分割

8.1 DBSCAN 聚类分割

Open3D实现了DBSCAN[1996]，这是一种基于密度的聚类算法。该算法在cluster_dbscan中实现，需要两个参数：eps 为同一簇内的最大点间距，min_points 定义有效聚类的最小点数。函数返回标签 label，其中label = -1表示噪声。

代码：

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

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("test.pcd")
print(pcd)

print("->正在DBSCAN聚类...")
eps = 0.5           # 同一聚类中最大点间距
min_points = 50     # 有效聚类的最小点数
with o3d.utility.VerbosityContextManager(o3d.utility.VerbosityLevel.Debug) as cm:
    labels = np.array(pcd.cluster_dbscan(eps, min_points, print_progress=True))
max_label = labels.max()    # 获取聚类标签的最大值 [-1,0,1,2,...,max_label]，label = -1 为噪声，因此总聚类个数为 max_label + 1
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  # labels = -1 的簇为噪声，以黑色显示
pcd.colors = o3d.utility.Vector3dVector(colors[:, :3])
o3d.visualization.draw_geometries([pcd])

输出结果：

->正在加载点云... 
PointCloud with 356478 points.
->正在DBSCAN聚类...
[Open3D DEBUG] Precompute neighbors.
[Open3D DEBUG] Done Precompute neighbors.
[Open3D DEBUG] Compute Clusters
Precompute neighbors.[========================================] 100%
[Open3D DEBUG] Done Compute Clusters: 49
point cloud has 49 clusters

可视化结果：

8.2 RANSAC 平面分割

代码：

import open3d as o3d

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("test.pcd")
print(pcd)

print("->正在RANSAC平面分割...")
distance_threshold = 0.2    # 内点到平面模型的最大距离
ransac_n = 3                # 用于拟合平面的采样点数
num_iterations = 1000       # 最大迭代次数

# 返回模型系数plane_model和内点索引inliers，并赋值
plane_model, inliers = pcd.segment_plane(distance_threshold, ransac_n, num_iterations)

# 输出平面方程
[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([0, 0, 1.0])
print(inlier_cloud)

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

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

输出结果：

->正在加载点云... 
PointCloud with 356478 points.
->正在RANSAC平面分割...
Plane equation: -0.00x + -0.00y + 1.00z + -0.27 = 0
PointCloud with 241099 points.
PointCloud with 115379 points.

可视化结果：

8.3 隐藏点剔除

假设您希望从给定的视点渲染点云，但背景中的点会泄漏到前景中，因为它们不会被其他点遮挡。为此，我们可以应用隐藏点移除算法。在Open3D中，实现了[Katz2007] 的方法，该方法从给定视图近似点云的可见性，无需曲面重建或法线估计。

代码：

import open3d as o3d
import numpy as np

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print(pcd)

print("->正在剔除隐藏点...")
diameter = np.linalg.norm(np.asarray(pcd.get_max_bound()) - np.asarray(pcd.get_min_bound()))
print("定义隐藏点去除的参数")
camera = [0, 0, diameter]       # 视点位置
radius = diameter * 100         # 噪声点云半径,The radius of the sperical projection
_, pt_map = pcd.hidden_point_removal(camera, radius)    # 获取视点位置能看到的所有点的索引 pt_map

# 可视点点云
pcd_visible = pcd.select_by_index(pt_map)
pcd_visible.paint_uniform_color([0, 0, 1]) # 可视点为蓝色
print("->可视点个数为：", pcd_visible)
# 隐藏点点云
pcd_hidden = pcd.select_by_index(pt_map, invert = True)
pcd_hidden.paint_uniform_color([1, 0, 0]) # 隐藏点为红色
print("->隐藏点个数为：", pcd_hidden)
print("->正在可视化可视点和隐藏点点云...")
o3d.visualization.draw_geometries([pcd_visible, pcd_hidden])

输出结果：

->正在加载点云... 
PointCloud with 35947 points.
->正在剔除隐藏点...
定义隐藏点去除的参数
->可视点个数为：PointCloud with 11490 points.
->正在可视化可视点和隐藏点点云...

可视化结果： 函数原型：

    def hidden_point_removal(self, camera_location, radius): # real signature unknown; restored from __doc__
        """
        hidden_point_removal(self, camera_location, radius)
        Removes hidden points from a point cloud and returns a mesh of the remaining points. Based on Katz et al. 'Direct Visibility of Point Sets', 2007. Additional information about the choice of radius for noisy point clouds can be found in Mehra et. al. 'Visibility of Noisy Point Cloud Data', 2010.
        
        Args:
            camera_location (numpy.ndarray[float64[3, 1]]): All points not visible from that location will be reomved
            radius (float): The radius of the sperical projection
        
        Returns:
            Tuple[open3d.geometry.TriangleMesh, List[int]]
        """
        pass

9 点云曲面重建

在许多情况下，我们希望生成密集的三维几何体，即三角形网格。然而，从多视图立体方法或深度传感器中，我们只能获得非结构化点云。要从非结构化输入中获得三角形网格，我们需要执行曲面重建。

9.1 Alpha shapes

Alpha shapes 是凸壳的推广。正如这里所描述的，我们可以直观地认为Alpha shapes如下：想象一个巨大的冰淇淋，其中包含点S作为硬巧克力块。使用其中一个球形冰淇淋勺，我们可以在不撞到巧克力块的情况下雕刻出冰淇淋块的所有部分，从而甚至在内部雕刻出孔（例如，仅从外部移动勺子无法触及的部分）。我们最终将得到一个（不一定是凸的）对象，该对象以帽、弧和点为边界。如果我们现在把所有的面拉直成三角形和线段，我们就可以直观地描述S的Alpha shapes。

Open3D实现了 create_from_point_cloud_alpha_shape 方法

代码：

import open3d as o3d
import numpy as np

# --------------------------- 加载点云 ---------------------------
print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print("原始点云：", pcd)
# ==============================================================

# ------------------------- Alpha shapes -----------------------
alpha = 0.03
print(f"alpha={alpha:.3f}")
mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_alpha_shape(pcd, alpha)
mesh.compute_vertex_normals()
o3d.visualization.draw_geometries([mesh], mesh_show_back_face=True)
# ==============================================================

可视化结果：|||

9.2 Ball pivoting

Ball pivoting（BPA）[Bernardini1999]是一种与Alpha shapes相关的曲面重建方法。直观地说，想象一个具有给定半径的3D球，我们将其落在点云上。如果它击中任何3个点（并且没有穿过这3个点），它将创建一个三角形。然后，该算法从现有三角形的边开始旋转，每当它击中球未落下的3个点时，我们创建另一个三角形。

open3D 对应的函数为 create_from_point_cloud_ball_pivoting

代码：

import open3d as o3d
import numpy as np

# ---------------------- 定义点云体素化函数 ----------------------
def get_mesh(_relative_path):
    mesh = o3d.io.read_triangle_mesh(_relative_path)
    mesh.compute_vertex_normals()
    return mesh
# =============================================================

# ------------------------- Ball pivoting --------------------------
print("->Ball pivoting...")
_relative_path = "bunny.ply"    # 设置相对路径
N = 2000                        # 将点划分为N个体素
pcd = get_mesh(_relative_path).sample_points_poisson_disk(N)
o3d.visualization.draw_geometries([pcd])

radii = [0.005, 0.01, 0.02, 0.04]
rec_mesh = o3d.geometry.TriangleMesh.create_from_point_cloud_ball_pivoting(pcd, o3d.utility.DoubleVector(radii))
o3d.visualization.draw_geometries([pcd, rec_mesh])
# ==============================================================

可视化结果：||

9.3 Poisson surface reconstruction

泊松曲面重建方法[Kazhdan2006]解决了一个正则化优化问题，以获得光滑曲面。因此，泊松曲面重建比上述方法更可取，因为它们会产生非平滑结果，因为点云的点也是生成的三角形网格的顶点，无需任何修改。

Open3D对应的方法为 create_from_point_cloud_poisson，这基本上是Kazhdan代码的包装。该函数的一个重要参数是depth，它定义了用于曲面重建的八叉树的深度，因此表示生成的三角形网格的分辨率。depth值越高，网格的细节就越多。

9.3.1 直接读取点云的方法

代码：

import open3d as o3d
import numpy as np

# --------------------------- 加载点云 ---------------------------
print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print("原始点云：", pcd)
# ==============================================================

# ------------------------- Ball pivoting --------------------------
print('run Poisson surface reconstruction')
radius = 0.001   # 搜索半径
max_nn = 30         # 邻域内用于估算法线的最大点数
pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius, max_nn))     # 执行法线估计

with o3d.utility.VerbosityContextManager(o3d.utility.VerbosityLevel.Debug) as cm:
    mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(pcd, depth=12)
print(mesh)
o3d.visualization.draw_geometries([mesh])
# ==============================================================

可视化结果：

9.3.2 mesh方法

代码：

import open3d as o3d
import numpy as np

# -------------------------- 定义点云体素化函数 ---------------------------
def get_mesh(_relative_path):
    mesh = o3d.io.read_triangle_mesh(_relative_path)
    mesh.compute_vertex_normals()
    return mesh
# =====================================================================

# -------------------- Poisson surface reconstruction ------------------
# 加载点云
print("->Poisson surface reconstruction...")
_relative_path = "bunny.ply"    # 设置相对路径
N = 5000                        # 将点划分为N个体素
pcd = get_mesh(_relative_path).sample_points_poisson_disk(N)
pcd.normals = o3d.utility.Vector3dVector(np.zeros((1, 3)))  # 使现有法线无效

# 法线估计
pcd.estimate_normals()
o3d.visualization.draw_geometries([pcd], point_show_normal=True)
pcd.orient_normals_consistent_tangent_plane(100)
o3d.visualization.draw_geometries([pcd], point_show_normal=True)

# 泊松重建
print('run Poisson surface reconstruction')
with o3d.utility.VerbosityContextManager(
        o3d.utility.VerbosityLevel.Debug) as cm:
    mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(
        pcd, depth=9)
print(mesh)
o3d.visualization.draw_geometries([mesh])
# =====================================================================

可视化结果：|||

10 点云空间变换

Open3D的几何体类型有多种转换方法。下面将介绍如何使用taranslate平移、rotate旋转、scale缩放和transform变换（旋转+平移）。

10.1 Translate 平移

我们要研究的第一种转换方法是translate。该函数接受两个输入参数：

第一个输入参数为一个三维行向量(tx,ty,tz)- 第二个参数为布尔值，默认relative = True，实现点云平移。若设置为relative = False，则将点云质心平移到第一个参数指定的位置

pcd.translate((tx,ty,tz),relative=True)

10.1.1 指定平移行向量，实现点云平移

只输入第一个参数，或者将第二个参数设置为relative = True或True，relative可以省略

代码：

import copy  # 点云深拷贝
import open3d as o3d

# -------------------------- 加载点云 ------------------------
print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print(pcd)
print(f'pcd质心：{pcd.get_center()}')
# ===========================================================

# -------------------------- 点云平移 ------------------------
print("\n->沿X轴平移0.2m")
pcd_tx = copy.deepcopy(pcd).translate((0.2, 0, 0))
pcd_tx.paint_uniform_color([1, 0, 0])
print(pcd_tx)
print(f'pcd_tx质心：{pcd_tx.get_center()}')

print("\n->沿Y轴平移0.2m")
pcd_ty = copy.deepcopy(pcd_tx).translate((0, 0.2, 0))
pcd_ty.paint_uniform_color([0, 1, 0])
print(pcd_ty)
print(f'pcd_ty质心：{pcd_ty.get_center()}')

print("\n->沿X轴平移-0.2m，再沿Y轴平移0.2m")
pcd_txy = copy.deepcopy(pcd).translate((-0.2, 0.2, 0))
pcd_txy.paint_uniform_color([0, 0, 1])
print(pcd_txy)
print('pcd_txy质心：', pcd_txy.get_center())
# ===========================================================

# -------------------------- 可视化 --------------------------
o3d.visualization.draw_geometries([pcd, pcd_tx, pcd_ty, pcd_txy])
# ===========================================================

输出结果：

->正在加载点云... 
PointCloud with 35947 points.
pcd质心：[-0.02675991  0.09521606  0.00894711]

->沿X轴平移0.2m
PointCloud with 35947 points.
pcd_tx质心：[0.17324009 0.09521606 0.00894711]

->沿Y轴平移0.2m
PointCloud with 35947 points.
pcd_ty质心：[0.17324009 0.29521606 0.00894711]

->沿X轴平移-0.2m，再沿Y轴平移0.2m
PointCloud with 35947 points.
pcd_txy质心： [-0.22675991  0.29521606  0.00894711]

结果展示：

10.1.2 将点云质心平移到指定位置

将第二个参数relative设置为False即可。

代码：

import copy  # 点云深拷贝
import open3d as o3d

# -------------------------- 加载点云 ------------------------
print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print(pcd)
print(f'pcd质心：{pcd.get_center()}')
# ===========================================================

# -------------------------- 点云平移 ------------------------
print("\n->沿X轴平移0.2m")
pcd_tx = copy.deepcopy(pcd).translate((0.2, 0.2, 0.2))
pcd_tx.paint_uniform_color([1, 0, 0])
print(pcd_tx)
print(f'pcd_tx质心：{pcd_tx.get_center()}')

print("\n->将点云质心平移到指定位置")
pcd_new = copy.deepcopy(pcd_tx).translate((0.2, 0.2, 0.2),relative=False)
#pcd_new = copy.deepcopy(pcd_tx).translate((0.2, 0.2, 0.2),False)  # relative 可以省略
pcd_new.paint_uniform_color([0, 1, 0])
print(pcd_new)
print(f'pcd_new：{pcd_new.get_center()}')
# ===========================================================

# -------------------------- 可视化 --------------------------
o3d.visualization.draw_geometries([pcd, pcd_tx, pcd_new])
# ===========================================================

输出结果：

->正在加载点云... 
PointCloud with 35947 points.
pcd质心：[-0.02675991  0.09521606  0.00894711]

->沿X轴平移0.2m
PointCloud with 35947 points.
pcd_tx质心：[0.17324009 0.29521606 0.20894711]

->将点云质心平移到指定位置
PointCloud with 35947 points.
pcd_new：[0.2 0.2 0.2]

结果展示：

10.2 Rotation 旋转

Open3D使用rotate进行旋转，接受两个输入参数：

第一个参数是旋转矩阵`R`。由于3D中的旋转可以通过多种方式进行参数化，Open3D提供了方便的功能，可以将不同的参数化转换为旋转矩阵：

1. 使用`get_rotation_matrix_from_xyz`从欧拉角  中转换（其中，`xyz`也可以是`yzx`, `zxy`, `xzy`, `zyx`, 和`yxz`的形式）1. 使用`get_rotation_matrix_from_axis_angle`从轴角表示法  中转换1. 使用`get_rotation_matrix_from_quaternion`从四参数  中转换

1.  第二个参数`center`，若不设置，则围绕点云质心旋转，**旋转前后点云质心位置不变**。

pcd.rotate(R)     # 不指定旋转中心

若指定center，则整个点云围绕指定的坐标中心(x,y,z)旋转。旋转前后点云质心位置发生改变。

pcd.rotate(R,center=(x,y,z)) # 指定旋转中心

10.2.1 使用欧拉角旋转

10.2.1.1 未指定旋转中心

未指定旋转中心，默认以点云质心为旋转中心，旋转前后点云执行不变。

代码：

import copy  # 点云深拷贝
import open3d as o3d
import numpy as np

# -------------------------- 加载点云 ------------------------
print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print(pcd)
pcd.paint_uniform_color([1,0,0])
print("->pcd质心：",pcd.get_center())
# ===========================================================

# -------------------------- 点云旋转 ------------------------
print("\n->采用欧拉角进行点云旋转")
pcd_EulerAngle = copy.deepcopy(pcd)
R1 = pcd.get_rotation_matrix_from_xyz((0,np.pi/2,0))
print("旋转矩阵：\n",R1)
pcd_EulerAngle.rotate(R1)    # 不指定旋转中心
pcd_EulerAngle.paint_uniform_color([0,0,1])
print("\n->pcd_EulerAngle质心：",pcd_EulerAngle.get_center())
# ===========================================================

# -------------------------- 可视化 --------------------------
o3d.visualization.draw_geometries([pcd, pcd_EulerAngle])
# ===========================================================

输出结果：

->正在加载点云... 
PointCloud with 35947 points.
->pcd质心： [-0.02675991  0.09521606  0.00894711]

->采用欧拉角进行点云旋转
旋转矩阵：
 [[ 6.123234e-17  0.000000e+00  1.000000e+00]
 [ 0.000000e+00  1.000000e+00  0.000000e+00]
 [-1.000000e+00  0.000000e+00  6.123234e-17]]

->pcd_EulerAngle质心： [-0.02675991  0.09521606  0.00894711]

结果展示：

10.2.1.2 指定旋转中心

指定旋转中心，旋转前后点云质心改变。

代码：

import copy  # 点云深拷贝
import open3d as o3d
import numpy as np

# -------------------------- 加载点云 ------------------------
print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print(pcd)
pcd.paint_uniform_color([1,0,0])
print("->pcd质心：",pcd.get_center())
# ===========================================================

# -------------------------- 点云旋转 ------------------------
print("\n->采用欧拉角进行点云旋转")
pcd_EulerAngle = copy.deepcopy(pcd)
R1 = pcd.get_rotation_matrix_from_xyz((0,np.pi/2,0))
print("旋转矩阵：\n",R1)
pcd_EulerAngle.rotate(R1,center = (0.1,0.1,0.1))    # 指定旋转中心
pcd_EulerAngle.paint_uniform_color([0,0,1])
print("\n->pcd_EulerAngle质心：",pcd_EulerAngle.get_center())
# ===========================================================

# -------------------------- 可视化 --------------------------
o3d.visualization.draw_geometries([pcd, pcd_EulerAngle])
# ===========================================================

输出结果：

->正在加载点云... 
PointCloud with 35947 points.
->pcd质心： [-0.02675991  0.09521606  0.00894711]

->采用欧拉角进行点云旋转
旋转矩阵：
 [[ 6.123234e-17  0.000000e+00  1.000000e+00]
 [ 0.000000e+00  1.000000e+00  0.000000e+00]
 [-1.000000e+00  0.000000e+00  6.123234e-17]]

->pcd_EulerAngle质心： [0.00894711 0.09521606 0.22675991]

结果展示：

10.2.2 使用轴向角旋转

10.2.3 使用四元数旋转

10.3 Scale 缩放

使用rotate函数实现点云缩放，接受两个输入参数：

pcd_scale.rotate(2,center=pcd.get_center()) # 缩放前后质心一致

pcd_scale.rotate(2,center=(x,y,z))   # 缩放后质心为(x,y,z)

代码：

import copy  # 点云深拷贝
import open3d as o3d
import numpy as np

# -------------------------- 加载点云 ------------------------
print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("circle.pcd")
print(pcd)
pcd.paint_uniform_color([1,0,0])
print("->pcd质心：",pcd.get_center())
# ===========================================================

# -------------------------- 点云缩放 ------------------------
print("\n->点云缩放")
pcd_scale1 = copy.deepcopy(pcd)
pcd_scale1.scale(1.5,center=pcd.get_center())
pcd_scale1.paint_uniform_color([0,0,1])
print("->pcd_scale1质心：",pcd_scale1.get_center()) # 缩放前后质心不变

pcd_scale2 = copy.deepcopy(pcd)
pcd_scale2.scale(0.5,center=(1,1,1)) # 自定义缩放后的质心
pcd_scale2.paint_uniform_color([0,1,0])
print("->pcd_scale2质心：",pcd_scale2.get_center())
# ===========================================================

# -------------------------- 可视化 --------------------------
o3d.visualization.draw_geometries([pcd, pcd_scale1,pcd_scale2])
# ===========================================================

输出结果：

->正在加载点云... 
PointCloud with 63 points.
->pcd质心： [1.01066824 1.99955513 0.        ]

->点云缩放
->pcd_scale1质心： [1.01066824 1.99955513 0.        ]
->pcd_scale2质心： [1.00533412 1.49977756 0.5       ]

结果展示：

10.4 General transformation 一般变换（平移+旋转）

使用transform实现点云的一般变换，接受1个输入参数，为4×4变换矩阵，前3行3列为旋转矩阵，第4列前3行为平移向量。

代码：

import copy  # 点云深拷贝
import open3d as o3d
import numpy as np

# -------------------------- 加载点云 ------------------------
print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("circle.pcd")
print(pcd)
pcd.paint_uniform_color([1,0,0])
print("->pcd质心：",pcd.get_center())
# ===========================================================

# -------------------------- transform ------------------------
print("\n->点云的一般变换")
pcd_T = copy.deepcopy(pcd)
T = np.eye(4)
T[ :3, :3] = pcd.get_rotation_matrix_from_xyz((np.pi/6,np.pi/4,0)) # 旋转矩阵
T[0,3] = 5.0 # 平移向量的dx
T[1,3] = 3.0 # 平移向量的dy
print("\n->变换矩阵：\n",T)
pcd_T.transform(T)
pcd_T.paint_uniform_color([0,0,1])
print("\n->pcd_scale1质心：",pcd_T.get_center())
# ===========================================================

# -------------------------- 可视化 --------------------------
o3d.visualization.draw_geometries([pcd, pcd_T])
# ===========================================================

输出结果：

->正在加载点云... 
PointCloud with 63 points.
->pcd质心： [1.01066824 1.99955513 0.        ]

->点云的一般变换

->变换矩阵：
 [[ 0.70710678  0.          0.70710678  5.        ]
 [ 0.35355339  0.8660254  -0.35355339  3.        ]
 [-0.61237244  0.5         0.61237244  0.        ]
 [ 0.          0.          0.          1.        ]]

->pcd_scale1质心： [5.71465037 5.08899072 0.38087219]

结果展示：

11 点云配准

12 其他常用算法

12.1 计算点云间的距离

函数介绍：

Open3D提供了compute_point_cloud_distance方法来计算从源点云到目标点云的距离。它为源点云中的每个点计算到目标点云中最近点的距离。

代码：

import open3d as o3d
import numpy as np

print("->正在加载点云1... ")
pcd1 = o3d.io.read_point_cloud("bunny.pcd")
print(pcd1)
print("->正在加载点云2...")
pcd2 = o3d.io.read_point_cloud("test.pcd")
print(pcd2)

print("->正在点云1每一点到点云2的最近距离...")
dists = pcd1.compute_point_cloud_distance(pcd2)
dists = np.asarray(dists)
print("->正在打印前10个点...")
print(dists[:10])

print("->正在提取距离大于3.56的点")
ind = np.where(dists > 3.56)[0]
pcd3 = pcd1.select_by_index(ind)
print(pcd3)
o3d.visualization.draw_geometries([pcd3])

输出结果：

->正在加载点云1... 
PointCloud with 35947 points.
->正在加载点云2...
PointCloud with 356478 points.
->正在点云1每一点到点云2的最近距离...
->正在打印前10个点...
[3.56774778 3.57472048 3.58922401 3.53018802 3.55273519 3.55542831
 3.56749706 3.49795738 3.49527627 3.55416983]
->正在提取距离大于3.56的点
PointCloud with 19935 points.

可视化结果：

12.2 计算点云最小包围盒

函数介绍：

代码：

import open3d as o3d
import numpy as np

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print(pcd)

print("->正在计算点云轴向最小包围盒...")
aabb = pcd.get_axis_aligned_bounding_box()
aabb.color = (1, 0, 0)
print("->正在计算点云最小包围盒...")
obb = pcd.get_oriented_bounding_box()
obb.color = (0, 1, 0)
o3d.visualization.draw_geometries([pcd, aabb, obb])

输出结果：

->正在加载点云... 
PointCloud with 35947 points.
->正在计算点云轴向最小包围盒...
->正在计算点云最小包围盒...

可视化结果：

12.3 计算点云凸包

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

代码：

import open3d as o3d
import numpy as np

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print(pcd)

print("->正在计算点云凸包...")
hull, _ = pcd.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([pcd, hull_ls])

输出结果：

->正在加载点云... 
PointCloud with 35947 points.
->正在计算点云凸包...

可视化结果：

12.4 点云体素化

12.4.1 一种简单的方法

代码：

import open3d as o3d
import numpy as np

# --------------------------- 加载点云 ---------------------------
print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.ply")
print("原始点云：", pcd)
# ==============================================================

# --------------------------- 体素化点云 -------------------------
print('执行体素化点云')
voxel_grid = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd, voxel_size=0.005)
print("正在可视化体素...")
o3d.visualization.draw_geometries([voxel_grid])
# ==============================================================

输出结果：

->正在加载点云... 
原始点云：PointCloud with 35947 points.
执行体素化点云
正在可视化体素...

可视化结果：

12.1.2 复杂方法

代码：

import open3d as o3d
import numpy as np

# ---------------------- 定义点云体素化函数 ----------------------
def get_mesh(_relative_path):
    mesh = o3d.io.read_triangle_mesh(_relative_path)
    mesh.compute_vertex_normals()
    return mesh
# =============================================================

# ------------------------- 点云体素化 --------------------------
print("->正在进行点云体素化...")
_relative_path = "bunny.ply"    # 设置相对路径
N = 2000        # 将点划分为N个体素
pcd = get_mesh(_relative_path).sample_points_poisson_disk(N)

# fit to unit cube
pcd.scale(1 / np.max(pcd.get_max_bound() - pcd.get_min_bound()),
          center=pcd.get_center())
pcd.colors = o3d.utility.Vector3dVector(np.random.uniform(0, 1, size=(N, 3)))
print("体素下采样点云：", pcd)
print("正在可视化体素下采样点云...")
o3d.visualization.draw_geometries([pcd])

print('执行体素化点云')
voxel_grid = o3d.geometry.VoxelGrid.create_from_point_cloud(pcd, voxel_size=0.05)
print("正在可视化体素...")
o3d.visualization.draw_geometries([voxel_grid])
# ===========================================================

输出结果：

->正在进行点云体素化...
体素下采样点云：PointCloud with 2000 points.
正在可视化体素下采样点云...
执行体素化点云
正在可视化体素...

可视化结果：||

get_mesh 函数参考代码*：

open3d_tutorial.py

# Helpers and monkey patches for ipynb tutorials

import open3d as o3d
import numpy as np
import PIL.Image
import IPython.display
import os
import urllib
import tarfile
import gzip
import zipfile
import shutil

interactive = True


def jupyter_draw_geometries(
        geoms,
        window_name="Open3D",
        width=1920,
        height=1080,
        left=50,
        top=50,
        point_show_normal=False,
        mesh_show_wireframe=False,
        mesh_show_back_face=False,
        lookat=None,
        up=None,
        front=None,
        zoom=None,
):
    vis = o3d.visualization.Visualizer()
    vis.create_window(
        window_name=window_name,
        width=width,
        height=height,
        left=left,
        top=top,
        visible=True,  # If false, capture_screen_float_buffer() won't work.
    )
    vis.get_render_option().point_show_normal = point_show_normal
    vis.get_render_option().mesh_show_wireframe = mesh_show_wireframe
    vis.get_render_option().mesh_show_back_face = mesh_show_back_face
    for geom in geoms:
        vis.add_geometry(geom)
    if lookat is not None:
        vis.get_view_control().set_lookat(lookat)
    if up is not None:
        vis.get_view_control().set_up(up)
    if front is not None:
        vis.get_view_control().set_front(front)
    if zoom is not None:
        vis.get_view_control().set_zoom(zoom)
    if interactive:
        vis.run()
    else:
        for geom in geoms:
            vis.update_geometry(geom)
        vis.poll_events()
        vis.update_renderer()
    im = vis.capture_screen_float_buffer()
    vis.destroy_window()
    im = (255 * np.asarray(im)).astype(np.uint8)
    IPython.display.display(PIL.Image.fromarray(im, "RGB"))


o3d.visualization.draw_geometries = jupyter_draw_geometries


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 _relative_path(path):
    script_path = os.path.realpath(__file__)
    script_dir = os.path.dirname(script_path)
    return os.path.join(script_dir, path)


def download_fountain_dataset():
    fountain_path = _relative_path("../TestData/fountain_small")
    fountain_zip_path = _relative_path("../TestData/fountain.zip")
    if not os.path.exists(fountain_path):
        print("downloading fountain dataset")
        url = "https://storage.googleapis.com/isl-datasets/open3d-dev/fountain.zip"
        urllib.request.urlretrieve(url, fountain_zip_path)
        print("extract fountain dataset")
        with zipfile.ZipFile(fountain_zip_path, "r") as zip_ref:
            zip_ref.extractall(os.path.dirname(fountain_path))
        os.remove(fountain_zip_path)
    return fountain_path


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


def get_armadillo_mesh():
    armadillo_path = _relative_path("../TestData/Armadillo.ply")
    if not os.path.exists(armadillo_path):
        print("downloading armadillo mesh")
        url = "http://graphics.stanford.edu/pub/3Dscanrep/armadillo/Armadillo.ply.gz"
        urllib.request.urlretrieve(url, armadillo_path + ".gz")
        print("extract armadillo mesh")
        with gzip.open(armadillo_path + ".gz", "rb") as fin:
            with open(armadillo_path, "wb") as fout:
                shutil.copyfileobj(fin, fout)
        os.remove(armadillo_path + ".gz")
    mesh = o3d.io.read_triangle_mesh(armadillo_path)
    mesh.compute_vertex_normals()
    return mesh


def get_bunny_mesh():
    bunny_path = _relative_path("../TestData/Bunny.ply")
    if not os.path.exists(bunny_path):
        print("downloading bunny mesh")
        url = "http://graphics.stanford.edu/pub/3Dscanrep/bunny.tar.gz"
        urllib.request.urlretrieve(url, bunny_path + ".tar.gz")
        print("extract bunny mesh")
        with tarfile.open(bunny_path + ".tar.gz") as tar:
            tar.extractall(path=os.path.dirname(bunny_path))
        shutil.move(
            os.path.join(
                os.path.dirname(bunny_path),
                "bunny",
                "reconstruction",
                "bun_zipper.ply",
            ),
            bunny_path,
        )
        os.remove(bunny_path + ".tar.gz")
        shutil.rmtree(os.path.join(os.path.dirname(bunny_path), "bunny"))
    mesh = o3d.io.read_triangle_mesh(bunny_path)
    mesh.compute_vertex_normals()
    return mesh


def get_knot_mesh():
    mesh = o3d.io.read_triangle_mesh(_relative_path("../TestData/knot.ply"))
    mesh.compute_vertex_normals()
    return mesh


def get_eagle_pcd():
    path = _relative_path("../TestData/eagle.ply")
    if not os.path.exists(path):
        print("downloading eagle pcl")
        url = "http://www.cs.jhu.edu/~misha/Code/PoissonRecon/eagle.points.ply"
        urllib.request.urlretrieve(url, path)
    pcd = o3d.io.read_point_cloud(path)
    return pcd

12.5 计算点云质心

使用get_center() 实现点云质心计算

代码：

import open3d as o3d

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("tree.pcd")
print(pcd)
print(f'pcd质心：{pcd.get_center()}')

输出结果：

->正在加载点云... 
PointCloud with 5746 points.
pcd质心：[ 66.36420378 -52.42476729 -15.28512276]

：

->点云质心为：(66.3642, -52.4248, -15.2851)

12.6 根据索引提取点云

函数原型：

select_by_index(self, indices, invert=False)

代码：

import open3d as o3d
import numpy as np

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("desk.pcd")
print(pcd)

"""------------------- 根据索引提取点云 --------------------"""
print("->正在根据索引提取点云...")
idx = list(range(20000))    # 生成 从0到19999的列表

# 索引对应的点云（内点）
inlier_pcd = pcd.select_by_index(idx)
inlier_pcd.paint_uniform_color([1, 0, 0])
print("内点点云：", inlier_pcd)

# 索引外的点云（外点）
outlier_pcd = pcd.select_by_index(idx, invert=True)     # 对索引取反
outlier_pcd.paint_uniform_color([0, 1, 0])
print("外点点云：", outlier_pcd)

o3d.visualization.draw_geometries([inlier_pcd, outlier_pcd])
"""========================================================"""

输出结果：

->正在加载点云... 
PointCloud with 41049 points.
->正在根据索引提取点云...
内点点云：PointCloud with 20000 points.
外点点云：PointCloud with 21049 points.

可视化结果：

源码：

    def select_by_index(self, indices, invert=False): # real signature unknown; restored from __doc__
        """
        select_by_index(self, indices, invert=False)
        Function to select points from input pointcloud into output pointcloud.
        
        Args:
            indices (List[int]): Indices of points to be selected.
            invert (bool, optional, default=False): Set to ``True`` to invert the selection of indices.
        
        Returns:
            open3d.geometry.PointCloud
        """
        pass

12.7 点云赋色

赋色函数：

paint_uniform_color ：将点云进行单一颜色赋值，为RGB颜色空间，R、G、B分量的范围为[0,1]，注意不是[0,255]

代码：

import open3d as o3d
import numpy as np

print("->正在加载点云... ")
pcd = o3d.io.read_point_cloud("bunny.pcd")
print(pcd)

print("->正在点云赋色...")
pcd.paint_uniform_color([1,0.706,0])
print("->正在可视化赋色后的点云...")
o3d.visualization.draw_geometries([pcd])

print("->正在保存赋色后的点云")
o3d.io.write_point_cloud("color.pcd", pcd, True) # 默认false，保存为Binarty；True 保存为ASICC形式

输出结果：

->正在加载点云... 
PointCloud with 35947 points.
->正在点云赋色...
->正在可视化赋色后的点云...
->正在保存赋色后的点云

可视化结果：

缩放倍数- 缩放后点云质心center，不可省略。