三维点云处理（十九）：Harris 角点检测 (2D/3D/6D)

特征点检测（Keypoint Detection） 是三维点云匹配、三维姿态估计与点云配准的第一步。通过筛选具有剧烈几何变化与高可重复性的少数关键点，能够大幅降低高维特征描述子的计算开销。

在二维图像处理中，Harris 角点检测器是最为经典的特征提取算法。本文深入探讨 Harris 算法从 2D 图像向 3D 点云及多维空间（3D 几何 + 反射强度/颜色/法线）推广的数学推导、算法实现与原理细节。

核心算法特点对比

算法方案	核心原理	筛选判定准则	优点与适用场景
Harris 2D	图像梯度结构张量		2D 像素网格，利用梯度检测角点
Harris 3D (曲面拟合法)	局部拟合切曲面导数	拟合曲面高程梯度的结构张量角点响应	借助切平面局部坐标系，几何意义明确
Harris 3D (协方差谱法)	3D 坐标局部协方差特征值	局部极大值	计算高效，检测三个正交方向变异均显著的三维角点
Harris 6D (联合特征)	拼接坐标与法线/强度构建联合协方差	6D 结构张量角点响应极值	可突破几何限制，在平面上跨维度检测颜色/强度或法向突变点

一、特征点的三大准则与几何分类

1.1 什么是好的特征点？

一个优秀的三维特征点（Keypoint）应当满足以下三大准则：

准则	物理含义	违背反例
可重复性（Repeatability）	同一物理场景在不同视角、扫描距离或噪声扰动下，能反复检测到相同的特征点	只在特定拍摄角度可检测出的孤立点
显著性（Saliency）	点的局部邻域具有与众不同的几何结构，各向异性强	大面积平坦区域上的平庸点
信息量（Informativeness）	包含足够丰富的几何/纹理信息，可用于后续特征描述与匹配	可稳定检测但缺乏几何判别力的点

1.2 局部几何形态与特征谱分类

从三维曲面的局部几何形态来看，表面上的点可分为三类。通过 PCA（协方差矩阵分析）导出的特征值谱（降序），可清晰地区分它们的几何属性：

点类型	特征值分布（降序）	几何物理表达	Harris 响应趋势
平面点 (Planar)		点集中在二维切平面内，法线方向变异极小	响应极低
边缘点 (Edge)		点集中在一维直线上	响应中等（仅一个维度变化）
3D 角点 (Corner)	且均显著	三维空间三个正交方向上均有剧烈几何突变	极高（三维角点）

二、Harris 2D 角点检测原理（回顾）

在 2D 灰度图像中，Harris 通过移动小窗口观察像素灰度的变化。在窗口发生微小位移时，自相关函数表达式为：

其中为 2D 结构张量（二阶矩矩阵）：

的两个特征值决定了局部区域在两个正交方向上的梯度强度：

平坦区域：均很小；
边缘区域：或；
角点区域：均很大。

为避免显式计算特征分解，Harris 定义角点响应函数：

三、Harris 3D：从 2D 到 3D 点云的推广

在无规则、无网格结构的三维点云中，将 2D Harris 推广到 3D 主要有两条经典路径：

3.1 路径一：曲面拟合与梯度张量法（Sipiran 路径）

确定局部切平面：对点及其半径支撑邻域内的点云，计算主法向量，以为原点、法线为轴正方向，构建局部坐标系。
拟合曲面：将邻域点投影到局部坐标系中，拟合二次曲面：
计算高程梯度：在原点处计算二次曲面的偏导数 , 。
构造二阶矩张量 ：
计算角点响应：同 2D Harris 计算响应。

3.2 路径二：3D 协方差矩阵特征谱法

在三维空间中，直接利用点支撑邻域内三维坐标的协方差矩阵：

设的三个特征值为。可以通过特征值乘积与和的运算构造三维 Harris 响应函数：

或者采用标准化更强、不依赖量纲参数的表达形式：

物理直觉：只有当三个特征值都显著大于 0 时（即点云在空间三个正交方向上都有剧烈变动），响应值才会取得极大值，精确对应三维尖锐角点（如立方体顶角、三平面交汇点）。

3.3 Harris 3D 参数选择的工程经验

在实际应用中，提取的三维角点质量极大程度上取决于以下两个超参数：

邻域搜索半径 (Scale / Radius)：
- 物理含义：决定了角点检测的感受野（尺度）。
- 调参法则：如果设得过小（接近点间距），提取的将全是点云的高频传感器噪声（把微小的粗糙表面当成角点）；如果设得过大，则会严重模糊几何细节，甚至多个相邻的小角点会被融合成一个无意义的大团。通常建议设置为点云平均点间距（Resolution）的 5~15 倍。
Harris 经验系数 ：
- 物理含义：控制着角点筛选的“严苛程度”。
- 调参法则：通常取值。越大，惩罚项越大，筛选条件越苛刻，提取到的角点数量越少但绝对尖锐度越高。若想提取更多潜在的关键点，可适当降低值。

四、Harris 6D 与强度扩展（Harris with Intensity / Normals）

4.1 几何角点的局限性

纯 3D 坐标算出的 Harris 3D 只能检测几何结构上剧烈交叠的“几何角点”。但在实际应用中（如搭载 RGB-D 传感器或带反射强度的 LiDAR 点云）：

如果一个物体表面在几何上是完全平坦的，但表面上有丰富的颜色图案纹理或反射强度变化（如地面白线、标志牌），Harris 3D 会完全失效；
如果曲面法向量变化剧烈（如折线或凸起边界），坐标分布可能变异不大，但法线方向有巨大突变。

4.2 Harris 6D 联合协方差

为了克服这一缺陷，经典文献中提出了扩展的 Harris 6D（或带 intensities 的 Harris 3D）：

将点的 3D 空间坐标与 3D 法向量（或色彩/强度梯度）进行拼接，构成 6 维联合特征向量：

在邻域内计算联合协方差矩阵：

计算特征分解的最小特征值或综合角点响应。

text

  Harris 3D vs Harris 6D 对比
  
  特征维度            检测类型                       典型失效例
  ─────────────────────────────────────────────────────────────────────────────
  Harris 3D (3D)      纯几何 3D 顶角与三面交点         平坦表面上的颜色/强度/法线突变点
  Harris 6D (6D)      几何角点 + 纹理/强度/法向突变点    全各向同性均匀散点

五、纯 Python 与 Open3D 代码实现

以下提供使用纯 Python (NumPy) 与 Open3D 实现 Harris 3D/6D 角点检测的完整代码。

python

import numpy as np
import open3d as o3d
from scipy.spatial import KDTree


def harris_3d_keypoints(pcd, radius=0.05, harris_k=0.04, top_k=500):
    """
    基于 3D 协方差矩阵特征谱法的 Harris 3D 角点检测实现
    
    :param pcd: Open3D PointCloud 对象
    :param radius: 邻域搜索半径
    :param harris_k: Harris 经验系数 k
    :param top_k: 保留响应值最大的前 top_k 个关键点
    :return: 关键点 open3d.geometry.PointCloud
    """
    points = np.asarray(pcd.points)
    N = len(points)
    tree = KDTree(points)
    
    responses = np.zeros(N)
    
    for i in range(N):
        # 搜索半径 r 内的所有邻域点
        idx = tree.query_ball_point(points[i], r=radius)
        if len(idx) < 5:
            continue
            
        neighbors = points[idx]
        # 计算 3x3 协方差矩阵
        centroid = np.mean(neighbors, axis=0)
        centered = neighbors - centroid
        cov = (centered.T @ centered) / len(neighbors)
        
        # 计算特征值 (降序 λ1 >= λ2 >= λ3)
        eigvals = np.linalg.eigvalsh(cov)
        l3, l2, l1 = eigvals[0], eigvals[1], eigvals[2]
        
        if l1 <= 1e-8:
            continue
            
        # 计算响应函数 R = λ1*λ2*λ3 / (λ1 + λ2 + λ3) 或 det(Σ) - k*tr(Σ)^3
        det_cov = l1 * l2 * l3
        trace_cov = l1 + l2 + l3
        responses[i] = det_cov / (trace_cov + 1e-8)
        
    # 非极大值抑制 (NMS) 保留 Top-K
    sorted_indices = np.argsort(responses)[::-1]
    selected_indices = []
    
    selected_mask = np.zeros(N, dtype=bool)
    for idx in sorted_indices:
        if responses[idx] <= 0:
            break
        if selected_mask[idx]:
            continue
            
        selected_indices.append(idx)
        if len(selected_indices) >= top_k:
            break
            
        # 抑制半径 r 内的其他点
        suppress_idx = tree.query_ball_point(points[idx], r=radius)
        selected_mask[suppress_idx] = True
        
    keypoints = pcd.select_by_index(selected_indices)
    return keypoints


if __name__ == "__main__":
    # 创建简单测试立方体点云
    mesh = o3d.geometry.TriangleMesh.create_box(width=1.0, height=1.0, depth=1.0)
    pcd = mesh.sample_points_uniformly(number_of_points=5000)
    
    keypoints = harris_3d_keypoints(pcd, radius=0.1, top_k=200)
    print(f"原始点数: {len(pcd.points)}, 提取 Harris 3D 关键点数: {len(keypoints.points)}")

六、总结

维度	Harris 2D	Harris 3D	Harris 6D
输入空间	2D 像素矩阵	3D 空间点集	6D 联合空间 / 强度
核心矩阵	梯度结构张量	空间坐标协方差	联合特征协方差
检测目标	图像角点与边缘突变	3D 空间立体几何角点	几何角点 + 纹理/反射强度/法线突变点
局限性	缺乏深度与 3D 空间不变性	对点云密度变化与尺度缩放敏感	计算开销较大，且依赖高质量法线/强度信息

Harris 角点检测器奠定了三维特征点检测的理论基础。为了进一步解决 Harris 依赖绝对尺度、对点云密度变化敏感的问题，在下一章中，我们将深入解析工业界更为通用的关键点检测算法——[20-iss-keypoints.md](file:///i:/Nutstore/1/blog/point-cloud/20-iss-keypoints.md)（ISS 固有形状特征）。

三维点云处理（十九）：Harris 角点检测 (2D/3D/6D) ​

核心算法特点对比 ​

一、特征点的三大准则与几何分类 ​

1.1 什么是好的特征点？ ​

1.2 局部几何形态与特征谱分类 ​

二、Harris 2D 角点检测原理（回顾） ​

三、Harris 3D：从 2D 到 3D 点云的推广 ​

3.1 路径一：曲面拟合与梯度张量法（Sipiran 路径） ​

3.2 路径二：3D 协方差矩阵特征谱法 ​

3.3 Harris 3D 参数选择的工程经验 ​

四、Harris 6D 与强度扩展（Harris with Intensity / Normals） ​

4.1 几何角点的局限性 ​

4.2 Harris 6D 联合协方差 ​

五、纯 Python 与 Open3D 代码实现 ​

六、总结 ​