17.2 主成分分析 (Principal Component Analysis) §

相关章节 §

• 17.1 最优包围盒
• 17.3 二次规划求解器
• 18.1 点云处理与分析

17.2.1 理论基础 §

数学基础 §

主成分分析(PCA)是一种统计方法，通过正交变换将一组可能相关的变量转换为一组线性不相关的变量（主成分）。在几何处理中，PCA用于：

1. 寻找主要变化方向：确定点集分布的主轴
2. 降维：将高维数据投影到低维子空间
3. 对齐：将物体对齐到标准坐标系

协方差矩阵 §

给定 $n$ 个 $d$ 维点 $\{p_1,p_2,...,p_n\}$，首先计算质心：

$c=\frac{1}{n}{\sum }_{i=1}^n{p}_i$

协方差矩阵 $C$ 定义为：

$C=\frac{1}{n}{\sum }_{i=1}^n(p_i-c)(p_i-c{)}^T$

对于三维点，协方差矩阵是对称的：

c_{xx} & c_{xy} & c_{xz} \\ c_{xy} & c_{yy} & c_{yz} \\ c_{xz} & c_{yz} & c_{zz} \end{bmatrix}$$ 其中： - $c_{xx} = \sum (x_i - c_x)^2$ - $c_{xy} = \sum (x_i - c_x)(y_i - c_y)$ - 以此类推 ### 特征分解 PCA通过对协方差矩阵进行特征分解来获得主成分： $$C = V \Lambda V^T$$ 其中： - $\Lambda = \text{diag}(\lambda_1, \lambda_2, \lambda_3)$ 是特征值对角矩阵（$\lambda_1 \geq \lambda_2 \geq \lambda_3$） - $V = [v_1, v_2, v_3]$ 是特征向量矩阵 特征向量 $v_1$ 对应最大特征值 $\lambda_1$，表示数据变化最大的方向（第一主成分）。 ### 拟合质量 CGAL的PCA实现返回一个拟合质量值 $q \in [0, 1]$： - **直线拟合**：$q = 1 - \frac{\lambda_2}{\lambda_3}$ - **平面拟合**：$q = 1 - \frac{\lambda_1}{\lambda_2}$ $q = 1$ 表示完美拟合（所有点共线或共面），$q = 0$ 表示各向同性分布。 ### 加权PCA CGAL支持基于几何原语的加权PCA： | 几何原语 | 维度 | 权重 | 二阶矩 | |----------|------|------|--------| | 点 | 0 | 1 | - | | 线段 | 1 | 长度 | $[1/3, 1/6; 1/6, 1/3]$ | | 三角形 | 2 | 面积 | $[1/12, 1/24; 1/24, 1/12]$ | | 四面体 | 3 | 体积 | $[1/60, 1/120; 1/120, 1/60]$ | | 球 | 3 | 体积 | $4/15 \cdot I$ | | 长方体 | 3 | 体积 | $[1/3, 1/4; 1/4, 1/3]$ | ## 17.2.2 架构分析 ### 模块结构 ``` Principal_component_analysis/ ├── include/CGAL/ │ ├── linear_least_squares_fitting_2.h # 2D主入口 │ ├── linear_least_squares_fitting_3.h # 3D主入口 │ ├── linear_least_squares_fitting_points_2.h │ ├── linear_least_squares_fitting_points_3.h │ ├── linear_least_squares_fitting_segments_2.h │ ├── linear_least_squares_fitting_segments_3.h │ ├── linear_least_squares_fitting_triangles_2.h │ ├── linear_least_squares_fitting_triangles_3.h │ ├── linear_least_squares_fitting_tetrahedra_3.h │ ├── linear_least_squares_fitting_cuboids_3.h │ ├── linear_least_squares_fitting_spheres_3.h │ ├── linear_least_squares_fitting_circles_2.h │ ├── linear_least_squares_fitting_rectangles_2.h │ ├── PCA_util.h # 核心工具函数 │ └── PCA_util_Eigen.h # Eigen对角化 └── examples/ ├── linear_least_squares_fitting_points_2.cpp └── linear_least_squares_fitting_triangles_3.cpp ``` ### 核心接口 **2D版本**： ```cpp template < typename InputIterator, typename Line, typename Tag > inline typename Kernel_traits<Line>::Kernel::FT linear_least_squares_fitting_2( InputIterator first, // 输入范围起点 InputIterator beyond, // 输入范围终点 Line& line, // 输出：拟合直线 const Tag& tag // 维度标签 ); ``` **3D版本**： ```cpp template < typename InputIterator, typename Object, typename Tag > inline typename Kernel_traits<Object>::Kernel::FT linear_least_squares_fitting_3( InputIterator first, InputIterator beyond, Object& object, // Plane_3 或 Line_3 typename Kernel::Point_3& centroid, // 输出质心 const Tag& tag ); ``` ### 维度标签 ```cpp CGAL::Dimension_tag<0>() // 点集 CGAL::Dimension_tag<1>() // 线段集 CGAL::Dimension_tag<2>() // 三角形/矩形/圆 CGAL::Dimension_tag<3>() // 四面体/长方体/球 ``` ### 对角化Traits ```cpp // 默认使用内部实现或Eigen typename Default_diagonalize_traits<typename Kernel::FT, 3> // 自定义对角化 template<typename FT, unsigned int dim> struct MyDiagonalizeTraits { typedef std::array<FT, dim> Vector; typedef std::array<FT, dim*(dim+1)/2> Covariance_matrix; static void diagonalize_selfadjoint_covariance_matrix( const Covariance_matrix& cov, Vector& eigenvalues, Vector& eigenvectors ); }; ``` ## 17.2.3 实现细节 ### 核心算法流程 ```cpp // 来自 linear_least_squares_fitting_points_3.h template < typename InputIterator, typename K, typename DiagonalizeTraits > typename K::FT linear_least_squares_fitting_3( InputIterator first, InputIterator beyond, typename K::Plane_3& plane, typename K::Point_3& c, const typename K::Point_3*, const K& k, const CGAL::Dimension_tag<0>& tag, const DiagonalizeTraits& diagonalize_traits) { typedef typename K::Point_3 Point; // 1. 计算质心 c = centroid(first, beyond, K(), tag); // 2. 组装协方差矩阵 typename DiagonalizeTraits::Covariance_matrix covariance = {{ 0., 0., 0., 0., 0., 0. }}; assemble_covariance_matrix_3(first, beyond, covariance, c, k, (Point*) nullptr, tag, diagonalize_traits); // 3. 计算拟合平面 return fitting_plane_3(covariance, c, plane, k, diagonalize_traits); } ``` ### 协方差矩阵组装 ```cpp // 来自 PCA_util.h // 点集的协方差矩阵 template < typename InputIterator, typename K, typename DiagonalizeTraits > void assemble_covariance_matrix_3( InputIterator first, InputIterator beyond, typename DiagonalizeTraits::Covariance_matrix& covariance, const typename K::Point_3& c, const K& k, const typename K::Point_3*, const CGAL::Dimension_tag<0>&, const DiagonalizeTraits&) { typedef typename K::FT FT; typedef typename K::Point_3 Point; typedef typename K::Vector_3 Vector; // 矩阵编号（上三角存储）： // 0 1 2 // 3 4 // 5 covariance[0] = covariance[1] = covariance[2] = covariance[3] = covariance[4] = covariance[5] = (FT)0.0; for(InputIterator it = first; it \!= beyond; ++it) { const Point& p = *it; Vector d = k.construct_vector_3_object()(c, p); covariance[0] += d.x() * d.x(); // xx covariance[1] += d.x() * d.y(); // xy covariance[2] += d.x() * d.z(); // xz covariance[3] += d.y() * d.y(); // yy covariance[4] += d.y() * d.z(); // yz covariance[5] += d.z() * d.z(); // zz } } ``` ### 三角形集的协方差矩阵 ```cpp // 来自 PCA_util.h template < typename InputIterator, typename K, typename DiagonalizeTraits > void assemble_covariance_matrix_3( InputIterator first, InputIterator beyond, typename DiagonalizeTraits::Covariance_matrix& covariance, const typename K::Point_3& c, const K&, const typename K::Triangle_3*, const CGAL::Dimension_tag<2>&, const DiagonalizeTraits&) { typedef typename DiagonalizeTraits::Vector::value_type FT; typedef typename K::Triangle_3 Triangle; typedef typename CGAL::Linear_algebraCd<FT> LA; typedef typename LA::Matrix Matrix; FT mass = FT(0.0); // 标准二阶矩（关于原点） FT temp[9] = {FT(1.0/12.0), FT(1.0/24.0), FT(1.0/24.0), FT(1.0/24.0), FT(1.0/12.0), FT(1.0/24.0), FT(1.0/24.0), FT(1.0/24.0), FT(1.0/12.0)}; Matrix moment = init_matrix<FT>(3, temp); for(InputIterator it = first; it \!= beyond; ++it) { const Triangle& t = *it; // 构造变换矩阵 FT delta[9] = {FT(t[0].x()), FT(t[1].x()), FT(t[2].x()), FT(t[0].y()), FT(t[1].y()), FT(t[2].y()), FT(t[0].z()), FT(t[1].z()), FT(t[2].z())}; Matrix transformation = init_matrix<FT>(3, delta); FT area = FT(std::sqrt(t.squared_area())); if(area == (FT)0.0) continue; // 仿射变换：2 * area * T * moment * T^T transformation = 2 * area * transformation * moment * LA::transpose(transformation); // 累加到协方差矩阵 covariance[0] += transformation[0][0]; covariance[1] += transformation[1][0]; covariance[2] += transformation[2][0]; covariance[3] += transformation[1][1]; covariance[4] += transformation[2][1]; covariance[5] += transformation[2][2]; mass += area; } // 将原点矩平移到质心 covariance[0] += mass * FT(-1.0 * c.x() * c.x()); covariance[1] += mass * FT(-1.0 * c.x() * c.y()); // ... 其他项 } ``` ### 平面拟合 ```cpp // 来自 PCA_util.h template < typename K, typename DiagonalizeTraits > typename K::FT fitting_plane_3( typename DiagonalizeTraits::Covariance_matrix& covariance, const typename K::Point_3& c, typename K::Plane_3& plane, const K&, const DiagonalizeTraits&) { typedef typename K::FT FT; typedef typename K::Plane_3 Plane; typedef typename K::Vector_3 Vector; // 求解特征值和特征向量 typename DiagonalizeTraits::Vector eigenvalues = {{ 0., 0., 0. }}; typename DiagonalizeTraits::Matrix eigenvectors = {{ 0., 0., 0., 0., 0., 0., 0., 0., 0. }}; DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix( covariance, eigenvalues, eigenvectors); // 退化情况：各向同性 if(eigenvalues[0] == eigenvalues[1] && eigenvalues[1] == eigenvalues[2]) { plane = Plane(c, Vector(FT(0), FT(0), FT(1))); return FT(0); } else { // 最小特征值对应的特征向量是法向 Vector normal(eigenvectors[0], eigenvectors[1], eigenvectors[2]); plane = Plane(c, normal); if (eigenvalues[1] == 0) return FT(0); // 共线情况 else return FT(1) - eigenvalues[0] / eigenvalues[1]; // 拟合质量 } } ``` ### 直线拟合 ```cpp // 来自 PCA_util.h template < typename K, typename DiagonalizeTraits > typename K::FT fitting_line_3( typename DiagonalizeTraits::Covariance_matrix& covariance, const typename K::Point_3& c, typename K::Line_3& line, const K&, const DiagonalizeTraits&) { typedef typename K::FT FT; typedef typename K::Line_3 Line; typedef typename K::Vector_3 Vector; typename DiagonalizeTraits::Vector eigenvalues = {{ 0., 0., 0. }}; typename DiagonalizeTraits::Matrix eigenvectors = {{ 0., 0., 0., 0., 0., 0., 0., 0., 0. }}; DiagonalizeTraits::diagonalize_selfadjoint_covariance_matrix( covariance, eigenvalues, eigenvectors); // 各向同性：返回默认x轴方向 if(eigenvalues[0] == eigenvalues[1] && eigenvalues[0] == eigenvalues[2]) { line = Line(c, Vector(FT(1), FT(0), FT(0))); return (FT)0.0; } else { // 最大特征值对应的特征向量是方向 Vector direction(eigenvectors[6], eigenvectors[7], eigenvectors[8]); line = Line(c, direction); return (FT)1.0 - eigenvalues[1] / eigenvalues[2]; // 拟合质量 } } ``` ## 17.2.4 使用示例 ### 2D点集直线拟合 ```cpp #include <CGAL/Simple_cartesian.h> #include <CGAL/linear_least_squares_fitting_2.h> #include <vector> typedef double FT; typedef CGAL::Simple_cartesian<FT> K; typedef K::Line_2 Line; typedef K::Point_2 Point; int main() { std::vector<Point> points; points.push_back(Point(1.0, 2.0)); points.push_back(Point(3.0, 4.0)); points.push_back(Point(5.0, 6.0)); Line line; // 拟合直线，Dimension_tag<0>表示点集 FT quality = linear_least_squares_fitting_2( points.begin(), points.end(), line, CGAL::Dimension_tag<0>() ); std::cout << "Fitted line: " << line << std::endl; std::cout << "Quality: " << quality << std::endl; return 0; } ``` ### 3D点集平面拟合 ```cpp #include <CGAL/Simple_cartesian.h> #include <CGAL/linear_least_squares_fitting_3.h> #include <vector> typedef double FT; typedef CGAL::Simple_cartesian<FT> K; typedef K::Plane_3 Plane; typedef K::Point_3 Point; int main() { std::vector<Point> points; // 添加近似共面的点 for(int i = 0; i < 10; ++i) { for(int j = 0; j < 10; ++j) { points.push_back(Point(i, j, 0.1 * (i + j))); // 近似z=0平面 } } Plane plane; Point centroid; FT quality = CGAL::linear_least_squares_fitting_3( points.begin(), points.end(), plane, centroid, CGAL::Dimension_tag<0>() ); std::cout << "Plane: " << plane << std::endl; std::cout << "Centroid: " << centroid << std::endl; std::cout << "Quality: " << quality << std::endl; return 0; } ``` ### 三角形网格PCA ```cpp #include <CGAL/Simple_cartesian.h> #include <CGAL/linear_least_squares_fitting_3.h> #include <CGAL/Triangle_3.h> #include <vector> typedef double FT; typedef CGAL::Simple_cartesian<FT> K; typedef K::Plane_3 Plane; typedef K::Point_3 Point; typedef K::Triangle_3 Triangle; int main() { std::vector<Triangle> triangles; // 添加三角形（例如从网格中提取） triangles.push_back(Triangle( Point(0, 0, 0), Point(1, 0, 0), Point(0.5, 1, 0) )); triangles.push_back(Triangle( Point(1, 0, 0), Point(1, 1, 0), Point(0.5, 1, 0) )); // ... 更多三角形 Plane plane; Point centroid; // Dimension_tag<2>表示三角形集（2维对象） FT quality = CGAL::linear_least_squares_fitting_3( triangles.begin(), triangles.end(), plane, centroid, CGAL::Dimension_tag<2>() ); std::cout << "Mesh PCA plane: " << plane << std::endl; std::cout << "Quality: " << quality << std::endl; return 0; } ``` ### 使用Eigen进行对角化 ```cpp #include <CGAL/Simple_cartesian.h> #include <CGAL/linear_least_squares_fitting_3.h> #include <CGAL/PCA_util_Eigen.h> // Eigen对角化支持 // 使用Eigen进行更高效的数值计算 typedef CGAL::Simple_cartesian<double> K; typedef K::Plane_3 Plane; typedef K::Point_3 Point; int main() { std::vector<Point> points; // ... 填充点 Plane plane; Point centroid; // 显式指定使用Eigen对角化 FT quality = CGAL::linear_least_squares_fitting_3( points.begin(), points.end(), plane, centroid, CGAL::Dimension_tag<0>(), K(), CGAL::Eigen_diagonalize_traits<double, 3>() ); return 0; } ``` ## 17.2.5 复杂度分析 ### 时间复杂度 | 操作 | 复杂度 | 说明 | |------|--------|------| | 质心计算 | $O(n)$ | 单次遍历 | | 协方差矩阵组装 | $O(n)$ | 点集；$O(m)$ 对于 $m$ 个三角形 | | 特征分解 | $O(d^3)$ | $d$ 为维度（2或3），常数时间 | | **总体** | $O(n)$ | 线性于输入规模 | ### 空间复杂度 - 输入存储：$O(n)$ - 协方差矩阵：$O(1)$（固定 $3 \times 3$） - 特征向量/值：$O(1)$ - **总体**：$O(n)$（主要由输入决定） ### 数值稳定性 1. **协方差计算**：使用双精度累加，避免单精度溢出 2. **特征分解**： - 默认使用CGAL内部实现（基于Jacobi迭代） - 可选Eigen后端（更稳定的QR算法） 3. **退化处理**： - 检测各向同性分布 - 检测共线/共面情况 - 返回默认方向而非失败 ### 精度考虑 - **浮点类型**：支持 `float`、`double`、精确数类型 - **大规模数据**：对于 $n > 10^6$，考虑使用Kahan求和 - **病态数据**：当点分布接近球对称时，特征值接近，方向可能不稳定 ## 17.2.6 应用场景 ### 1. 物体对齐 ```cpp // 将物体对齐到主坐标系 Plane plane; Point centroid; CGAL::linear_least_squares_fitting_3(points.begin(), points.end(), plane, centroid, CGAL::Dimension_tag<0>()); // 构造变换矩阵 Vector n = plane.orthogonal_vector(); Vector u = plane.base1(); // 平面内第一个方向 Vector v = plane.base2(); // 平面内第二个方向 // 将物体平移到原点并旋转对齐 // ... ``` ### 2. 法向估计 ```cpp // 从点云估计法向 for each point p in point_cloud { // 获取k近邻 std::vector<Point> neighbors = kdtree.nearest_neighbors(p, k); Plane local_plane; CGAL::linear_least_squares_fitting_3(neighbors.begin(), neighbors.end(), local_plane, CGAL::Dimension_tag<0>()); Vector normal = local_plane.orthogonal_vector(); // 存储法向... } ``` ### 3. 对称性检测 ```cpp // 通过PCA分析物体对称性 Line axis; CGAL::linear_least_squares_fitting_3(points.begin(), points.end(), axis, CGAL::Dimension_tag<0>()); // 分析沿主轴的分布... ``` ### 4. 与OBB的关系 PCA可用于计算快速但非最优的OBB： ```cpp // PCA-based OBB（快速但非最优） Plane plane; CGAL::linear_least_squares_fitting_3(points.begin(), points.end(), plane, CGAL::Dimension_tag<0>()); // 使用PCA主轴作为OBB方向 // 比CGAL::oriented_bounding_box快但包围盒体积更大 ``` ## 17.2.7 扩展与变体 ### 鲁棒PCA 对于含噪声或离群点的数据，可考虑： - RANSAC-based拟合 - 加权PCA（根据置信度加权） - 迭代重加权最小二乘 ### 核PCA 对于非线性结构，CGAL当前未直接支持，但可通过预处理实现： - 将数据映射到高维特征空间 - 在特征空间执行标准PCA