6. 较全的Open3D点云数据处理(python)

注意：以下内容来自博客爆肝5万字❤️Open3D 点云数据处理基础（Python版）_python 点云焊缝-CSDN博客，这篇博客写的全且详细，在这里是为了记笔记方便查看，并非抄袭。

1.点云的读写

代码如下：

import open3d as o3d

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd")  
    print(pcd)

输出结果如下：

如下代码：

pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd", format='xyz')

2.点云可视化

2.1 单个点云的可视化

代码如下：

import open3d as o3d

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd")
    print(pcd)
    o3d.visualization.draw_geometries([pcd])

输出结果如下：

可视化结果如下：

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

代码如下：

import open3d as o3d

if __name__ == '__main__':
    pcd1 = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd")
    pcd2 = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_1.pcd")
    #可视化代码如下
    o3d.visualization.draw_geometries([pcd1, pcd2])

可视化结果如下：

可视化属性设置：

函数原型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)

函数原型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)

代码如下：

import open3d as o3d

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.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=1000,
                                      height=800,
                                      left=300,
                                      top=300,
                                      point_show_normal=True)

可视化结果如下：

3. k_d tree 和 Octree

3.1 k_d tree

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd")
    print(pcd)

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

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

    # 将第1500个点设置为紫色
    pcd.colors[15000] = [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[15000], k)  # 返回邻域点的个数和索引
    np.asarray(pcd.colors)[idx_k[1:], :] = [0, 0, 1]  # 跳过最近邻点（查询点本身）进行赋色
    print("k邻域内的点数为：", num_k)

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

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

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

结果如下：

可视化结果如下：

3.2 Octree

3.2.1 从点云中构建Octree

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.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])

可视化结果如下：

3.2.2 从体素网格中构建Octree

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.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])

输出结果如下：

可视化结果如下：

4.点云滤波

4.1 体素下采样

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.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])

输出结果如下：

可视化结果如下：

4.2 半径滤波

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.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])

可视化结果如下：

5.点云特征提取

5.1 法线估计

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.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, :])

结果输出如下：

可视化结果如下：

6. 点云分割

6.1 DBSCAN算法

代码如下：

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

if __name__ == '__main__':
    # pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd")
    pcd = o3d.io.read_point_cloud("datas/1.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])

输出结果如下：

可视化结果如下：

6.2 RANSAC平面分割

代码如下：

import open3d as o3d

if __name__ == '__main__':
    # pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd")
    pcd = o3d.io.read_point_cloud("datas/1.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])

可视化结果如下：

6.3 隐藏点剔除

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd")
    #pcd = o3d.io.read_point_cloud("datas/1.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])

输出结果如下：

可视化结果如下：

7.点云曲面重建

7.1 Alpha shapes

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd")
    #pcd = o3d.io.read_point_cloud("datas/1.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)

可视化结果如下：

alpha=0.5

alpha=0.01

7.2 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])
# ==============================================================

可视化结果如下：

8.点云空间变换

8.1 translate 平移

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

9.点云配准

点云配准看我的另一篇博客4.点云数据的配准_点云叠加配准-CSDN博客。

10. 其他点云计算方法

10.1 计算点云间的距离

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd1 = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd")
    pcd2 = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_1.pcd")

    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 > 0.1)[0]
    pcd3 = pcd1.select_by_index(ind)
    print(pcd3)
    o3d.visualization.draw_geometries([pcd3])

输出结果如下：

可视化结果如下：

10.2 计算点云最小包围盒

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.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])

可视化结果如下：

10.3 计算点云凸包

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.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])

可视化结果如下：

10.4 点云体素化

10.4.1 简单方法

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd")
    print(pcd)

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

可视化结果如下：

10.4.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])
# ===========================================================

可视化结果如下：

10.5 计算点云质心

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.pcd")
    print(pcd)

    print(f'pcd质心：{pcd.get_center()}')

输出结果如下：

10.6 根据索引提取点云

select_by_index(self, indices, invert=False)

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.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])
    """========================================================"""

可视化结果如下：

10.7 点云赋色

代码如下：

import open3d as o3d
import numpy as np

if __name__ == '__main__':
    pcd = o3d.io.read_point_cloud("D:\AI\yq\clouds\datas\DemoICPPointClouds\\cloud_bin_0.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形式

可视化结果如下：