|
三维点云学习（7）5-Feature 实现ISS
参考博客
ISS思路参考 图像跟踪与识别-NMS非极大值抑制学习笔记 ISS笔记回顾：三维点云学习（7）3-Feature Detection-Intrinsic Shape Signatures(ISS)
数据集下载
为40种物体的三维点云数据集 链接：https://pan.baidu.com/s/1LX9xeiXJ0t-Fne8BCGSjlQ 提取码：es14
效果图：
eg.文末有其他物体的效果图展示
airplane_0001
红色圆圈为 提取的特征点 
  可以看出特征点大致呈对称分布或分布于物体中轴线上： 
主要coding思路：
特征点（关键点）
step1 使用radius NN 得到n个初始关键点, 确定 threshold 阈值
step2 随机选取一个关键点，寻找这个关键点的nearest 中的其他关键点，构成一个关键点群。 思路类似(DBSCAN)
step3 NMS 非极大抑制，比较 在一个关键点群中，每个关键点的 linda3，最大的linda3对应的点即定为这个关键点群的关键点
step4 iearation step2-step3
代码部分模块
构建点云坐标的协方差矩阵、求解cov特征值，并把特征值降序输出
def compute_cov_eigval(point_cloud):
x = np.asarray(point_cloud[:,0])
y = np.asarray(point_cloud[:,1])
z = np.asarray(point_cloud[:,2])
M = np.vstack((x, y, z)) # 每行表示一个属性， 每列代表一个点
cov = np.cov(M) # 使用每个点的坐标求解cov
eigval, eigvec = np.linalg.eigh(cov) # 求解三个特征值，升序排列 linda1 < linda2 < linda3
eigval = eigval[np.argsort(-eigval)] # 改为降序排列 linda1 > linda2 > linda3
return eigval #返回特征值
matplotlib显示点云函数
# matplotlib显示点云函数
def Point_Cloud_Show(point_cloud,feature_point):
fig = plt.figure(dpi=150)
ax = fig.add_subplot(111, projection='3d')
ax.scatter(point_cloud[:, 0], point_cloud[:, 1], point_cloud[:, 2], cmap='spectral', s=2, linewidths=0, alpha=1, marker=".")
ax.scatter(feature_point[:, 0], feature_point[:, 1], feature_point[:, 2], cmap='spectral', s=2, linewidths=5, alpha=1,marker=".",color='red')
plt.title('Point Cloud')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
主体 iss 函数
#本次使用三个数据集分别是 airplane_0001.txt piano_0001.txt cahir_0001.txt；效果最好是airplane_0001.txt，其他效果基于参数调整不好，效果不尽人意
注意：个人觉得关键点->两个参数
linda3_threshold 阈值参数设置，方法1：中值方法
阈值参数的设置需要根据点云的输入不断调整，因为每个点云的数据点密集程度都不尽相同。 这里我的处理方法是，得到所有点的最小特征值(linda3),取所有特征值的中值5，至于为什么5这个也是因为多次调试(数据集：airplane_0001.txt)观看结果得出个人觉得比较好的结果。
# 根据linda3的数值 确定linda3_threshold(linda的阈值)
linda3_threshold = np.median(eigvals,axis=0)[2]*5 #阈值取大约 是所有linda3的 中值得5倍， eg 为什么取5倍是个人调试决定，也可取1倍
linda3_threshold 阈值参数设置，方法2：估计方法
在调整linda3_threshold 时，可输出特征值，大致估计阈值设置多少
#特征值输出，可大致确定 阈值设置为 e-04级别
[[2.28817630e-03 1.64286391e-03 9.33870633e-04]
[1.19472537e-03 8.01334793e-04 3.84488142e-05]
[2.00268312e-03 7.17766696e-04 1.19501607e-05]
...
[2.07454874e-03 1.33158344e-03 5.35126375e-04]
[2.19574935e-03 7.43165606e-04 1.37694390e-05]
[2.51036763e-03 1.44480340e-03 8.17724964e-04]]
radius NN radius的设置
radius的设置也十分重要，radius太大会导致把大部分点云覆盖，每个点nearest 无限接近 全数据集，所有点的特征值几乎相同 如下所示，故意设置radius = 10.0（too big！！！）
#构建 kd_tree
leaf_size = 4
radius = 10 # 各文件参数 airplane_00001:0.1; chair_0001:0.1
tree = KDTree(data,leaf_size)
#step1 使用radius NN 得到n个初始关键点, threshold 阈值 ：每个radius内的linda大于某个数值
nearest_idx = tree.query_radius(data,radius)
radius太大，输出结果所有点的radius nn 几乎重合，所有点的局部特征值相同
#out:
[[0.13978321 0.13373699 0.01789397]
[0.13978321 0.13373699 0.01789397]
[0.13978321 0.13373699 0.01789397]
...
[0.13978321 0.13373699 0.01789397]
[0.13978321 0.13373699 0.01789397]
[0.13978321 0.13373699 0.01789397]]
NMS
NMS 个人觉得编写思路和DBSCAN贼像，所以沿用之前DBSCAN的代码
def iss(data):
#parameters
eigvals = []
feature = []
T = set() #T 关键点的集合
linda3_threshold = None #阈值，初步筛选 ,各文件参数 airplane_0001:0.001; chair_0001:0.0001
#构建 kd_tree
leaf_size = 4
radius = 0.1 # 各文件参数 airplane_00001:0.1; chair_0001:0.1
tree = KDTree(data,leaf_size)
#step1 使用radius NN 得到n个初始关键点, threshold 阈值 ：每个radius内的linda大于某个数值
nearest_idx = tree.query_radius(data,radius)
for i in range(len(nearest_idx)):
eigvals.append(compute_cov_eigval(data[nearest_idx[i]]))
eigvals = np.asarray(eigvals) # 求解每个点在各自的 radius 范围内的linda
print(eigvals) #打印所有的 特征值，供调试用
# 根据linda3的数值 确定linda3_threshold(linda的阈值)
linda3_threshold = np.median(eigvals,axis=0)[2]*5 #阈值取大约 是所有linda3的 中值得5倍， eg 为什么取5倍是个人调试决定，也可取1倍
print(linda3_threshold)
for i in range(len(nearest_idx)):
if eigvals[i,2] > linda3_threshold: # compute_cov_eigval(data[nearest_idx[i]])[2] -> 每个radius 里的最小的特征值 linda3
T.add(i) #获得初始关键点的索引
print(T) #输出 初始关键点
#step2 有 重叠(IOU)的 关键点群
unvisited = T #未访问集合
while len(T):
unvisited_old = unvisited #更新访问集合
core = list(T)[np.random.randint(0,len(T))] #从 关键点集T 中随机选取一个 关键点core
unvisited = unvisited - set([core]) #把核心点标记为 visited,从 unvisited 集合中剔除
visited = []
visited.append(core)
while len(visited): #遍历所有初始关键点
new_core = visited[0]
if new_core in T:
S = unvisited & set(nearest_idx[new_core]) #S : 当前 关键点(core) 的范围内所包含的其他关键点
# print(T)
# print(S)
visited += (list(S))
unvisited = unvisited - S
visited.remove(new_core) #new core 已做检测，去掉new core
cluster = unvisited_old - unvisited #cluster, 有 重叠(IOU)的 关键点群
T = T - cluster #去掉该类对象里面包含的核心对象,差集
#step3 NMS 非极大抑制，求解 一个关键点群的linda3最大 为 关键点
cluster_linda3 = []
for i in list(cluster):
cluster_linda3.append(eigvals[i][2]) #获取 每个关键点 的 linda3
cluster_linda3 = np.asarray(cluster_linda3)
NMS_OUTPUT = np.argmax(cluster_linda3)
feature.append(list(cluster)[NMS_OUTPUT]) #添加到 feature 特征点数组中
#output
return feature
主要coding思路： 其他效果图
table_0001.txt 可以看出，四角的地方明显可以规划为特征点   piano_0001.txt 可以看出，四角的地方明显可以规划为特征点   chair_0001.txt  
全部代码：
#ISS.py
import open3d as o3d
import os
import numpy as np
from pyntcloud import PyntCloud
import matplotlib.pyplot as plt
from sklearn.neighbors import KDTree # KDTree 进行搜索
import random
from pandas import DataFrame
# matplotlib显示点云函数
def Point_Cloud_Show(point_cloud,feature_point):
fig = plt.figure(dpi=150)
ax = fig.add_subplot(111, projection='3d')
ax.scatter(point_cloud[:, 0], point_cloud[:, 1], point_cloud[:, 2], cmap='spectral', s=2, linewidths=0, alpha=1, marker=".")
ax.scatter(feature_point[:, 0], feature_point[:, 1], feature_point[:, 2], cmap='spectral', s=2, linewidths=5, alpha=1,marker=".",color='red')
plt.title('Point Cloud')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()
def compute_cov_eigval(point_cloud):
x = np.asarray(point_cloud[:,0])
y = np.asarray(point_cloud[:,1])
z = np.asarray(point_cloud[:,2])
M = np.vstack((x, y, z)) # 每行表示一个属性， 每列代表一个点
cov = np.cov(M) # 使用每个点的坐标求解cov
eigval, eigvec = np.linalg.eigh(cov) # 求解三个特征值，升序排列 linda1 < linda2 < linda3
eigval = eigval[np.argsort(-eigval)] # 改为降序排列 linda1 > linda2 > linda3
return eigval #返回特征值
def iss(data):
#parameters
eigvals = []
feature = []
T = set() #T 关键点的集合
linda3_threshold = None #阈值，初步筛选 ,各文件参数 airplane_0001:0.001; chair_0001:0.0001
#构建 kd_tree
leaf_size = 4
radius = 0.1 # 各文件参数 airplane_00001:0.1; chair_0001:0.1
tree = KDTree(data,leaf_size)
#step1 使用radius NN 得到n个初始关键点, threshold 阈值 ：每个radius内的linda大于某个数值
nearest_idx = tree.query_radius(data,radius)
for i in range(len(nearest_idx)):
eigvals.append(compute_cov_eigval(data[nearest_idx[i]]))
eigvals = np.asarray(eigvals) # 求解每个点在各自的 radius 范围内的linda
print(eigvals) #打印所有的 特征值，供调试用
# 根据linda3的数值 确定linda3_threshold(linda的阈值)
linda3_threshold = np.median(eigvals,axis=0)[2]*5 #阈值取大约 是所有linda3的 中值得5倍， eg 为什么取5倍是个人调试决定，也可取1倍
print(linda3_threshold)
for i in range(len(nearest_idx)):
if eigvals[i,2] > linda3_threshold: # compute_cov_eigval(data[nearest_idx[i]])[2] -> 每个radius 里的最小的特征值 linda3
T.add(i) #获得初始关键点的索引
print(T) #输出 初始关键点
#step2 有 重叠(IOU)的 关键点群
unvisited = T #未访问集合
while len(T):
unvisited_old = unvisited #更新访问集合
core = list(T)[np.random.randint(0,len(T))] #从 关键点集T 中随机选取一个 关键点core
unvisited = unvisited - set([core]) #把核心点标记为 visited,从 unvisited 集合中剔除
visited = []
visited.append(core)
while len(visited): #遍历所有初始关键点
new_core = visited[0]
if new_core in T:
S = unvisited & set(nearest_idx[new_core]) #S : 当前 关键点(core) 的范围内所包含的其他关键点
# print(T)
# print(S)
visited += (list(S))
unvisited = unvisited - S
visited.remove(new_core) #new core 已做检测，去掉new core
cluster = unvisited_old - unvisited #cluster, 有 重叠(IOU)的 关键点群
T = T - cluster #去掉该类对象里面包含的核心对象,差集
#step3 NMS 非极大抑制，求解 一个关键点群的linda3最大 为 关键点
cluster_linda3 = []
for i in list(cluster):
cluster_linda3.append(eigvals[i][2]) #获取 每个关键点 的 linda3
cluster_linda3 = np.asarray(cluster_linda3)
NMS_OUTPUT = np.argmax(cluster_linda3)
feature.append(list(cluster)[NMS_OUTPUT]) #添加到 feature 特征点数组中
#output
return feature
if __name__ == '__main__':
point_cloud = np.genfromtxt(r"airplane_0001.txt", delimiter=",")
point_cloud = point_cloud[:, 0:3] # 为 xyz的 N*3矩阵
feature_idx = iss(point_cloud)
feature_point = point_cloud[feature_idx]
print(feature_point)
Point_Cloud_Show(point_cloud,feature_point)
|