三维点云处理-特征点描述 8.1 PFHFPFH

一、特征点的描述1) PFHPFH方法核心特点实现旋转不变性捕捉局部表面几何变化基于点对法向量与相对位置计算参数复杂度缺陷n个特征点需计算nk²次运算敏感性问题结果严重依赖法向量精度2) FPFHFPFH改进方案简化计算每个点仅计算与直接邻居的连线k次运算直方图构建采用3b维拼接方案非三维投票加权融合 中心点SPFH简化直方图邻居点SPFH的加权平均复杂度优化从nk²降至nk4) fph与pfh的区别对比维度PFHFPFH邻居连接方式全连接部分连接影响范围固定半径rr至2r参数重复计算无存在时间复杂度O(nk²)O(nk)描述符长度b³3b实际应用中优先选择FPFH因其在相近效果下具备显著效率优势。二. 代码实践main.py#!/opt/conda/envs/feature-detection/bin/python# main.py# 1. load point cloud in modelnet40 normal format# 2. calculate ISS keypoints# 3. calculate FPFH or SHOT for detected keypoints# 3. visualize the resultsimportargparse# IO utils:fromutils.ioimportread_modelnet40_normal# detector:fromdetection.issimportdetect# descriptor:fromdescription.fpfhimportdescribeimportnumpyasnpimportpandasaspdimportopen3daso3dimportseabornassnsimportmatplotlib.pyplotaspltdefget_arguments(): Get command-line arguments # init parser:parserargparse.ArgumentParser(Detect ISS keypoints on ModelNet40 dataset.)# add required and optional groups:requiredparser.add_argument_group(Required)optionalparser.add_argument_group(Optional)# add required:required.add_argument(-i,destinput,helpInput path of ModelNet40 sample.,requiredTrue)required.add_argument(-r,destradius,helpRadius for radius nearest neighbor definition.,requiredTrue,typefloat)# parse arguments:returnparser.parse_args()if__name____main__:# parse arguments:argumentsget_arguments()# load point cloud:point_cloudread_modelnet40_normal(arguments.input)# compute surface normals:# point_cloud.estimate_normals(# search_paramo3d.geometry.KDTreeSearchParamHybrid(radius0.1, max_nn30)# )# build search tree:search_treeo3d.geometry.KDTreeFlann(point_cloud)# detect keypoints:keypointsdetect(point_cloud,search_tree,arguments.radius)# visualize:# paint background as grey:point_cloud.paint_uniform_color([0.50,0.50,0.50])# show roi:max_boundpoint_cloud.get_max_bound()min_boundpoint_cloud.get_min_bound()center(min_boundmax_bound)/2.0min_bound[1]max_bound[1]-0.1max_bound[1]max_bound[1]min_bound[2]center[2]max_bound[2]max_bound[2]bounding_boxo3d.geometry.AxisAlignedBoundingBox(min_boundmin_bound,max_boundmax_bound)roipoint_cloud.crop(bounding_box)roi.paint_uniform_color([1.00,0.00,0.00])# paint keypoints as red:keypoints_in_roikeypoints.loc[(((keypoints[x]min_bound[0])(keypoints[x]max_bound[0]))((keypoints[y]min_bound[1])(keypoints[y]max_bound[1]))((keypoints[z]min_bound[2])(keypoints[z]max_bound[2]))),:]np.asarray(point_cloud.colors)[keypoints_in_roi[id].values,:][1.0,0.0,0.0]o3d.visualization.draw_geometries([point_cloud])# describe keypoints:df_signature_visualization[]forkeypoint_idinkeypoints_in_roi[id].values:signaturedescribe(point_cloud,search_tree,keypoint_id,arguments.radius,6)df_pd.DataFrame.from_dict({index:np.arange(len(signature)),feature:signature})df_[keypoint_id]keypoint_id# f{keypoint_id:06d}df_signature_visualization.append(df_)df_signature_visualizationpd.concat(df_signature_visualization)# limit the number:df_signature_visualizationdf_signature_visualization.head(5)# draw the plot:plt.figure(numNone,figsize(16,9))sns.lineplot(xindex,yfeature,huekeypoint_id,stylekeypoint_id,markersTrue,dashesFalse,datadf_signature_visualization)plt.title(Description Visualization for Keypoints)plt.show()detection/iss.pyimportheapqimportnumpyasnpimportpandasaspdimportopen3daso3ddefdetect(point_cloud,search_tree,radius): Detect point cloud key points using Intrinsic Shape Signature(ISS) Parameters ---------- point_cloud: Open3D.geometry.PointCloud input point cloud search_tree: Open3D.geometry.KDTree point cloud search tree radius: float radius for ISS computing Returns ---------- point_cloud: numpy.ndarray Velodyne measurements as N-by-3 numpy ndarray # points handler:pointsnp.asarray(point_cloud.points)# keypoints container:keypoints{id:[],x:[],y:[],z:[],lambda_0:[],lambda_1:[],lambda_2:[]}# cache for number of radius nearest neighbors:num_rnn_cache{}# heapq for non-maximum suppression:pq[]foridx_center,centerinenumerate(points):# find radius nearest neighbors:[k,idx_neighbors,_]search_tree.search_radius_vector_3d(center,radius)# for each point get its nearest neighbors count:w[]deviation[]foridx_neighborinnp.asarray(idx_neighbors[1:]):# check cache:ifnotidx_neighborinnum_rnn_cache:[k_,_,_]search_tree.search_radius_vector_3d(points[idx_neighbor],radius)num_rnn_cache[idx_neighbor]k_# update:w.append(num_rnn_cache[idx_neighbor])deviation.append(points[idx_neighbor]-center)# calculate covariance matrix:wnp.asarray(w)deviationnp.asarray(deviation)cov(1.0/w.sum())*np.dot(deviation.T,np.dot(np.diag(w),deviation))# get eigenvalues:w,_np.linalg.eig(cov)ww[w.argsort()[::-1]]# add to pq:heapq.heappush(pq,(-w[2],idx_center))# add to dataframe:keypoints[id].append(idx_center)keypoints[x].append(center[0])keypoints[y].append(center[1])keypoints[z].append(center[2])keypoints[lambda_0].append(w[0])keypoints[lambda_1].append(w[1])keypoints[lambda_2].append(w[2])# non-maximum suppression:suppressedset()whilepq:_,idx_centerheapq.heappop(pq)ifnotidx_centerinsuppressed:# suppress its neighbors:[_,idx_neighbors,_]search_tree.search_radius_vector_3d(points[idx_center],radius)foridx_neighborinnp.asarray(idx_neighbors[1:]):suppressed.add(idx_neighbor)else:continue# format:keypointspd.DataFrame.from_dict(keypoints)# first apply non-maximum suppression:keypointskeypoints.loc[keypoints[id].apply(lambdaid:notidinsuppressed),keypoints.columns]# then apply decreasing ratio test:keypointskeypoints.loc[(keypoints[lambda_0]keypoints[lambda_1])(keypoints[lambda_1]keypoints[lambda_2]),keypoints.columns]returnkeypointsdescription/fpfh.pyimportnumpyasnpimportpandasaspdimportopen3daso3ddefget_spfh(point_cloud,search_tree,keypoint_id,radius,B): Describe the selected keypoint using Simplified Point Feature Histogram(SPFH) Parameters ---------- point_cloud: Open3D.geometry.PointCloud input point cloud search_tree: Open3D.geometry.KDTree point cloud search tree keypoint_id: ind keypoint index radius: float nearest neighborhood radius B: float number of bins for each dimension Returns ---------- # points handler:pointsnp.asarray(point_cloud.points)# get keypoint:keypointnp.asarray(point_cloud.points)[keypoint_id]# find radius nearest neighbors:[k,idx_neighbors,_]search_tree.search_radius_vector_3d(keypoint,radius)# remove query point:idx_neighborsidx_neighbors[1:]# get normalized diff:diffpoints[idx_neighbors]-keypoint diff/np.linalg.norm(diff,ord2,axis1)[:,None]# get n1:n1np.asarray(point_cloud.normals)[keypoint_id]# get u:un1# get v:vnp.cross(u,diff)# get w:wnp.cross(u,v)# get n2:n2np.asarray(point_cloud.normals)[idx_neighbors]# get alpha:alpha(v*n2).sum(axis1)# get phi:phi(u*diff).sum(axis1)# get theta:thetanp.arctan2((w*n2).sum(axis1),(u*n2).sum(axis1))# get alpha histogram:alpha_histogramnp.histogram(alpha,binsB,range(-1.0,1.0))[0]alpha_histogramalpha_histogram/alpha_histogram.sum()# get phi histogram:phi_histogramnp.histogram(phi,binsB,range(-1.0,1.0))[0]phi_histogramphi_histogram/phi_histogram.sum()# get theta histogram:theta_histogramnp.histogram(theta,binsB,range(-np.pi,np.pi))[0]theta_histogramtheta_histogram/theta_histogram.sum()# build signature:signaturenp.hstack((# alpha:alpha_histogram,# phi:phi_histogram,# theta:phi_histogram))returnsignaturedefdescribe(point_cloud,search_tree,keypoint_id,radius,B): Describe the selected keypoint using Fast Point Feature Histogram(FPFH) Parameters ---------- point_cloud: Open3D.geometry.PointCloud input point cloud search_tree: Open3D.geometry.KDTree point cloud search tree keypoint_id: ind keypoint index radius: float nearest neighborhood radius B: float number of bins for each dimension Returns ---------- # points handler:pointsnp.asarray(point_cloud.points)# get keypoint:keypointnp.asarray(point_cloud.points)[keypoint_id]# find radius nearest neighbors:[k,idx_neighbors,_]search_tree.search_radius_vector_3d(keypoint,radius)ifk1:returnNone# remove query point:idx_neighborsidx_neighbors[1:]# weights:w1.0/np.linalg.norm(points[idx_neighbors]-keypoint,ord2,axis1)# SPFH from neighbors:Xnp.asarray([get_spfh(point_cloud,search_tree,i,radius,B)foriinidx_neighbors])# neighborhood contribution:spfh_neighborhood1.0/(k-1)*np.dot(w,X)# query point spfh:spfh_queryget_spfh(point_cloud,search_tree,keypoint_id,radius,B)# finally:spfhspfh_queryspfh_neighborhood# normalize again:spfhspfh/np.linalg.norm(spfh)returnspfh这份代码中的特征点检测与特征点描述这段代码实现了一个完整的点云特征提取流程先检测关键点再对关键点做描述。它由三个部分组合而成iss.py特征点检测采用 ISSIntrinsic Shape Signaturefpfh.py特征描述采用 FPFHFast Point Feature Histogrammain.py整体流程控制、可视化与说明1. 特征点检测ISS输入检测模块的输入是带法线的点云点云的 KD-tree半径参数radius检测思想ISS 的核心思想是从点云中寻找那些局部几何变化最显著的点通常是“角点”或“曲率高”的位置。主要处理步骤对每个点做半径邻域搜索找到这个点周围一定范围内的邻居。对邻居中的每个点再分别统计它自己的半径邻居数量作为权重。这一步目的是让更“稠密”区域的邻居对局部描述贡献更大。计算目标点邻域的加权协方差矩阵反映邻域点云的空间分布。求协方差矩阵的特征值并按降序排序lambda_0是最大特征值lambda_1是中间特征值lambda_2是最小特征值用lambda_2作为点的兴趣度指标。因为lambda_2代表最小方向上的变化程度越大说明该点邻域不是平面/边而是更“角点式”的结构。非极大值抑制先把所有点按兴趣度排序构成一个优先队列。依次取出当前最大的点保留它并把它半径范围内的所有邻居都标记为“已抑制”。这样可以避免在同一局部区域里保留多个相近关键点。Ratio Test在非极大值抑制后还做一个特征值顺序筛选只保留满足lambda_0 lambda_1 lambda_2的点这一步的目的在于进一步排除那些邻域结构更像平面或线性的点保留更具“角点特征”的点。输出最终输出的是关键点集合每个关键点包含点索引坐标对应的三个特征值2. 特征点描述FPFH / SPFH描述目的描述器的目标是为每个关键点生成一个向量能够表达其局部几何和法线关系便于后续的匹配或分类。基本单元SPFHSPFH 是 FPFH 的基础是对一个点邻域中几何关系的直方图描述。计算过程对关键点的邻域点做半径搜索。对每个邻居计算关键点到邻居的相对方向向量并归一化。用关键点的法线和邻居的法线构建局部参考系用关键点法线作为uv和w由叉乘得到计算三个角度量alpha邻居法线在v方向上的投影phi关键点法线与相对方向的投影theta邻居法线相对于关键点法线的方位角对这三个角度分别构建直方图并做归一化。最终 SPFH把alpha、phi、theta的直方图拼接起来得到一个描述向量这个向量反映了关键点邻域内法线方向与空间分布的统计特征快速版本FPFHFPFH 的想法是让描述更稳定、更具上下文感先计算关键点自身的 SPFH再计算邻居每个点的 SPFH用距离加权方式把邻居 SPFH 的贡献聚合起来把关键点自身的 SPFH 与邻居加权 SPFH 相加对结果做归一化这样得到的描述向量既包含关键点自身邻域的信息也融合了邻居点的邻域信息因此比单纯的 SPFH 更鲁棒。3. 代码中的实际流程在main.py中整体流程是读取带法线的点云建立 KD-tree调用 ISS 检测器得到关键点在可视化中将关键点标成红色对 ROI 内的关键点计算其 FPFH 描述最后把若干关键点描述曲线绘制出来注意点关键点检测依赖radius参数半径大小直接影响邻域范围进而影响检测结果。描述也依赖同一radius因此整个流程对尺度敏感。FPFH计算过程中用的是 “距离倒数” 作为邻居权重距离越近的邻居对描述贡献越大。4. 关键算法区别与作用ISS 特征点检测作用从点云中挑选出“显著”、“稳定”的点原理基于邻域协方差矩阵特征值找几何结构显著点结果只产生一组关键点不需要提前知道点数FPFH 特征描述作用为关键点生成具有几何意义的向量描述原理统计法线之间的角度关系构建直方图结果得到一个向量可用于匹配、检索、分类5. 额外说明这段代码里ISS 检测完全不直接依赖点云法线而是依赖邻域几何分布但点云本身带的 normals 对 FPFH 描述非常重要。main.py的可视化部分还会对一个 ROI 区域做红色高亮用来观察局部关键点。如果你想我还可以继续帮你把“ISS 和 SIFT/ Harris 关键点”或“FPFH 和 SHOT 描述子”之间的差异做对比解释。