Open3D 点云分割与聚类

Open3D 点云分割与聚类：

RANSAC 平面分割：RANSAC是随机抽样一致算法。它是根据一组包含异常数据的样本数据集，计算出数据的数学模型参数，得到有效样本数据的算法。RANSAC算法的基本理论基础是大数定律，也就是当采样数达到一定数量后采样的数据就会符合它自身原有的概率属性。这是一种通过概率的方式来进行拟合。

结合前面知识点的理解：
1.以RANSAC平面分割为例，由于三个点可以确定一个平面，因此RANSAC会随机选择三个点来构建一个平面
2.并用点云中实际上有多少个点落到这个平面上来作为评估这个平面的正确程度。
3.当随机抽样的次数足够多时，我们有较大概率获得所需要的平面：A x + B y + C z + D = 0

一、基于RANSAC平面分割的过程：
1：确定拟合模型
2：确定拟合该模型每次需要的点个数
3：随机选点计算模型参数
4：根据距离d计算inlier个数
5：重复随机，保留inlier个数最多的那个模型

1、首先确定拟合数学模型：Ax+By+Cz+D=0，拟合模型需要三个点。

1. 确定拟合模型所需要的必要点个数n
   直线：空间两点确定一条直线 n=2
   球：空间四个点确定外接球心和半径 n=4
   正方体：空间了四个点确定正方体形状位置 n=4
   平面：空间三点确定平面 n=3
2. 随机选择n个点计算数学模型

2、计算inlier数量，Inlier判定方式：计算所有点与随机点计算的平面之间的距离d，如果d<D(距离阈值)，判断为内点。

3、随机次数达到给定次数，保留inlier个数最多的那个平面模型。

segment_plane函数解析：

plane_model, inliers = pcd.segment_plane(distance_threshold=0.01ransac_n=3,
num_iterations=1000)

参数：（距离d：判定是否是inlier的依据 随机次数m：停止随机的阈值）
distance_threshold：inlier的最大距离阈值
ransac_n：随机采样的平面点数
num_iterations：表示最小迭代次数。

返回值：plane_model：平面模型，即个平面方程系数（a,b,c,d），作为一个平面，对于平面上每个点(x,y,z)，我们有ax+by+cz+d=0。
inliers ：内点索引

代码：

import open3d as o3d

pc = o3d.io.read_point_cloud("pointcloud.pcd")
pc = pc.uniform_down_sample(5) #采样
o3d.visualization.draw_geometries([pc])
# 返回模型系数plane_model和内点索引inliers，并赋值
plane_model, inliers = pc.segment_plane(0.01, 3, 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 = pc.select_by_index(inliers)
inlier_cloud.paint_uniform_color([0, 0, 1.0])
# 外点
outlier_cloud = pc.select_by_index(inliers, invert=True)
# outlier_cloud.paint_uniform_color([1.0, 0, 0])
# 可视化
o3d.visualization.draw_geometries([inlier_cloud, outlier_cloud])
#保存结果
o3d.io.write_point_cloud("plane.pcd",inlier_cloud)

分割案例代码：

import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt
# Open3D 0.16版本的 api

pcd = o3d.io.read_point_cloud("sample.ply")
# 计算法线
pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=16), fast_normal_computation=True)
pcd.paint_uniform_color([0,0,1])
# #RANSAC
# [a, b, c, d], inliers = pcd.segment_plane(distance_threshold=0.01,ransac_n=3,num_iterations=1000)
# print(f"平面方程: {a:.2f}x + {b:.2f}y + {c:.2f}z + {d:.2f} = 0")
# inlier_cloud = pcd.select_by_index(inliers)
# outlier_cloud = pcd.select_by_index(inliers, invert=True)
# inlier_cloud.paint_uniform_color([1.0, 0, 0])
# outlier_cloud.paint_uniform_color([0.6, 0.6, 0.6])
# o3d.visualization.draw_geometries([inlier_cloud, outlier_cloud])

# # DBSCAN
labels = np.array(pcd.cluster_dbscan(eps=0.05, min_points=10))
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])

segment_models={}
segments={}
max_plane_idx=40
rest=pcd
d_threshold=0.01
for i in range(max_plane_idx):
    colors = plt.get_cmap("tab20")(i)
    segment_models[i], inliers = rest.segment_plane(distance_threshold=0.01,ransac_n=3,num_iterations=1000)
    segments[i]=rest.select_by_index(inliers)
    labels = np.array(segments[i].cluster_dbscan(eps=d_threshold*10, min_points=10))
    candidates=[len(np.where(labels==j)[0]) for j in np.unique(labels)]
    best_candidate=int(np.unique(labels)[np.where(candidates==np.max(candidates))[0]])
    print("the best candidate is: ", best_candidate)
    rest = rest.select_by_index(inliers, invert=True)+segments[i].select_by_index(list(np.where(labels!=best_candidate)[0]))
    segments[i]=segments[i].select_by_index(list(np.where(labels==best_candidate)[0]))
    segments[i].paint_uniform_color(list(colors[:3]))
    print("pass",i+1,"/",max_plane_idx,"done.")

    labels = np.array(rest.cluster_dbscan(eps=0.05, min_points=5))
    max_label = labels.max()
    print(f"point cloud has {max_label + 1} clusters")

    colors = plt.get_cmap("tab10")(labels / (max_label if max_label > 0 else 1))
    colors[labels < 0] = 0
    rest.colors = o3d.utility.Vector3dVector(colors[:, :3])

o3d.visualization.draw_geometries([segments[i] for i in range(max_plane_idx)]+[rest])

代码函数语法难点：
np.max 与 np.maximum

>> np.max([-2, -1, 0, 1, 2])
2

>> np.maximum([-2, -1, 0, 1, 2], 0)
array([0, 0, 0, 1, 2])

        # 当然 np.maximum 接受的两个参数，也可以大小一致
        # 或者更为准确地说，第二个参数只是一个单独的值时，其实是用到了维度的 broadcast 机制；

    # 当然 np.maximum 接受的两个参数，也可以大小一致
    # 或者更为准确地说，第二个参数只是一个单独的值时，其实是用到了维度的 broadcast 机制；

Broadcast机制(还不够了解，CSDN链接解析)
大神总结https://blog.csdn.net/qq_38683460/article/details/126331052

np.unique( )的用法
该函数是去除数组中的重复数字，并进行排序之后输出。

拓展：
np.array(png)，这个函数处理图片之后得到的结果是什么？
可以得到图像中存在的灰度值（或者像素值），然后分别统计图片里这些灰度值出现的个数，就可以得到灰度分布直方图了

numpy.where() 有两种用法：
网址链接:
大神链接https://www.cnblogs.com/massquantity/p/8908859.html
numpy网址https://numpy.org/doc/stable/reference/generated/numpy.where.html

1. np.where(condition, x, y)
   满足条件(condition)，输出x，不满足输出y。

2. np.where(condition)
   只有条件 (condition)，没有x和y，则输出满足条件 (即非0) 元素的坐标 (等价于numpy.nonzero)。这里的坐标以tuple的形式给出，通常原数组有多少维，输出的tuple中就包含几个数组，分别对应符合条件元素的各维坐标。

列表生成式:
链接：
大神链接https://www.liaoxuefeng.com/wiki/1016959663602400/1017317609699776

即List Comprehensions，是Python内置的非常简单却强大的可以用来创建list的生成式。
举个例子，要生成list [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
可以用list(range(1, 11))：