PCA（主成分分析）获取BoundingBox代码分析

引言

最近要用到PCA获取目标点云的BoundingBox，但是网上给出的有关PCA的代码大都太简洁了，我觉得可能是大佬觉得比较简单，没有详细描述。这里记录一下自己的探究结果，方便大家理解。欢迎留言讨论，如有错误，请批评指正。

理论

PCA(Principal Component Analysis)是一种十分常用的数据降维方法，其主要思想是通过线性变换的方法，将一组高维数据投影成互不相关的低维数据，保留原始数据最关键、最具代表性的数据，以压缩数据量、加速运算过程。

通过问题概述PCA推导过程：

PCA的目标是对原始数据进行降维，怎么降维呢？
由于数据是基于一组N维正交基表示的，因此可以将原始数据进行重新投影，投影到新的一组K维正交基上，当K<N时，即可达到降维的目的。
那么新的正交基需要满足哪些条件呢？
为了最大程度保留原始数据的信息，需要使降维以后的数据尽可能的分散。也就是说，新的正交基需要保证降维数据的分散程度达到最大。而数据的分散程度通常使用方差来描述，即新的正交基要保证降维后数据的方差和达到最大。当然，还需要保证基之间相互正交。
具体怎么求呢？
在实际求目标正交基的时候，需要同时保证降维后数据的方差和最大、基之间相互正交。在数学上，数据的协方差矩阵的对角线元素描述的是数据的方差，非对角线元素描述的是数据的协方差，即变量之间的相互关系。因此上述两条件可以统一到协方差矩阵里，而且降维前后协方差矩阵满足关系：D=PCPTD=PCP^TD=PCPT，DDD为降维后数据协方差矩阵，CCC为原数据协方差矩阵。因此，目标就变成了求解矩阵PPP，使得上述等式中DDD为对角阵（对CCC进行对角化），且PPP为PCPTPCP^TPCPT中最大的前KKK个（保证方差和最大）对角线元素对应的特征向量组成的矩阵。
算法步骤
1）求原数据协方差矩阵C=1/m∗XXTC=1/m*XX^TC=1/m∗XXT；
2）求出协方差矩阵CCC的特征值及对应的特征向量；
3）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前KKK行组成矩阵；
4）Y=PXY=PXY=PX即为降维到KKK维后的数据.

PCA详细理论介绍，可参考如下博客，讲解的十分清楚：

PCA的数学原理
主成分分析（PCA）原理总结
线性代数之主成分分析(PCA)算法
主成分分析(PCA)原理及推导–参考其拉格朗日求解法

代码

环境：Win10+VS2015+PCL1.8

代码流程：

导入点云–>计算质心和协方差矩阵–>计算主方向(即特征向量)；
将点云转到参考坐标系，保证点云主方向与参考坐标系坐标轴重合；
计算转换后点云的边界值(min_pt、max_pt)–>计算BoundingBox形心(mean_diag)–>将mean_diag转换到原点云；
最后根据变换关系及形心位置addCube()，在原点云形心处添加BoundingBox。

其实，在流程中的第一步就已经完成了点云主方向的计算，后续过程只是为了保证BoundingBox能完整的包围点云。

#include <vtkAutoInit.h>
VTK_MODULE_INIT(vtkRenderingOpenGL);
VTK_MODULE_INIT(vtkInteractionStyle);
VTK_MODULE_INIT(vtkRenderingFreeType);#include <iostream>
#include <string>
#include <pcl/io/pcd_io.h>
#include <pcl/io/ply_io.h>
#include <pcl/point_cloud.h>
#include <pcl/point_types.h>
#include <Eigen/Core>
#include <pcl/common/transforms.h>
#include <pcl/common/common.h>
#include <pcl/visualization/pcl_visualizer.h>using namespace std;
typedef pcl::PointXYZ PointType;int main(int argc, char **argv)
{// 导入点云pcl::PointCloud<PointType>::Ptr original_cloud(new pcl::PointCloud<PointType>());std::string fileName("lamppost.pcd");pcl::io::loadPCDFile(fileName, *original_cloud);// PCA：计算主方向Eigen::Vector4f centroid;                            // 质心pcl::compute3DCentroid(*original_cloud, centroid); // 齐次坐标，（c0,c1,c2,1）Eigen::Matrix3f covariance;    computeCovarianceMatrixNormalized(*original_cloud, centroid, covariance);       // 计算归一化协方差矩阵// 计算主方向：特征向量和特征值Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigen_solver(covariance, Eigen::ComputeEigenvectors);Eigen::Matrix3f eigen_vectors = eigen_solver.eigenvectors();//Eigen::Vector3f eigen_values = eigen_solver.eigenvalues();eigen_vectors.col(2) = eigen_vectors.col(0).cross(eigen_vectors.col(1));  // 校正主方向间垂直（特征向量方向： (e0, e1, e0 × e1) --- note: e0 × e1 = +/- e2）// 转到参考坐标系，将点云主方向与参考坐标系的坐标轴进行对齐Eigen::Matrix4f transformation(Eigen::Matrix4f::Identity());transformation.block<3, 3>(0, 0) = eigen_vectors.transpose();                                      // R^(-1) = R^Ttransformation.block<3, 1>(0, 3) = -1.f * (transformation.block<3, 3>(0, 0) * centroid.head<3>());   // t^(-1) = -R^T * tpcl::PointCloud<PointType> transformed_cloud;    // 变换后的点云pcl::transformPointCloud(*original_cloud, transformed_cloud, transformation);PointType min_pt, max_pt;                     // 沿参考坐标系坐标轴的边界值pcl::getMinMax3D(transformed_cloud, min_pt, max_pt);const Eigen::Vector3f mean_diag = 0.5f*(max_pt.getVector3fMap() + min_pt.getVector3fMap());   // 形心// 参考坐标系到主方向坐标系的变换关系const Eigen::Quaternionf qfinal(eigen_vectors);const Eigen::Vector3f tfinal = eigen_vectors * mean_diag + centroid.head<3>();// 显示结果pcl::visualization::PCLVisualizer viewer;viewer.addPointCloud(original_cloud);viewer.addCoordinateSystem();// 显示点云主方向Eigen::Vector3f whd;      // 3个方向尺寸：宽高深whd = max_pt.getVector3fMap() - min_pt.getVector3fMap();// getVector3fMap:返回Eigen::Map<Eigen::Vector3f> float scale = (whd(0) + whd(1) + whd(2)) / 3;         // 点云平均尺度，用于设置主方向箭头大小// std::cout << "width/heigth/depth：" << whd << endl;PointType cp;           // 箭头由质心分别指向pirncipal_dir_X、pirncipal_dir_Y、pirncipal_dir_Zcp.x = centroid(0);cp.y = centroid(1);cp.z = centroid(2);PointType principal_dir_X;principal_dir_X.x = scale * eigen_vectors(0, 0) + cp.x;principal_dir_X.y = scale * eigen_vectors(1, 0) + cp.y;principal_dir_X.z = scale * eigen_vectors(2, 0) + cp.z;PointType principal_dir_Y;principal_dir_Y.x = scale * eigen_vectors(0, 1) + cp.x;principal_dir_Y.y = scale * eigen_vectors(1, 1) + cp.y;principal_dir_Y.z = scale * eigen_vectors(2, 1) + cp.z;PointType principal_dir_Z;principal_dir_Z.x = scale * eigen_vectors(0, 2) + cp.x;principal_dir_Z.y = scale * eigen_vectors(1, 2) + cp.y;principal_dir_Z.z = scale * eigen_vectors(2, 2) + cp.z;viewer.addArrow(principal_dir_X, cp, 1.0, 0.0, 0.0, false, "arrow_x");      // 箭头附在起点上viewer.addArrow(principal_dir_Y, cp, 0.0, 1.0, 0.0, false, "arrow_y");viewer.addArrow(principal_dir_Z, cp, 0.0, 0.0, 1.0, false, "arrow_z");// 显示包围盒，并设置包围盒属性，以显示透明度viewer.addCube(tfinal, qfinal, whd(0), whd(1), whd(2), "bbox");viewer.setShapeRenderingProperties(pcl::visualization::PCL_VISUALIZER_REPRESENTATION, pcl::visualization::PCL_VISUALIZER_REPRESENTATION_WIREFRAME, "bbox");viewer.setShapeRenderingProperties(pcl::visualization::PCL_VISUALIZER_COLOR, 1.0, 0.0, 0.0, "bbox");while (!viewer.wasStopped()){viewer.spinOnce(100);}return 0;
}

另外，还可以使用pcl中的pca接口计算点云的主方向：

// 使用pcl中的pca接口
pcl::PointCloud<pcl::PointXYZ>::Ptr cloudPCAprojection (new pcl::PointCloud<pcl::PointXYZ>);
pcl::PCA<pcl::PointXYZ> pca;
pca.setInputCloud(cloudSegmented);
pca.project(*cloudSegmented, *cloudPCAprojection);
std::cerr << std::endl << "EigenVectors: " << pca.getEigenVectors() << std::endl; // 计算特征向量
std::cerr << std::endl << "EigenValues: " << pca.getEigenValues() << std::endl;       // 计算特征值

结果

点云下载地址：lamppost.pcd

分析

为什么要将点云变换到原点，而不是直接对原点云进行操作？

把结论写在前面：如果只求原点云的主方向，则直接对原点云进行操作就可以了。之所以要将点云转换到原点，是为了后续能够正确的显示包围盒。原因在于getMinMax3D()计算得到的min_pt与max_pt代表的是原点坐标系下的边界值，而不是点云主方向坐标系下的值。

首先需要明确，原点云的特征向量组成了其主方向坐标系(称为EEE系)，该坐标系原点即是点云质心。同时，特征向量和质心坐标组成的3∗43*43∗4的变换矩阵，即代表由参考坐标系(称为WWW系)到主方向坐标系的变换矩阵TEWT_{EW}TEW。

将点云从EEE系变换到WWW系，此时的点云的3个主方向已经与WWW系的坐标轴重合，同时其质心也与WWW系原点重合；在这种状态下，求点云3个方向上的边界值才是准确的。

因为pcl::getMinMax3D(*cloud, min_pt, max_pt)计算的是点云沿WWW系3个坐标轴方向上的最小、最大值，所以如果点云本身的EEE系没有与WWW系重合，即未将点云从EEE系变换到WWW系，则通过getMinMax3D()计算得到的min_pt与max_pt代表的就不是点云主方向上的边界值。

因此，如果直接对原点云进行操作，而不变换到WWW系进行操作，则得到的BoundingBox的位置和尺寸都会有偏差，其原因就在于pcl::getMinMax3D(*point_cloud_ptr, min_pt, max_pt)得到的min_pt、max_pt结果不对，对应的代码如下：

PointType min_pt, max_pt;
pcl::getMinMax3D(*point_cloud_ptr, min_pt, max_pt); // 对原点云进行操作
const Eigen::Vector3f mean_diag = 0.5f*(max_pt.getVector3fMap() + min_pt.getVector3fMap());// T_EW
const Eigen::Quaternionf qfinal(eigen_vectors);
const Eigen::Vector3f tfinal = mean_diag ;

参考

PCA最小包围盒：Find minimum oriented bounding box of point cloud (C++ AND PCL)
网上很多相似代码的起源：Nicola Fioraio
AABB&OBB包围盒：Moment of inertia and eccentricity based descriptors
C++_Eigen函数库用法笔记——Block Operations

PCA（主成分分析）获取BoundingBox代码分析相关推荐

中药材鉴别-方法：聚类；PCA 主成分分析；线性判别式分析；判别式检验
基于线性判别式的中药材鉴别问题的数学模型摘要本文旨在讨论如何利用中药材的光谱特征鉴别药材的种类及产地,主要运用系统聚类,PCA 主成分分析,线性判别,判别式运用等方法,使用了 MATLAB,Ex ...
PCA主成分分析教程（origin分析绘制，无须R语言）
PCA主成分分析教程(origin分析&绘制,无须R语言) 相关性分析,相关的介绍内容大家自行搜索资料即可,这里不给大家过多阐述. 案例解读 PCA作为常见的一种聚类分析方法,在很多SCI论文 ...
这次终于理解了PCA主成分分析（附代码）
在降维过程中,我们会减少特征的数量,这意味着删除数据,数据量变少则表示模型可以获取的信息会变少,模型的表现可能会因此受影响.同时,在高维数据中,必然有一些特征是不带有有效的信息的(比如噪音),或者有一 ...
20189200余超 2018-2019-2 移动平台应用开发实践作项目代码分析
20189200余超 2018-2019-2 移动平台应用开发实践作项目代码分析项目名称小说阅读器项目功能注册登录用户信息.用户密码.用户图像修改书籍分类书架书籍搜索(作者名或书籍名) ...
数据分析实战：python热门音乐分析附代码+数据 +论文（PCA 主成分分析，sklearn 机器学习，pytorch 神经网络，k-means 聚类，Librosa 音频处理，midi 音序)
项目概述: 本选取了抖音当下最热门的 400 首音乐,通过一系列方法提取每首歌的波形特征,再经过降维以及机器学习等手段,进行无监督学习对音乐数据进行聚类的同时训练并使用监督学习分类器进行音乐流派分类, ...
复现nature communication PCA原图|代码分析（一）
复现PCA原图之bulk RNA-seq NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程) ...
PCA主成分分析实战和可视化 | 附R代码和测试数据
一文看懂PCA主成分分析中介绍了PCA分析的原理和分析的意义(基本简介如下,更多见博客),今天就用数据来实际操练一下.(注意:用了这么多年的PCA可视化竟然是错的!!!) 在公众号后台回复**&quo ...
【Android 逆向】Android 进程注入工具开发 ( 注入代码分析 | 获取注入的 libbridge.so 动态库中的 load 函数地址并通过远程调用执行该函数 )
文章目录一.dlsym 函数简介二.获取目标进程 linker 中的 dlsym 函数地址三.远程调用目标进程 linker 中的 dlsym 函数获取注入的 libbridge.so ...
【Android 逆向】Android 进程注入工具开发 ( 注入代码分析 | 获取 linker 中的 dlopen 函数地址并通过远程调用执行该函数 )
文章目录一.dlopen 函数简介二.获取目标进程 linker 中的 dlopen 函数地址三.远程调用目标进程 linker 中的 dlopen 函数一.dlopen 函数简介 dlo ...

PCA（主成分分析）获取BoundingBox代码分析

引言

理论

代码

结果

分析

参考

PCA（主成分分析）获取BoundingBox代码分析相关推荐

最新文章

热门文章