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

一、特征点的描述

1) PFH

PFH方法核心特点：

• 实现旋转不变性
• 捕捉局部表面几何变化
• 基于点对法向量与相对位置计算参数
• 复杂度缺陷：n个特征点需计算nk²次运算
• 敏感性问题：结果严重依赖法向量精度
  

2) FPFH

FPFH改进方案：

• 简化计算：每个点仅计算与直接邻居的连线（k次运算）

• 直方图构建：采用3b维拼接方案（非三维投票）

• 加权融合 ：

  • 中心点SPFH（简化直方图）
  • 邻居点SPFH的加权平均
  • 复杂度优化：从nk²降至nk
    

4) fph与pfh的区别

二. 代码实践

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:parser=argparse.ArgumentParser("Detect ISS keypoints on ModelNet40 dataset.")# add required and optional groups:required=parser.add_argument_group('Required')optional=parser.add_argument_group('Optional')# add required:required.add_argument("-i",dest="input",help="Input path of ModelNet40 sample.",required=True)required.add_argument("-r",dest="radius",help="Radius for radius nearest neighbor definition.",required=True,type=float)# parse arguments:returnparser.parse_args()if__name__=='__main__':# parse arguments:arguments=get_arguments()# load point cloud:point_cloud=read_modelnet40_normal(arguments.input)# compute surface normals:# point_cloud.estimate_normals(# search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.1, max_nn=30)# )# build search tree:search_tree=o3d.geometry.KDTreeFlann(point_cloud)# detect keypoints:keypoints=detect(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_bound=point_cloud.get_max_bound()min_bound=point_cloud.get_min_bound()center=(min_bound+max_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_box=o3d.geometry.AxisAlignedBoundingBox(min_bound=min_bound,max_bound=max_bound)roi=point_cloud.crop(bounding_box)roi.paint_uniform_color([1.00,0.00,0.00])# paint keypoints as red:keypoints_in_roi=keypoints.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:signature=describe(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_visualization=pd.concat(df_signature_visualization)# limit the number:df_signature_visualization=df_signature_visualization.head(5)# draw the plot:plt.figure(num=None,figsize=(16,9))sns.lineplot(x="index",y="feature",hue="keypoint_id",style="keypoint_id",markers=True,dashes=False,data=df_signature_visualization)plt.title('Description Visualization for Keypoints')plt.show()

detection/iss.py

importheapqimportnumpyasnpimportpandasaspdimportopen3daso3ddefdetect(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:points=np.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:w=np.asarray(w)deviation=np.asarray(deviation)cov=(1.0/w.sum())*np.dot(deviation.T,np.dot(np.diag(w),deviation))# get eigenvalues:w,_=np.linalg.eig(cov)w=w[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:suppressed=set()whilepq:_,idx_center=heapq.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:keypoints=pd.DataFrame.from_dict(keypoints)# first apply non-maximum suppression:keypoints=keypoints.loc[keypoints['id'].apply(lambdaid:notidinsuppressed),keypoints.columns]# then apply decreasing ratio test:keypoints=keypoints.loc[(keypoints['lambda_0']>keypoints['lambda_1'])&(keypoints['lambda_1']>keypoints['lambda_2']),keypoints.columns]returnkeypoints

description/fpfh.py

importnumpyasnpimportpandasaspdimportopen3daso3ddefget_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:points=np.asarray(point_cloud.points)# get keypoint:keypoint=np.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_neighbors=idx_neighbors[1:]# get normalized diff:diff=points[idx_neighbors]-keypoint diff/=np.linalg.norm(diff,ord=2,axis=1)[:,None]# get n1:n1=np.asarray(point_cloud.normals)[keypoint_id]# get u:u=n1# get v:v=np.cross(u,diff)# get w:w=np.cross(u,v)# get n2:n2=np.asarray(point_cloud.normals)[idx_neighbors]# get alpha:alpha=(v*n2).sum(axis=1)# get phi:phi=(u*diff).sum(axis=1)# get theta:theta=np.arctan2((w*n2).sum(axis=1),(u*n2).sum(axis=1))# get alpha histogram:alpha_histogram=np.histogram(alpha,bins=B,range=(-1.0,+1.0))[0]alpha_histogram=alpha_histogram/alpha_histogram.sum()# get phi histogram:phi_histogram=np.histogram(phi,bins=B,range=(-1.0,+1.0))[0]phi_histogram=phi_histogram/phi_histogram.sum()# get theta histogram:theta_histogram=np.histogram(theta,bins=B,range=(-np.pi,+np.pi))[0]theta_histogram=theta_histogram/theta_histogram.sum()# build signature:signature=np.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:points=np.asarray(point_cloud.points)# get keypoint:keypoint=np.asarray(point_cloud.points)[keypoint_id]# find radius nearest neighbors:[k,idx_neighbors,_]=search_tree.search_radius_vector_3d(keypoint,radius)ifk<=1:returnNone# remove query point:idx_neighbors=idx_neighbors[1:]# weights:w=1.0/np.linalg.norm(points[idx_neighbors]-keypoint,ord=2,axis=1)# SPFH from neighbors:X=np.asarray([get_spfh(point_cloud,search_tree,i,radius,B)foriinidx_neighbors])# neighborhood contribution:spfh_neighborhood=1.0/(k-1)*np.dot(w,X)# query point spfh:spfh_query=get_spfh(point_cloud,search_tree,keypoint_id,radius,B)# finally:spfh=spfh_query+spfh_neighborhood# normalize again:spfh=spfh/np.linalg.norm(spfh)returnspfh

这份代码中的特征点检测与特征点描述

这段代码实现了一个完整的点云特征提取流程：先检测关键点，再对关键点做描述。它由三个部分组合而成：

• iss.py：特征点检测，采用 ISS（Intrinsic Shape Signature）
• fpfh.py：特征描述，采用 FPFH（Fast Point Feature Histogram）
• main.py：整体流程控制、可视化与说明

1. 特征点检测：ISS

输入

检测模块的输入是：

• 带法线的点云
• 点云的 KD-tree
• 半径参数radius

检测思想

ISS 的核心思想是：从点云中寻找那些局部几何变化最显著的点，通常是“角点”或“曲率高”的位置。

主要处理步骤

1. 对每个点做半径邻域搜索，找到这个点周围一定范围内的邻居。
2. 对邻居中的每个点，再分别统计它自己的半径邻居数量，作为权重。这一步目的是让更“稠密”区域的邻居对局部描述贡献更大。
3. 计算目标点邻域的加权协方差矩阵，反映邻域点云的空间分布。
4. 求协方差矩阵的特征值，并按降序排序：
   • lambda_0是最大特征值
   • lambda_1是中间特征值
   • lambda_2是最小特征值
5. 用lambda_2作为点的兴趣度指标。因为lambda_2代表最小方向上的变化程度，越大说明该点邻域不是平面/边，而是更“角点式”的结构。

非极大值抑制

• 先把所有点按兴趣度排序，构成一个优先队列。
• 依次取出当前最大的点，保留它，并把它半径范围内的所有邻居都标记为“已抑制”。
• 这样可以避免在同一局部区域里保留多个相近关键点。

Ratio Test

• 在非极大值抑制后，还做一个特征值顺序筛选：
  • 只保留满足lambda_0 > lambda_1 > lambda_2的点
• 这一步的目的在于进一步排除那些邻域结构更像平面或线性的点，保留更具“角点特征”的点。

输出

最终输出的是关键点集合，每个关键点包含：

• 点索引
• 坐标
• 对应的三个特征值

2. 特征点描述：FPFH / SPFH

描述目的

描述器的目标是为每个关键点生成一个向量，能够表达其局部几何和法线关系，便于后续的匹配或分类。

基本单元：SPFH

SPFH 是 FPFH 的基础，是对一个点邻域中几何关系的直方图描述。

计算过程

1. 对关键点的邻域点做半径搜索。
2. 对每个邻居，计算关键点到邻居的相对方向向量，并归一化。
3. 用关键点的法线和邻居的法线构建局部参考系：
   • 用关键点法线作为u
   • v和w由叉乘得到
4. 计算三个角度量：
   • alpha：邻居法线在v方向上的投影
   • phi：关键点法线与相对方向的投影
   • theta：邻居法线相对于关键点法线的方位角
5. 对这三个角度分别构建直方图，并做归一化。

最终 SPFH

• 把alpha、phi、theta的直方图拼接起来，得到一个描述向量
• 这个向量反映了关键点邻域内法线方向与空间分布的统计特征

快速版本：FPFH

FPFH 的想法是让描述更稳定、更具上下文感：

1. 先计算关键点自身的 SPFH
2. 再计算邻居每个点的 SPFH
3. 用距离加权方式，把邻居 SPFH 的贡献聚合起来
4. 把关键点自身的 SPFH 与邻居加权 SPFH 相加
5. 对结果做归一化

这样得到的描述向量既包含关键点自身邻域的信息，也融合了邻居点的邻域信息，因此比单纯的 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 描述子”之间的差异做对比解释。