三维点云处理（二十一）：PFH 与 FPFH 点特征直方图描述子

在完成特征点检测（Keypoint Detection）后，下一步是提取局部特征描述子（Feature Descriptor）。描述子用高维向量刻画特征点邻域的几何构型，是点云配准（Registration）、三维姿态估计和语义分割的基石。

本文深入解析最经典的两大局部特征描述子：PFH（Point Feature Histogram） 与其加速算法 FPFH（Fast Point Feature Histogram）。

---

一、特征描述子引言与应用背景

1.1 特征描述子的核心任务

在三维点云中，特征提取分为两个独立但紧密联系的步骤：

1. 检测（Detection）：在点云中识别出具有显著几何特征的关关键点（Keypoints，如 Harris 3D/6D、ISS）。
2. 描述（Description）：以每个关键点为中心，在其局部邻域内提取一个能够表征几何形状的特征向量（Descriptor）。
3. 匹配（Matching）：通过计算不同视角点云描述子之间的距离（如欧氏距离），寻找点与点之间的对应关系（Correspondences）。

1.2 理想描述子的四大属性

一个优秀的三维局部特征描述子必须具备以下性质：

属性	物理含义	违背反例
6D 位姿不变性	无论点云如何旋转或平移，同一位置的描述子保持不变	直接利用绝对世界坐标
采样密度鲁棒性	点云稀疏度或扫描距离发生变化时，描述子保持稳定	强依赖固定点数或绝对距离的统计
抗噪与遮挡能力	传感器噪声或局部遮挡时不发生剧烈变动	直接使用二阶高阶微分或原始点阵
高区分度（Distinctiveness）	不同几何结构（平面、角点、球极）对应截然不同的向量	维度过低或丢弃空间结构信息

---

二、PFH（Point Feature Histogram）数学原理

PFH（点特征直方图） 是 Rusu 等人于 2008 年提出的经典几何描述子。其核心思想是：利用局部邻域内**所有点对（Pairwise）**的法向量相对偏转角，构建高维概率统计直方图。

2.1 局部 Darboux 坐标系（UVW 框架）

为了保证刚体变换不变性（6D Pose Invariance），PFH 不直接使用全局坐标，而是为邻域内的每一对点 建立一个随局部法向量自适应旋转的 Darboux 局部参考坐标系。

假定 为源点（Source Point）， 为目标点（Target Point），且 分别为两点处的单位法向量。选择法向量与两点连线夹角较小的一点作为 （确保坐标系确定），Darboux 坐标系的正交基 定义为：

2.2 四元组角度特征（Quadruplet）

在 Darboux 框架下，点对 及其法向量 的相对几何姿态可完全由 3 个角度和 1 个距离构成的四元组 唯一表征：

  四元组 [α, φ, θ, d] 的几何意义
  
  α (alpha):  目标法向量 n_t 在 v-w 平面上的投影 → 描述两切平面的"面外"倾斜程度
  φ (phi):    连线 (p_t - p_s) 与 u 轴 (n_s) 的夹角 → 描述两点间的相对"高度与倾角"
  θ (theta):  目标法向量 n_t 在 u-w 平面上的偏转角 → 描述绕主法向的旋转角度
  d:          两点之间的三维欧氏距离

注意：在实践中，通常忽略距离参数 。因为点云采样密度在不同视角下会剧烈改变，忽略 可以使描述子对点云密度变化和缩放具备强鲁棒性。

2.3 PFH 的直方图离散化与维度计算

对关键点 及其半径 支撑邻域内包含的 个邻域点：

1. 组合全连接点对：在 个邻域点中取遍所有点对，共计产生 对组合。
2. 计算角度三元组：对每对组合计算 。
3. 三维体素网格（Voxel Grid）直方图统计：
   • 将每个角度的取值范围区间均匀划分为 个区间（Bin，通常取 ）。
   • 三个角度组合形成一个大小为 的三维直方图网格。
   • 对每一个点对的三元组，将其投射到对应的体素网格中并累加频次。
4. 归一化：将长度为 的直方图向量进行 或 归一化，即得到 的 PFH 描述子（当 时，维度为 维）。

2.4 PFH 的时间复杂度与局限

对于包含 个点的点云，若每个点的半径邻域平均包含 个点：

• 单个点的计算复杂度为 ；
• 整幅点云的计算复杂度为 。

当点云密度较高（如 ）或点云规模达百万级时， 的二次方复杂度会导致计算极度缓慢，难以满足实时处理的要求。

---

三、FPFH（Fast Point Feature Histogram）加速机制

为了解决 PFH 计算复杂度过高的痛点，Rusu 等人在 2009 年提出了 FPFH（快速点特征直方图），将计算复杂度从 成功降低至 ，同时保持了与 PFH 相当的几何区分度。

3.1 简化版点特征直方图（SPFH）

FPFH 不再计算邻域内所有点对的两两关系，而是引入了简化版 PFH（Simplified Point Feature Histogram, SPFH）：

1. 对于查询关键点 ，仅计算它与邻域内 个邻域点 之间直接相连的单向点对。
2. 对这 个点对计算角度三元组 。
3. 分别为 建立 1D 直方图（设 Bin 数为 ，如 ）。
4. 将 3 个 1D 直方图首尾拼接，得到该点的 SPFH 向量（维度为 维）。

3.2 邻域二次加权聚合

仅靠 SPFH 丢失了邻域点之间的相互几何关联。为了弥补这一信息，FPFH 通过二次邻域加权聚合，将邻域点的 SPFH 混合到中心点中：

权重 定义为中心点 与邻域点 之间三维欧氏距离倒数的反比（距离反比加权）：

3.3 拓扑覆盖范围扩展与跨阶交互

通过二次加权聚合，FPFH 产生了独特的几何效应：

• 空间拓扑：计算从星形连接（Star Topology）扩展到了半径为 的跨阶邻域（Range ）。
• 二次统计信息：使得距离为 范围内的点能够通过中间邻域点间接地向中心点传递相对几何姿态。

---

四、PFH vs FPFH 详尽对比

下面的对比表总结了 PFH 与 FPFH 的核心差异：

维度	PFH (Point Feature Histogram)	FPFH (Fast Point Feature Histogram)
拓扑连接方式	全连接网格 (Fully Connected Topology)	星形连接 + 邻域加权聚合 (Star Topology + Aggregation)
计算复杂度	（二次方量级，极慢）	（线性量级，极快）
直方图编码方式	3D 联合体素直方图 ()	3 个独立 1D 直方图拼接 ()
典型特征维度	125 维 ()	33 维 ()
邻域影响半径	支撑半径	扩展支撑半径
重复边计算	每条边无重复计算	部分邻域节点之间的连接被重复计算（Thick Edges）
适用场景	小规模高精度语义分类/局部识别	大规模实时点云配准 (Registration / SLAM)

4.1 为什么 FPFH 是全局粗配准的“工业界黄金标准”？

在三维视觉工程落地中（如基于 RANSAC 或 Fast Global Registration 的点云粗配准），FPFH 几乎成为了各大开源库（PCL, Open3D）的默认描述子，原因在于它在以下三个维度的完美平衡：

1. 极致的并行化潜力 (GPU/多线程友好)： 计算 SPFH 时，每个点与其 个邻居的计算是完全独立的；在第二步邻域聚合时，依然只是单纯的查表与加权求和。这种完全解耦的两阶段（Two-pass）算法极其适合在 CPU 上通过 OpenMP 进行多线程加速，或在 GPU（CUDA）上进行成千上万个核心的高效并行映射。
2. 紧凑的维度与极速的特征空间检索： 描述子匹配需要在特征空间中利用 KD-Tree 进行 K-NN 检索。33 维的 FPFH 相比于 125 维的 PFH 或高达 352 维的 SHOT，大大缓解了高维空间搜索的“维数灾难”。它在确保几何区分度（足以区分边缘、角点与平面）的同时，将特征匹配的速度提升了一个数量级。
3. 隐式的多尺度特征融合： 通过第二步的邻域二次加权，FPFH 实际上隐式地将感受野从 扩展到了 。这种“中心清晰、边缘模糊”的描述模式天然具备了对抗遮挡和传感器噪声的强鲁棒性。

---

五、Python 代码实现与 Open3D 实践

5.1 纯 NumPy 实现 FPFH (原理演示)

import numpy as np
from scipy.spatial import KDTree


def compute_darboux_triplet(p_s, n_s, p_t, n_t):
    """
    计算两点及其法向量在 Darboux 框架下的角度三元组 [alpha, phi, theta]。
    """
    dp = p_t - p_s
    d = np.linalg.norm(dp)
    if d < 1e-10:
        return 0.0, 0.0, 0.0

    # 1. 构建 UVW 框架
    u = n_s
    v = np.cross(u, dp / d)
    v_norm = np.linalg.norm(v)
    if v_norm < 1e-10:
        return 0.0, 0.0, 0.0
    v /= v_norm
    w = np.cross(u, v)

    # 2. 计算三角度特征
    alpha = np.dot(v, n_t)
    phi = np.dot(u, dp / d)
    theta = np.arctan2(np.dot(w, n_t), np.dot(u, n_t))

    return alpha, phi, theta


def compute_spfh(points, normals, query_idx, tree, radius, num_bins=11):
    """
    计算单个点的 SPFH 向量 (33维)。
    """
    p_q = points[query_idx]
    n_q = normals[query_idx]

    neighbor_indices = tree.query_ball_point(p_q, radius)
    if len(neighbor_indices) <= 1:
        return np.zeros(3 * num_bins), neighbor_indices

    alphas, phis, thetas = [], [], []

    for idx in neighbor_indices:
        if idx == query_idx:
            continue
        p_k = points[idx]
        n_k = normals[idx]

        # 保证 u 轴选择夹角小的一方
        if np.dot(n_q, p_k - p_q) >= 0:
            a, p, t = compute_darboux_triplet(p_q, n_q, p_k, n_k)
        else:
            a, p, t = compute_darboux_triplet(p_k, n_k, p_q, n_q)

        alphas.append(a)
        phis.append(p)
        thetas.append(t)

    # 建立 1D 直方图
    hist_alpha, _ = np.histogram(alphas, bins=num_bins, range=(-1.0, 1.0))
    hist_phi, _ = np.histogram(phis, bins=num_bins, range=(-1.0, 1.0))
    hist_theta, _ = np.histogram(thetas, bins=num_bins, range=(-np.pi, np.pi))

    spfh = np.concatenate([hist_alpha, hist_phi, hist_theta]).astype(float)
    sum_val = np.sum(spfh)
    if sum_val > 0:
        spfh /= sum_val

    return spfh, neighbor_indices


def compute_fpfh_numpy(points, normals, radius=0.05, num_bins=11):
    """
    纯 NumPy 计算点云所有点的 FPFH 描述子。
    :return: (N, 33) FPFH 特征矩阵
    """
    N = points.shape[0]
    tree = KDTree(points)

    # 步骤 1: 计算所有点的 SPFH
    spfh_array = np.zeros((N, 3 * num_bins))
    neighbors_list = []

    for i in range(N):
        spfh, neighbors = compute_spfh(points, normals, i, tree, radius, num_bins)
        spfh_array[i] = spfh
        neighbors_list.append(neighbors)

    # 步骤 2: 邻域二次加权求和
    fpfh_array = np.zeros((N, 3 * num_bins))

    for i in range(N):
        neighbors = neighbors_list[i]
        if len(neighbors) <= 1:
            fpfh_array[i] = spfh_array[i]
            continue

        p_q = points[i]
        weights = []
        spfh_neighbors = []

        for k in neighbors:
            if k == i:
                continue
            dist = np.linalg.norm(p_q - points[k])
            if dist > 1e-10:
                weights.append(1.0 / dist)
                spfh_neighbors.append(spfh_array[k])

        if len(weights) == 0:
            fpfh_array[i] = spfh_array[i]
        else:
            weights = np.array(weights)
            weights /= np.sum(weights)
            weighted_spfh = np.sum(weights[:, None] * np.array(spfh_neighbors), axis=0)
            fpfh_array[i] = spfh_array[i] + weighted_spfh
            
            # L1 归一化
            norm_val = np.sum(fpfh_array[i])
            if norm_val > 0:
                fpfh_array[i] /= norm_val

    return fpfh_array

5.2 Open3D 内置 FPFH API 调用

在实际工程项目中，使用 Open3D 提供的 C++ 加速 API 进行 FPFH 计算：

import open3d as o3d
import numpy as np


def extract_fpfh_open3d(pcd, voxel_size=0.05):
    """
    使用 Open3D 高性能提取 FPFH 描述子。
    """
    # 1. 法向量估计
    radius_normal = voxel_size * 2
    pcd.estimate_normals(
        o3d.geometry.KDTreeSearchParamHybrid(radius=radius_normal, max_nn=30)
    )

    # 2. 计算 FPFH 描述子
    radius_feature = voxel_size * 5
    fpfh = o3d.pipelines.registration.compute_fpfh_feature(
        pcd,
        o3d.geometry.KDTreeSearchParamHybrid(radius=radius_feature, max_nn=100)
    )

    # fpfh.data 的形状为 (33, N)
    features = fpfh.data.T
    print(f"提取点数: {len(pcd.points)}, FPFH 特征矩阵形状: {features.shape}")
    return features

---

六、总结与知识脉络

概念	核心要点
Darboux 坐标系	基于法向量 和两点连线建立的局部正交基 ，消除 6D 位姿依赖
四元组特征	表征相对几何姿态；忽略 可获得采样密度鲁棒性
PFH 核心	全连接点对 + 维三维体素直方图，表征能力强但计算极慢
SPFH	仅计算中心点到邻域点星形连接的 1D 直方图拼接（ 维）
FPFH 核心	SPFH + 邻域距离反比二次加权，将复杂度降至 ，支撑半径扩展至

PFH 与 FPFH 完美解答了“如何高效统计局部几何表面曲率变异”的问题。在实际工程与算法体系中，FPFH 结合 RANSAC 被广泛应用于粗配准（Coarse Registration）。下一章将探讨另一种经典描述子——SHOT（Signature of Histograms of OrienTations）。