三维点云法线估计与几何特征计算

点云法向量（Point Cloud Normal）记录了三维表面局部切平面的正交朝向，是三维几何处理、光照渲染、点云配准（如 Point-to-Plane ICP）以及曲面重建（如 Poisson 泊松重建）中的核心几何特征。本文系统讲解基于主成分分析（PCA）的点云法线估计原理、局部曲率计算、法线朝向一致化传播算法以及通用 C++ 代码实现。

---

1. 点云法线与曲率的几何意义

无网格拓扑信息的原始三维点云仅包含独立的空间坐标 。在进行点云配准与几何分析时，单一的坐标点无法提供表面的弯曲与朝向趋势。

通过在每个点的局部邻域内拟合一个微观切平面（Tangential Plane），切平面的法向量 即可表征该点处的表面朝向；而邻域点沿法向的偏离程度则反映了该处的曲率（Curvature）。

---

2. 基于 PCA 的局部切平面拟合

2.1 质心中心化与协方差矩阵构造

对点云中任意一点 ，利用 KD-Tree 提取其半径 邻域内的 个最近邻点集 ：

1. 计算邻域几何质心：

2. 构建 样本协方差矩阵 (Covariance Matrix)：

   协方差矩阵 是一个对称正定矩阵，描述了邻域点集在空间三个主轴方向上的方差分布。

2.2 特征值分解与几何解释

对矩阵 进行特征值分解，求得三个升序排列的特征值与对应的单位特征向量：

• 特征向量 ：对应最大与次特征值 ，代表邻域点分布方差最大的两个正交方向，张成了局部的拟合切平面。

• 特征向量 ：对应最小特征值 ，代表偏离切平面方差最小的方向。因此，最小特征向量 即为该点处的单位法向量 ：

2.3 局部曲率 (Curvature) 估计

局部表面曲率变异度（Surface Variation）可以通过最小特征值占总方差的比例来定量估计：

• 当 时：说明 ，邻域点完全分布在严格的平面上。
• 当 较大时：说明表面起伏剧烈、存在锐利边缘或高频噪点。

---

3. 法线朝向一致化 (Orientation Consistency)

3.1 符号二义性问题

由于特征值分解中 与 等价，计算出的无向法向量 存在 的 符号二义性。如果相邻点的法向量朝向混乱（一部分指向模型外，一部分指向模型内），依赖法向投影的算法（如 Point-to-Plane ICP、Poisson 重建）将彻底失效。

3.2 视点对齐法 (Viewpoint Alignment)

若已知相机的采集视点坐标 （例如视点位于传感器原点 ），可通过向量内积判定并翻转法向，使其始终朝向视点：

若点乘结果小于 0，则执行法向翻转：。

3.3 基于最小生成树 (MST) 的拓扑传播

对于无视点信息的闭合体点云，需使用图论中的**最小生成树（Minimum Spanning Tree, MST）**算法进行法向平滑传播：

1. 构建以点云为顶点的无向图，边权值设定为相邻点法向量正交程度的倒数：

2. 求解该图的最小生成树。从边界或最外侧起始点开始，沿着生成树边缘进行广度优先遍历（BFS）： 若相邻节点满足 ，则将 翻转为 。

---

4. 基于 Eigen 的 C++ 通用代码实现

#include <iostream>
#include <vector>
#include <Eigen/Dense>

struct Point3D {
    Eigen::Vector3d position;
    Eigen::Vector3d normal{0.0, 0.0, 0.0};
    double curvature{0.0};
};

// 单点 PCA 法线估计函数
void compute_point_normal_pca(Point3D& target_point, const std::vector<Eigen::Vector3d>& neighbors) {
    if (neighbors.size() < 3) {
        return;
    }

    // 1. 计算质心
    Eigen::Vector3d centroid = Eigen::Vector3d::Zero();
    for (const auto& pt : neighbors) {
        centroid += pt;
    }
    centroid /= static_cast<double>(neighbors.size());

    // 2. 构造 $3 \times 3$ 协方差矩阵
    Eigen::Matrix3d covariance = Eigen::Matrix3d::Zero();
    for (const auto& pt : neighbors) {
        Eigen::Vector3d de_centered = pt - centroid;
        covariance += de_centered * de_centered.transpose();
    }
    covariance /= static_cast<double>(neighbors.size());

    // 3. 特征值分解
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> solver(covariance);
    if (solver.info() != Eigen::Success) {
        return;
    }

    // 特征值按升序排列
    Eigen::Vector3d eigenvalues = solver.eigenvalues();
    Eigen::Matrix3d eigenvectors = solver.eigenvectors();

    // 4. 提取法向量与曲率
    target_point.normal = eigenvectors.col(0).normalized(); // 最小特征值对应的特征向量
    
    double sum_evals = eigenvalues.sum();
    if (sum_evals > 1e-12) {
        target_point.curvature = eigenvalues(0) / sum_evals;
    }
}

int main() {
    Point3D target;
    target.position = Eigen::Vector3d(0.0, 0.0, 0.0);

    // 构造模拟平面邻域点集
    std::vector<Eigen::Vector3d> neighbors = {
        {-1.0, -1.0,  0.02},
        { 1.0, -1.0, -0.01},
        {-1.0,  1.0,  0.01},
        { 1.0,  1.0, -0.02},
        { 0.0,  0.0,  0.00}
    };

    compute_point_normal_pca(target, neighbors);

    std::cout << "Computed Normal: (" 
              << target.normal.x() << ", " 
              << target.normal.y() << ", " 
              << target.normal.z() << ")\n";
    std::cout << "Estimated Curvature: " << target.curvature << "\n";

    return 0;
}

---

5. 核心原理与高频追问

FAQ 1：邻域搜索半径 的选取对法线计算有什么影响？如何自适应确定？

答：

• 半径 过小：邻域内包含的点数量不足（小于 3 个），协方差矩阵欠秩；且对扫描设备的随机噪声极度敏感，导致计算出的法向量剧烈抖动。
• 半径 过大：邻域会跨越模型的锐利边缘（Sharp Edges）或薄壁结构，导致拟合出的平面“模糊化”，平滑掉微小的几何细节。
• 自适应确定方法：通常基于八叉树或点云包围盒（Bounding Box）对角线长度 与点云总数 启发式估计基准半径：；或使用 近邻搜索（如取 ）确保各点局部邻域密度均衡。

FAQ 2：为什么拟合切平面的代价目标恰好等价于协方差矩阵的特征值分解？

答： 设切平面单位法向量为 ，平面方程为 。邻域点 到切平面的垂直欧氏距离为 。 最小二乘切平面拟合目标为最小化垂直距离平方和：

根据瑞利商（Rayleigh Quotient）定理，在 的约束下，二次型 取最小值时的向量 恰好是矩阵 最小特征值对应的特征向量。