利用PCA计算点云的法线

我们知道PCA可以用来降维,并使降维后的数据尽可能保持原来的特征。比如二维散乱的点：

经过PCA降维后,变成了一维直线,而该直线保证点尽可能分散,变成如下图(跟最小二乘是一样的)：

具体原理可参考点击打开链接

前面说的是二维降到一维时的情况，假如我们有一堆散乱的三维点云,则可以这样计算法线：

1）对每一个点，取临近点，比如取最临近的50个点，当然会用到K-D树

        2）对临近点做PCA降维，把它降到二维平面上,可以想象得到这个平面一定是它的切平面(在切平面上才可以尽可能分散）

        3）切平面的法线就是该点的法线了，而这样的法线有两个，取哪个还需要考虑临近点的凸包方向，这里就不叙述了，详情查看：点击打开链接

顺便贴上我做三维直线提取的代码，跟最上面的提取二维线一样：

std::vector<Eigen::Vector3d> LeastSquare::Generate() const{std::vector<Eigen::Vector3d> res;std::size_t count = _vertexs.size();double xmean = 0.0, ymean = 0.0;for (int i = 0; i < count; i++){xmean += _vertexs.at(i)(0);ymean += _vertexs.at(i)(1);}xmean = xmean / count;ymean = ymean / count;Eigen::MatrixXd a(2, count);for (int i = 0; i < count; i++){a(0, i) = _vertexs.at(i)(0) - xmean;a(1, i) = _vertexs.at(i)(1) - ymean;}Eigen::Matrix2d c = a * a.transpose() * 1 / count;Eigen::EigenSolver<Eigen::MatrixXd> es(c);std::complex<double> ev[2];ev[0] = es.eigenvalues()[0];ev[1] = es.eigenvalues()[1];int iMax = 0; double dMin = -DBL_MAX;for (int i = 0; i < 2; i++){if (ev[i].real() > dMin){iMax = i;dMin = ev[i].real();}}Eigen::VectorXcd V1 = es.eigenvectors().col(iMax);V1.normalized();Eigen::VectorXcd V2 = iMax == 0 ? es.eigenvectors().col(1) : es.eigenvectors().col(0);V2.normalized();Eigen::MatrixXd YY = V1.real().transpose() * a;Eigen::MatrixXd Y(2, count);for (int i = 0; i < count; i++){Y(0, i) = YY(0, i);Y(1, i) = 0.0;}Eigen::Matrix2d P;P(0, 0) = V1(0).real();P(0, 1) = V1(1).real();P(1, 0) = V2(0).real();P(1, 1) = V2(1).real();Eigen::MatrixXd xyf = P.inverse() * Y;for (int i = 0; i < count; i++){xyf(0, i) = xyf(0, i) + xmean;xyf(1, i) = xyf(1, i) + ymean;}for (int i = 0; i < count; i++){res.push_back(Eigen::Vector3d(xyf(0, i), xyf(1, i), 0.0));}return res;}