Open3d(二)——点云数据的操作

 亲测代码程序可运行使用，open3d版本0.13.0。

open3d数据资源下载:​​​​​​GitHub - Cobotic/Open3D: Open3D: A Modern Library for 3D Data Processing

代码执行功能有：点云可视化、体素降采样、顶点法线估计、裁剪点云、点云上色、包围框、凸包、DBSCAN 聚类、平面分割、隐点移除，详情请见代码

'''
Author: dongcidaci
Date: 2021-09-08 13:33:56
LastEditTime: 2021-09-10 14:04:46
LastEditors: Please set LastEditors
Description: In User Settings Edit
FilePath: \open3d_code\pointcloud.py
'''
import open3d as o3d
import numpy as np
#import open3d_tutorial as o3dtut
import matplotlib.pyplot as plt


#点云数据
pcd = o3d.io.read_point_cloud("test_data/fragment.ply")
print(pcd)#输出点云的个数
print(np.asarray(pcd.points))#输出点的三维坐标
#o3d.io.write_point_cloud("table_scene_mug_stereo_textured.pcd",pcd)#写入点云数据
o3d.visualization.draw_geometries([pcd],width=800,height=800)#显示点云

#体素降采样通过使用规则体素网格从输入点云创造一致化降采样点云，步骤为：
#1.把点云装进体素网格
#2.把每个被占据的体素中的点做平均，取一个精确的点
downpcd = pcd.voxel_down_sample(voxel_size=0.05)
o3d.visualization.draw_geometries([downpcd])

#顶点法线估计
#定点法线估计，按 n 查看点云法线。使用 - 和 + 可以缩放法线长度
downpcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1,max_nn=30))
#estimate_normals计算每个点的法线。这个函数通过对邻点采用协方差分析来计算主轴
#俩个参数radius=0.1和max_nn=30设置搜索半径和邻域最大点数
o3d.visualization.draw_geometries([downpcd],
                                  point_show_normal=True)
print(downpcd.normals[0])#检索估计的顶点法线
print(np.asarray(downpcd.normals)[:10,:])#通过np.asarry来转化成numpy数组，显示前十个

#裁剪点云
vol = o3d.visualization.read_selection_polygon_volume("test_data/Crop/cropped.json")#读取一个指定多边形区域的json文件
chair = vol.crop_point_cloud(pcd)#过滤掉点，只保留椅子部分
o3d.visualization.draw_geometries([chair])

#点云上色
chair.paint_uniform_color([1,0.706,0])#所有的点上统一的颜色，颜色取值[0,1]
o3d.visualization.draw_geometries([chair])

#包围框
#open3d实现了俩个包围框接口
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])


#凸包
#点云得凸包是包里所有点的最小凸集
vol = o3d.visualization.read_selection_polygon_volume("test_data/Crop/cropped.json")
chair = vol.crop_point_cloud(chair)
hull,_=chair.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([chair,hull_ls])

#DBSCAN 聚类,一种基于密度的聚类算法
#算法接口为cluster_dbscan
#俩个参数：eps表示聚类的领域距离，min_points表示聚类的最小点数
#该函数返回一个label,label为-1表示为噪声
cd = o3d.io.read_point_cloud("test_data/fragment.ply")

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])

#平面分割
#使用segement_plane函数。这个函数需要三个参数，
# destance_threshold定义了一个点到一个估计平面的最大距离，这些距离内的点被认为是内点（inlier），
# ransac_n定义了使用随机抽样估计一个平面的点的个数，
# num_iterations定义了随机平面采样和验证的频率（迭代次数）。
#函数返回（a,b,c,d）作为一个平面，对于平面上每个点(x,y,z)我们有ax+by+cz+d=0。这个函数还会返回内点索引的列表。
pcd = o3d.io.read_point_cloud("test_data/fragment.pcd")
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])

#隐点移除
#想从一个给定视角去渲染一个点云，但是背景的点在前景也会被看到，这是因为并没有其他点把它遮挡住。
#pcd = o3dtut.get_armadillo_mesh().sample_points_poisson_disk(5000)
diameter = np.linalg.norm(np.asarray(chair.get_max_bound()) - np.asarray(chair.get_min_bound()))
o3d.visualization.draw_geometries([chair])

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])