• 点云法向量计算
千次阅读
2019-07-29 10:46:46

试了PCL中点云的法向量计算：

// 待计算的点云
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_xyz(new pcl::PointCloud<pcl::PointXYZ>)
int cld_sz = cloud_xyz->points.size();
pcl::search::Search<pcl::PointXYZ>::Ptr kdtree(new pcl::search::KdTree<pcl::PointXYZ>());
int search_k = 10;
// 整体一起算
pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne;
ne.setInputCloud(cloud_xyz);
ne.setSearchMethod(kdtree);
ne.setKSearch(search_k);
ne.setViewPoint(std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max());
pcl::PointCloud<pcl::Normal>::Ptr normals_1(new pcl::PointCloud<pcl::Normal>);
ne.compute(*normals_1);
// 单个计算: 用pcl的函数，逐点计算；自己组建协方差矩阵，用Eigen求解计算。
kdtree->setInputCloud(cloud_xyz);
pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne1;
std::vector< std::vector< int > > idx_all(cld_sz);
std::vector< std::vector< float > > dis_all(cld_sz);
std::vector< Eigen::Vector4f > ppara(cld_sz);
std::vector< float > curva(cld_sz);
std::vector< float > difangle(cld_sz), difcurva(cld_sz);
pcl::PointCloud<pcl::Normal>::Ptr normals(new pcl::PointCloud<pcl::Normal>);
normals->points.resize(cld_sz);
float max_dif = -1, max_difc = -10;
for (int i = 0; i < cld_sz; i++)
{
idx_all[i].resize(search_k);
dis_all[i].resize(search_k);
if (kdtree->nearestKSearch(i, search_k, idx_all[i], dis_all[i]) > 0)
{
//计算中心点坐标
double curCenterX = 0, curCenterY = 0, curCenterZ = 0;
//去除自身,索引从1开始
for (int j = 1; j != search_k; ++j)
{
curCenterX += cloud_xyz->points[idx_all[i][j]].x;
curCenterY += cloud_xyz->points[idx_all[i][j]].y;
curCenterZ += cloud_xyz->points[idx_all[i][j]].z;
}
curCenterX = curCenterX / (search_k - 1);
curCenterY = curCenterY / (search_k - 1);
curCenterZ = curCenterZ / (search_k - 1);
//计算协方差矩阵
double xx = 0, xy = 0, xz = 0, yy = 0, yz = 0, zz = 0;
for (int j = 1; j != search_k; ++j)
{
xx += (cloud_xyz->points[idx_all[i][j]].x - curCenterX)*(cloud_xyz->points[idx_all[i][j]].x - curCenterX);
xy += (cloud_xyz->points[idx_all[i][j]].x - curCenterX)*(cloud_xyz->points[idx_all[i][j]].y - curCenterY);
xz += (cloud_xyz->points[idx_all[i][j]].x - curCenterX)*(cloud_xyz->points[idx_all[i][j]].z - curCenterZ);
yy += (cloud_xyz->points[idx_all[i][j]].y - curCenterY)*(cloud_xyz->points[idx_all[i][j]].y - curCenterY);
yz += (cloud_xyz->points[idx_all[i][j]].y - curCenterY)*(cloud_xyz->points[idx_all[i][j]].z - curCenterZ);
zz += (cloud_xyz->points[idx_all[i][j]].z - curCenterZ)*(cloud_xyz->points[idx_all[i][j]].z - curCenterZ);
}
Eigen::Matrix3f covMat(3, 3);
covMat(0, 0) = xx / (search_k - 1);
covMat(0, 1) = covMat(1, 0) = xy / (search_k - 1);
covMat(0, 2) = covMat(2, 0) = xz / (search_k - 1);
covMat(1, 1) = yy / (search_k - 1);
covMat(1, 2) = covMat(2, 1) = yz / (search_k - 1);
covMat(2, 2) = zz / (search_k - 1);
Eigen::EigenSolver<Eigen::Matrix3f> es(covMat);
Eigen::Matrix3f val = es.pseudoEigenvalueMatrix();
Eigen::Matrix3f vec = es.pseudoEigenvectors();
//统计最小特征值及对应的位置,首先默认为第一个特征值为最小值
double t1 = val(0, 0);
unsigned short ii = 0;
//如果第一个特征值比第二个特征值大，则第二个特征值为最小
if (t1 > val(1, 1))
{
ii = 1;
t1 = val(1, 1);
}
//如果第二个特征值比第三个特征值大，则第三个特征值为最小
if (t1 > val(2, 2))
{
ii = 2;
t1 = val(2, 2);
}
Eigen::Vector3f v_n(vec(0, ii), vec(1, ii), vec(2, ii));
v_n = v_n / v_n.norm();
if (v_n(2) < 0) v_n = -v_n;
normals->points[i].normal_x = v_n(0);
normals->points[i].normal_y = v_n(1);
normals->points[i].normal_z = v_n(2);
normals->points[i].curvature = t1 / (val(0, 0) + val(1, 1) + val(2, 2));
ne1.computePointNormal(*cloud_xyz, idx_all[i], ppara[i], curva[i]);
Eigen::Vector3f v_n1(ppara[i](0), ppara[i](1), ppara[i](2));
v_n1 = v_n1 / v_n1.norm();
if (v_n1(2) < 0) v_n1 = -v_n1;
difangle[i] = acos(abs(v_n1.dot(v_n))) * 180 / M_PI;
difcurva[i] = abs(abs(normals->points[i].curvature) - abs(normals_1->points[i].curvature));
if (difangle[i] > max_dif) max_dif = difangle[i];
if (difcurva[i] > max_difc) max_difc = difcurva[i];
}
}

计算后作差，发现PCL整体一起算和单点计算的法向量夹角差在10^-2级别，可以忽略不计，但是自己组建协方差和PCL组建协方差的方式不一样（PCL协方差矩阵组建见源码，相关推导见https://blog.csdn.net/lming_08/article/details/21171491），计算出来的法向夹角差别最大能达到40多°，曲率差能达到0.08，但理论上两种方式都是正确的，为什么差别这么大还未知，如果有大神知道，麻烦留言告知，谢谢。

更多相关内容
• 基于matlab对三维点云法向量进行求取，并进行朝向统一。
• 点云法向量,点云法向量计算原理,matlab源码
• [PCL] 点云法向量计算NormalEstimation 转载来源 从GitHub的代码版本库下载源代码https://github.com/PointCloudLibrary/pcl，用CMake生成VS项目，查看PCL的源码位于pcl_features项目下 1.Feature类： template <...

## [PCL] 点云法向量计算NormalEstimation

转载来源
从GitHub的代码版本库下载源代码https://github.com/PointCloudLibrary/pcl，用CMake生成VS项目，查看PCL的源码位于pcl_features项目下
1.Feature类：

template <typename PointInT, typename PointOutT>   class Feature : public PCLBase<PointInT>


注意 Feature是一个泛型类，有一个compute方法。

template <typename PointInT, typename PointOutT> void pcl::Feature<PointInT, PointOutT>::compute (PointCloudOut &output)
{
if (!initCompute ())
{
output.width = output.height = 0;
output.points.clear ();
return;
}
// Resize the output dataset
if (output.points.size () != indices_->size ())
output.points.resize (indices_->size ());
// Check if the output will be computed for all points or only a subset
// If the input width or height are not set, set output width as size
if (indices_->size () != input_->points.size () || input_->width * input_->height == 0)
{
output.width = static_cast<uint32_t> (indices_->size ());
output.height = 1;
}
else
{
output.width = input_->width;
output.height = input_->height;
}
output.is_dense = input_->is_dense;
// Perform the actual feature computation
computeFeature (output);
deinitCompute ();
}


2.注意computeFeature (output);方法 ，可以知道这是一个私有的虚方法。

private:
/** \brief Abstract feature estimation method.
* \param[out] output the resultant features    */
virtual void    computeFeature (PointCloudOut &output) = 0;


3.查看Feature的继承关系可以知道

template <typename PointInT, typename PointOutT>   class NormalEstimation: public Feature<PointInT, PointOutT>


NormalEstimation类是Feature模板类的子类，因此执行的是NormalEstimation类的computeFeature方法。查看computeFeature方法：

template <typename PointInT, typename PointOutT> void pcl::NormalEstimation<PointInT, PointOutT>::computeFeature (PointCloudOut &output)
{
// Allocate enough space to hold the results
// \note This resize is irrelevant for a radiusSearch ().
std::vector< int> nn_indices (k_);
std::vector< float> nn_dists (k_);
output.is_dense = true;
// Save a few cycles by not checking every point for NaN/Inf values if the cloud is set to dense
if (input_->is_dense)
{
// Iterating over the entire index vector
for (size_t idx = 0; idx < indices_->size (); ++idx)
{
if (this ->searchForNeighbors ((*indices_)[idx], search_parameter_, nn_indices, nn_dists) == 0 ||
!computePointNormal (*surface_, nn_indices, output.points[idx].normal[0], output.points[idx].normal[1], output.points[idx].normal[2], output.points[idx].curvature))
{
output.points[idx].normal[0] = output.points[idx].normal[1] = output.points[idx].normal[2] = output.points[idx].curvature = std::numeric_limits<float >::quiet_NaN ();
output.is_dense = false;
continue;
}

flipNormalTowardsViewpoint (input_->points[(*indices_)[idx]], vpx_, vpy_, vpz_, output.points[idx].normal[0], output.points[idx].normal[1], output.points[idx].normal[2]);
}
}
else
{
// Iterating over the entire index vector
for (size_t idx = 0; idx < indices_->size (); ++idx)
{
if (!isFinite ((*input_)[(*indices_)[idx]]) ||
this->searchForNeighbors ((*indices_)[idx], search_parameter_, nn_indices, nn_dists) == 0 ||
!computePointNormal (*surface_, nn_indices, output.points[idx].normal[0], output.points[idx].normal[1], output.points[idx].normal[2], output.points[idx].curvature))
{
output.points[idx].normal[0] = output.points[idx].normal[1] = output.points[idx].normal[2] = output.points[idx].curvature = std::numeric_limits<float >::quiet_NaN ();

output.is_dense = false;
continue;
}
flipNormalTowardsViewpoint (input_->points[(*indices_)[idx]], vpx_, vpy_, vpz_,  output.points[idx].normal[0], output.points[idx].normal[1], output.points[idx].normal[2]);
}
}
}


4.因此分析NormalEstimation的算法流程如下：
（1）进行点云的初始化initCompute
（2）初始化计算结果输出对象output
（3）计算点云法向量，具体由子类的computeFeature方法实现。先搜索近邻searchForNeighbors ，然后计算computePointNormal
采用的方法是PCA主成分分析法。
参考文献：《Semantic 3D Object Maps for Everyday Manipulation in Human Living Environments》 P45-49
点云的法向量主要是通过点所在区域的局部拟合的表面进行计算。平面通过一个点和法向量进行表示。对于每一个点Pi,对应的协方差矩阵C

关于主成份分析的基本原理和算法流程参考：http://blog.csdn.net/lming_08/article/details/21335313
（4）flipNormalTowardsViewpoint 法向量定向，采用方法是：使法向量的方向朝向viewpoint。
5.NormalEstimation模板类的重载方法computeFeature分析，computePointNormal分析。

inline bool computePointNormal (const pcl::PointCloud<PointInT> &cloud, const std::vector<int> &indices,
float &nx, float &ny, float &nz, float &curvature)
{
if (indices.size () < 3 ||
computeMeanAndCovarianceMatrix (cloud, indices, covariance_matrix_, xyz_centroid_) == 0)
{
nx = ny = nz = curvature = std::numeric_limits<float>::quiet_NaN ();
return false;
}

// Get the plane normal and surface curvature
solvePlaneParameters (covariance_matrix_, nx, ny, nz, curvature);
return true;
}


computeMeanAndCovarianceMatrix主要是PCA过程中计算平均值和协方差矩阵，在类centroid.hpp中。
而solvePlaneParameters方法则是为了求解特征值和特征向量。代码见feature.hpp。具体实现时采用了pcl::eigen33方法。

inline void pcl::solvePlaneParameters (const Eigen::Matrix3f &covariance_matrix,
float &nx, float &ny, float &nz, float &curvature)
{
// Avoid getting hung on Eigen's optimizers
//  for (int i = 0; i < 9; ++i)
//    if (!pcl_isfinite (covariance_matrix.coeff (i)))
//    {
//      //PCL_WARN ("[pcl::solvePlaneParameteres] Covariance matrix has NaN/Inf values!\n");
//      nx = ny = nz = curvature = std::numeric_limits<float>::quiet_NaN ();
//      return;
//    }
// Extract the smallest eigenvalue and its eigenvector
EIGEN_ALIGN16 Eigen::Vector3f::Scalar eigen_value;
EIGEN_ALIGN16 Eigen::Vector3f eigen_vector;
pcl::eigen33 (covariance_matrix, eigen_value, eigen_vector);

nx = eigen_vector [0];
ny = eigen_vector [1];
nz = eigen_vector [2];

// Compute the curvature surface change
float eig_sum = covariance_matrix.coeff (0) + covariance_matrix.coeff (4) + covariance_matrix.coeff (8);
if (eig_sum != 0)
curvature = fabsf (eigen_value / eig_sum);
else
curvature = 0;
}


6.法向量定向

见normal_3d.h文件中，有多个覆写方法。摘其一：

/** \brief Flip (in place) the estimated normal of a point towards a given viewpoint
* \param point a given point
* \param vp_x the X coordinate of the viewpoint
* \param vp_y the X coordinate of the viewpoint
* \param vp_z the X coordinate of the viewpoint
* \param nx the resultant X component of the plane normal
* \param ny the resultant Y component of the plane normal
* \param nz the resultant Z component of the plane normal
* \ingroup features
*/
template <typename PointT> inline void
flipNormalTowardsViewpoint (const PointT &point, float vp_x, float vp_y, float vp_z,
float &nx, float &ny, float &nz)
{
// See if we need to flip any plane normals
vp_x -= point.x;
vp_y -= point.y;
vp_z -= point.z;

// Dot product between the (viewpoint - point) and the plane normal
float cos_theta = (vp_x * nx + vp_y * ny + vp_z * nz);

// Flip the plane normal
if (cos_theta < 0)
{
nx *= -1;
ny *= -1;
nz *= -1;
}
}


运行的实例结果：

展开全文
• MATLAB点云处理：2点云最近点查询和法向量计算、通过索引提取点云 这次主要是点云法向量计算和最近点查询 法向量计算 MATLAB提供了函数pcnormals用于计算点云的法向量，下面用经典的兔子展示一下 clc,clear rabbit...

上一次记录了

MATLAB点云处理1

这次继续写

## MATLAB点云处理：2点云最近点查询和法向量计算、通过索引提取点云

这次主要是点云的法向量计算最近点查询

## 法向量计算

MATLAB提供了函数pcnormals用于计算点云的法向量，下面用经典的兔子展示一下

clc,clear


通过上述代码，点开加载的兔子，可以看到里面是没有Normal属性的

然后通过

rabbit_normals = pcnormals(rabbit)


运行之后，就可以看到rabbit_normals的结果保存了计算出的法向量，这个结果不会保存到rabbit的Normal属性里，但是可以通过

rabbit.Normal = rabbit_normals


直接进行赋值；
为了让法向量也显示出来，参考官网提供的示例，可以通过quiver函数实现

rabbit = pcread('rabbit.pcd')
rabbit_normals = pcnormals(rabbit)
rabbit.Normal = rabbit_normals;
pcshow(rabbit)
hold on
x = rabbit.Location(:,1);
y = rabbit.Location(:,2);
z = rabbit.Location(:,3);
u = rabbit.Normal(:,1);
v = rabbit.Normal(:,2);
w = rabbit.Normal(:,3);
quiver3(x,y,z,u,v,w);


最终显示结果如下图所示

由于点云法向量的计算是通过找到临近点拟合平面计算法向量得来的，这就涉及到临近点数量的选取，pcnormals还提供了一个参数的设置，例如

rabbit_normals = pcnormals(rabbit,10)


这就是选取了当前点最近的10个点来计算法向量

## 最近点查询

刚才说到在计算法向量的时候，其实利用到最近点查询，那么就关注一下，MATLAB是怎么实现最近点查询的？

[indices,dists] = findNearestNeighbors(ptCloud,point,K)


这是这两种函数最普通也是最常用的调用方式
返回的indices是目标点云中的对应点的索引，dists是对应的距离，下面通过兔子展示一下：

clc,clear
% rabbit_normals = pcnormals(rabbit,10)
% rabbit.Normal = rabbit_normals;
current_point_indic = unidrnd(rabbit.Count);
current_point = rabbit.Location(current_point_indic,:)
[indices, dists] = findNearestNeighbors(rabbit,current_point,50);
pcshow(rabbit)
hold on
plot3(current_point(1),current_point(2),current_point(3),'*r')
plot3(rabbit.Location(indices,1),rabbit.Location(indices,2),rabbit.Location(indices,3),'*')


最终的结果如下，可以看到兔子的身体黄了一坨

半径范围查询代码和结果如下

clc,clear
% rabbit_normals = pcnormals(rabbit,10)
% rabbit.Normal = rabbit_normals;
current_point_indic = unidrnd(rabbit.Count);
current_point = rabbit.Location(current_point_indic,:)
pcshow(rabbit)
hold on
plot3(current_point(1),current_point(2),current_point(3),'*r')
plot3(rabbit.Location(indices,1),rabbit.Location(indices,2),rabbit.Location(indices,3),'*')


这里可以看到，由于点的距离比较近，这次直接黄了一大块

## 通过索引提取点云

刚才通过最近点搜索返回的是点云里对应的索引，我们是通过复制的方式，比较笨的提取需要的点，其实MATLAB中提供了函数select函数

ptCloudOut = select(ptCloud,indices)


展开全文
• 基于局部平面拟合求点云法向量，是我看到过的比较容易懂的，算法推导也写得很清楚
• 因此分析NormalEstimation的算法流程如下： （1）进行点云的初始化initCompute　（2）初始化计算结果输出对象output　（3）计算点云法向量，具体由子类的computeFeature方法实现。先搜索近邻searchForNeighbors ，...

从GitHub的代码版本库下载源代码https://github.com/PointCloudLibrary/pcl，用CMake生成VS项目，查看PCL的源码位于pcl_features项目下1.Feature类：template <typename PointInT, typename PointOutT> class Feature : public PCLBase注意 Feature是一个泛型类，有一个compute方法。 1 template <typename PointInT, typename PointOutT> void pcl::Feature<PointInT, PointOutT>::compute (PointCloudOut &output)
2 {
3 if (!initCompute ())
4 {
5 output.width = output.height = 0;
6 output.points.clear ();
7 return;
8 }
11 // Resize the output dataset
12 if (output.points.size () != indices_->size ())
13 output.points.resize (indices_->size ());
14 // Check if the output will be computed for all points or only a subset
15 // If the input width or height are not set, set output width as size
16 if (indices_->size () != input_->points.size () || input_->width * input_->height == 0)
17 {
18 output.width = static_cast<uint32_t> (indices_->size ());
19 output.height = 1;
20 }
21 else
22 {
23 output.width = input_->width;
24 output.height = input_->height;
25 }
26 output.is_dense = input_->is_dense;
27 // Perform the actual feature computation
28 computeFeature (output);
29 deinitCompute ();
30 }2.注意computeFeature (output);方法 ，可以知道这是一个私有的虚方法。 private: /** \brief Abstract feature estimation method. * \param[out] output the resultant features */ virtual void computeFeature (PointCloudOut &output) = 0;3.查看Feature的继承关系可以知道template <typename PointInT, typename PointOutT> class NormalEstimation: public Feature<PointInT, PointOutT>NormalEstimation类是Feature模板类的子类，因此执行的是NormalEstimation类的computeFeature方法。查看computeFeature方法： 1 template <typename PointInT, typename PointOutT> void pcl::NormalEstimation<PointInT, PointOutT>::computeFeature (PointCloudOut &output)
2 {
3 // Allocate enough space to hold the results
4 // \note This resize is irrelevant for a radiusSearch ().
5 std::vector< int> nn_indices (k_);
6 std::vector< float> nn_dists (k_);
7 output.is_dense = true;
8 // Save a few cycles by not checking every point for NaN/Inf values if the cloud is set to dense
9 if (input_->is_dense)
10 {
11 // Iterating over the entire index vector
12 for (size_t idx = 0; idx < indices_->size (); ++idx)
13 {
14 if (this ->searchForNeighbors ((*indices_)[idx], search_parameter_, nn_indices, nn_dists) == 0 ||
15 !computePointNormal (*surface_, nn_indices, output.points[idx].normal[0], output.points[idx].normal[1], output.points[idx].normal[2], output.points[idx].curvature))
16 {
17 output.points[idx].normal[0] = output.points[idx].normal[1] = output.points[idx].normal[2] = output.points[idx].curvature = std::numeric_limits::quiet_NaN ();
18 output.is_dense = false;
19 continue;
20 }
21
22 flipNormalTowardsViewpoint (input_->points[(*indices_)[idx]], vpx_, vpy_, vpz_, output.points[idx].normal[0], output.points[idx].normal[1], output.points[idx].normal[2]);
23 }
24 }
25 else
26 {
27 // Iterating over the entire index vector
28 for (size_t idx = 0; idx < indices_->size (); ++idx)
29 {
30 if (!isFinite ((*input_)[(*indices_)[idx]]) ||
31 this->searchForNeighbors ((*indices_)[idx], search_parameter_, nn_indices, nn_dists) == 0 ||
32 !computePointNormal (surface_, nn_indices, output.points[idx].normal[0], output.points[idx].normal[1], output.points[idx].normal[2], output.points[idx].curvature))
33 {
34 output.points[idx].normal[0] = output.points[idx].normal[1] = output.points[idx].normal[2] = output.points[idx].curvature = std::numeric_limits::quiet_NaN ();
35
36 output.is_dense = false;
37 continue;
38 }
39 flipNormalTowardsViewpoint (input_->points[(indices_)[idx]], vpx_, vpy_, vpz_, output.points[idx].normal[0], output.points[idx].normal[1], output.points[idx].normal[2]);
40 }
41 }
42 }4.因此分析NormalEstimation的算法流程如下： 　 （1）进行点云的初始化initCompute　　（2）初始化计算结果输出对象output　　（3）计算点云法向量，具体由子类的computeFeature方法实现。先搜索近邻searchForNeighbors ，然后计算computePointNormal　　　　采用的方法是PCA主成分分析法。　　　　参考文献：《Semantic 3D Object Maps for Everyday Manipulation in Human Living Environments》 P45-49　　　　　　　　　　http://www.pclcn.org/study/shownews.php?lang=cn&id=105　　　　点云的法向量主要是通过点所在区域的局部拟合的表面进行计算。平面通过一个点和法向量进行表示。对于每一个点Pi,对应的协方差矩阵C　　　　　　　　关于主成份分析的基本原理和算法流程参考：http://blog.csdn.net/lming_08/article/details/21335313　　（4）flipNormalTowardsViewpoint 法向量定向，采用方法是：使法向量的方向朝向viewpoint。5.NormalEstimation模板类的重载方法computeFeature分析，computePointNormal分析。 1 inline bool computePointNormal (const pcl::PointCloud &cloud, const std::vector &indices,
2 float &nx, float &ny, float &nz, float &curvature)
3 {
4 if (indices.size () < 3 ||
5 computeMeanAndCovarianceMatrix (cloud, indices, covariance_matrix_, xyz_centroid_) == 0)
6 {
7 nx = ny = nz = curvature = std::numeric_limits::quiet_NaN ();
8 return false;
9 }
10
11 // Get the plane normal and surface curvature
12 solvePlaneParameters (covariance_matrix_, nx, ny, nz, curvature);
13 return true;
14 }computeMeanAndCovarianceMatrix主要是PCA过程中计算平均值和协方差矩阵，在类centroid.hpp中。
而solvePlaneParameters方法则是为了求解特征值和特征向量。代码见feature.hpp。具体实现时采用了pcl::eigen33方法。 1 inline void pcl::solvePlaneParameters (const Eigen::Matrix3f &covariance_matrix,
2 float &nx, float &ny, float &nz, float &curvature)
3 {
4 // Avoid getting hung on Eigen’s optimizers
5 // for (int i = 0; i < 9; ++i)
6 // if (!pcl_isfinite (covariance_matrix.coeff (i)))
7 // {
8 // //PCL_WARN ("[pcl::solvePlaneParameteres] Covariance matrix has NaN/Inf values!\n");
9 // nx = ny = nz = curvature = std::numeric_limits::quiet_NaN ();
10 // return;
11 // }
12 // Extract the smallest eigenvalue and its eigenvector
13 EIGEN_ALIGN16 Eigen::Vector3f::Scalar eigen_value;
14 EIGEN_ALIGN16 Eigen::Vector3f eigen_vector;
15 pcl::eigen33 (covariance_matrix, eigen_value, eigen_vector);
16
17 nx = eigen_vector [0];
18 ny = eigen_vector [1];
19 nz = eigen_vector [2];
20
21 // Compute the curvature surface change
22 float eig_sum = covariance_matrix.coeff (0) + covariance_matrix.coeff (4) + covariance_matrix.coeff (8);
23 if (eig_sum != 0)
24 curvature = fabsf (eigen_value / eig_sum);
25 else
26 curvature = 0;
27 } 6.法向量定向见normal_3d.h文件中，有多个覆写方法。摘其一： 1 /
\brief Flip (in place) the estimated normal of a point towards a given viewpoint
2 * \param point a given point
3 * \param vp_x the X coordinate of the viewpoint
4 * \param vp_y the X coordinate of the viewpoint
5 * \param vp_z the X coordinate of the viewpoint
6 * \param nx the resultant X component of the plane normal
7 * \param ny the resultant Y component of the plane normal
8 * \param nz the resultant Z component of the plane normal
9 * \ingroup features
10 */
11 template inline void
12 flipNormalTowardsViewpoint (const PointT &point, float vp_x, float vp_y, float vp_z,
13 float &nx, float &ny, float &nz)
14 {
15 // See if we need to flip any plane normals
16 vp_x -= point.x;
17 vp_y -= point.y;
18 vp_z -= point.z;
19
20 // Dot product between the (viewpoint - point) and the plane normal
21 float cos_theta = (vp_x * nx + vp_y * ny + vp_z * nz);
22
23 // Flip the plane normal
24 if (cos_theta < 0)
25 {
26 nx *= -1;
27 ny *= -1;
28 nz *= -1;
29 }
30 }运行的实例结果：

展开全文
• 最小二乘计算点云法向量原理 对于任意一点 p(x,y,z)p(x, y, z)p(x,y,z)查找其一定领域内的点{pi}\{{p_i}\}{pi​}； 求得一个平面 ∏:a0x+a1y+a2z+1=0\prod: a_0x + a_1y+a_2z + 1 = 0∏:a0​x+a1​y+a2​z+1=0使得其...
• 利用Matlab计算三维点云中每个点的单位法向量与特征值。 通过计算当前点与其邻域点的协方差矩阵，进而求得法向量
• ## 点云法向量

千次阅读 2022-04-23 11:02:13
python open3d 计算法向量
• 核心代码 void pcaNormal(const PointCloud<PointXYZ>::Ptr& cloud, const int& k_neighbors, PointCloud<Normal>::Ptr& normals) //cloud为所求取的点云，k_neighbors为输入的Kd树的k值，normals为点云法向量 { /...
• 法向量是点云的一个重要特征，我们可以应用点云法向量到平面分割、模板匹配目标检测等领域中。下面介绍如何在PCL中计算法向量。 3.数据 本例中使用的点云数据（test.pcd）请见百度网盘分享。 链接：...
• 点云法向量与点云平面拟合的关系（PCA） Estimating Surface Normals in a PointCloud 3D【24】PCA点云法向量估计 利用PCA计算点云的法线 3D点云法向量估计(最小二乘拟合平面) 为什么用PCA做点云法线估计？ ...
• ## PCL 计算点云法向量并显示

万次阅读 多人点赞 2020-04-16 16:11:27
PCL计算法向量并显示，用OMP进行计算加速以及法线定向的深层次理解、操作。
• Meshlab是一个强大的三维模型处理的软件，在三维点云的处理方面也十分强大，不仅仅可以用三维点云重构出三维模型，还可以...
• import pcl import pcl.pcl_visualization def main(): cloud = pcl.load('filename.pcd') # load the pcd file ne = cloud.make_NormalEstimation() # instance the NormalEstimation ne.set_KSearch(20) # ...
• ## 计算点云的法向量

千次阅读 2022-04-08 15:26:54
用pcl::NormalEstimation简直就是坑爹，计算出的点云法向量有40~50%都是有问题的 pcl::search::KdTree::Ptr kdtree_submap(new pcl::search::KdTree); kdtree_submap->setInputCloud(cloud_submap);// Make sure the...
• 该代码用于点云数据中每一点的法向量的估算，计算速度较快，带测试数据
• ## pcl计算点云法向量

千次阅读 2020-12-03 20:46:16
最近因为项目，需要计算点云的法向量，所以在网上看了一些资料，然后知道pcl库里面有这些功能，pcl的法向量计算的原理： pcl里面计算点云（自己的理解） 根据顶点采样最近的局部点云(k个)，根据自己的点云拟合出一个...
• 目录 一、法向量 1、概述 2、主要函数 2.1 计算法线 2.2 自定义方向定向 2.3 与相机位置定向 2.4 最小生成树定向 二、代码实现 三、结果展示 1、自定义方向定向 2、与相机位置定向 3、最小代价生成树定向 四、法线...
• 对于三维点云的每一个点利用邻域点拟合最小二乘平面并求得法向量
• 该博客主要进行： （1）点云读取可视化 （2）下采样可视化 （3）法向量三种估计方式（K近邻估计，半径近邻估计，混合搜索估计） （4）点云法向量可视化 （5）点云每个点对应的法向量点存储 （6）点云法向量点可视化
• Output:主成分向量 ，principle vectors ，仅仅是一个向量，代表一个方向。主成分的个数k<=原始空间维度数d Q&A 什么是最主要的成分？ 点的投影后分布的方差最大的方向。(在该方向上的点云分布的最分散) ...
• 文章目录 1 点云法向量（法线）估计 2 法线估计与可视化 3 OpenMP加速法线估计 1 点云法向量（法线）估计 PCL法线估计是通过寻找待估计点与其邻域点共同构成的局部平面的法向量实现的。通过主成分分析，最小特征向量...
• ## 3D【24】PCA点云法向量估计

万次阅读 多人点赞 2018-06-05 19:50:09
基于PCA的点云法向量估计，其实是从最小二乘法推导出来的。假设我们要估计某一点的法向量，我们需要通过利用该点的近邻点估计出一个平面，然后我们就能计算出该点的法向量。或者可以这么说，通过最小化一个目标函数...
• 计算点云法向量的夹角
• 法向量为投影到某方向后，信息量最小的方向，因此需要PCA变换，将点云投影到特征值最小的方向。由于需要PCA算法，因此先把昨天的PCA算法拷贝过来。` import open3d as o3d import numpy as np from pyntcloud ...
• 上代码： #include #include #include #include #include #include #include int main(int argc, char** argv) { //读取点云 pcl::PointCloud::Ptr cloud(new pcl::PointCloud); if (pcl::io::loadPCDFile("E://1PCD...

...