估算点云表面法向
表面法线是几何表面的重要性质,在许多领域,如计算机图形应用程序中被广泛应用,用于产生光源的阴影和其他视觉效果。
给定一个几何曲面,通常用垂直于该点的向量来推断曲面上某一点法线的方向是琐细的。然而,由于我们获取的点云数据集代表了真实表面上的一组点样本,所以有两种可能性:
- 利用曲面网格划分技术,从获取的点云数据集中获取下垫面,然后从网格中计算曲面法线;
- 利用近似法直接从点云数据集中推断表面法线。
本教程将讨论后者,即给定一个点云数据集,直接计算云中每个点的表面法线。
理论基础
尽管存在许多不同的常规估计方法,但我们将在本教程中重点介绍的方法是最简单的方法之一,其公式如下。确定曲面上某一点法线的问题近似于估计与曲面相切的平面法线的问题,进而成为一个最小二乘平面拟合估计问题。
注意
为获得更多信息,包括数学方程的最小二乘问题,请查阅Rusu论文
因此,估计表面法线的解决方案被简化为对由查询点的最近邻点创建的协方差矩阵的特征向量和特征值(或PCA -主成分分析)进行分析。更具体地说,对于每个点
p
i
p_i
pi,我们组装协方差矩阵
C
C
C如下:
C
=
1
k
⋅
∑
i
=
0
k
(
p
i
−
p
ˉ
)
⋅
(
p
i
−
p
ˉ
)
T
,
C
⋅
v
→
j
=
λ
j
⋅
v
→
j
,
j
∈
0
,
1
,
2
C = \frac{1}{k}\cdot\sum_{i=0}^k(p_i-\bar p)\cdot(p_i-\bar p)^T, C\cdot \overrightarrow v_j = \lambda_j\cdot \overrightarrow v_j, j\in{0, 1, 2}
C=k1⋅i=0∑k(pi−pˉ)⋅(pi−pˉ)T,C⋅vj=λj⋅vj,j∈0,1,2
其中
k
k
k为在
p
i
p_i
pi邻域内考虑的点数,
p
ˉ
\bar p
pˉ表示最近邻的三维质心,
λ
j
\lambda_j
λj为协方差矩阵的第
j
j
j个特征值,
v
→
j
\overrightarrow v_j
vj为第
j
j
j个特征向量。
要从PCL中的一组点估计协方差矩阵,可以使用:
// 在每个表面块上的3x3协方差矩阵
Eigen::Matrix3f covariance_matrix;
// 16字节对齐的占位符,用于表面块的XYZ质心
Eigen::Vector4f xyz_centroid;
// 估计XYZ质心
compute3DCentroid (cloud, xyz_centroid);
// 计算3x3协方差矩阵
computeCovarianceMatrix (cloud, xyz_centroid, covariance_matrix);
一般来说,由于没有数学方法来解法向准确值,因此上面所示的通过主成分分析(PCA)计算出的法线方向是模糊的,并且在整个点云数据集上的方向并不一致。下图显示了在厨房环境一部分的较大数据集中的两个部分上的这些效果。图的右侧是扩展高斯图像(Extended Gaussian Image, EGI),也称为法向球,它描述了从点云出发的所有法线的方向。由于数据集是2.5D的,因此是从一个单一的角度获取的,因此在EGI中,法线应该只出现在球体的一半上。然而,由于方向不一致,它们分布在整个球面上。
如果视图
V
p
V_p
Vp实际上是已知的,那么这个问题的解决方案就是简单的。为了使所有法线
n
→
j
\overrightarrow n_j
nj始终指向视点,需要满足以下条件:
n
→
j
⋅
(
V
p
−
p
i
)
>
0
\overrightarrow n_j\cdot (V_p-p_i) > 0
nj⋅(Vp−pi)>0
下图显示的结果是,来自上图的数据集中的所有法线都一致指向视点之后的结果。
要在PCL中手动重新定位给定的法线点,可以使用:
flipNormalTowardsViewpoint (const PointT &point, float vp_x, float vp_y, float vp_z, Eigen::Vector4f &normal);
注意
如果数据集具有多个采集视点,则上述常规的重定向方法不成立,需要实现更复杂的算法。更多信息请参见Rusu论文。
选择合适的比例
如前所述,需要从点的周围临近点(也称为k邻域)估计点处的表面法线。
近邻估计的具体问题提出合适的比例因子的问题:给定一个采样点云数据集,什么是正确的k(鉴于通过pcl::Feature::setKSearch)或者r(鉴于通过pcl::Feature::setRadiusSearch)值应该在决定使用最近的邻居的一个点的集合?
这个问题非常重要,是自动评估(即用户没有给定阈值)的点特征表示。为了更好地说明这个问题,下图展示了选择较小尺度(即,小r或k)与大尺度。左边的图描绘了一个合理的精心选择的比例因子,两个平面的估计表面法线近似垂直,整个桌子上的小边都可见。如果缩放系数太大(右)的一部分,因此从相邻的邻居更大的集合覆盖点表面,估计特征点表征得到扭曲,旋转曲面法线在边缘的两个平面的表面,并涂抹边缘和隐含的细节。
在不涉及太多细节的情况下,我们可以假设,目前必须根据应用程序所需的详细程度来选择确定点的邻域的范围。简单地说,如果杯子把手和圆柱形部分之间的边缘曲率很重要,那么比例因子需要足够小才能捕捉到这些细节,否则就需要足够大。
估计法线
虽然在特性中已经给出了一个正常估计的示例,但是我们将在这里修改其中一个示例,以便更好地解释幕后发生了什么。
下面的代码片段将为输入数据集中的所有点估计一组表面法线。
#include <pcl/point_types.h>
#include <pcl/features/normal_3d.h>
{
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
... read, pass in or create a point cloud ...
// 创建normal estimate类,并将输入数据集传递给它
pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne;
ne.setInputCloud (cloud);
// 创建一个空的kdtree表示,并将其传递给普通的评估对象。
// 它的内容将根据给定的输入数据集填充到对象中(因为没有其他搜索表面)。
pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ> ());
ne.setSearchMethod (tree);
// 输出数据集
pcl::PointCloud<pcl::Normal>::Ptr cloud_normals (new pcl::PointCloud<pcl::Normal>);
// 使用半径为3cm的球体中的所有相邻点
ne.setRadiusSearch (0.03);
// 计算特征
ne.compute (*cloud_normals);
// cloud_normals->points.size () 应该 cloud->points.size ()的大小一样
}
normalestimate类的实际计算调用在内部什么都不做,但是:
对于点云中的每一点p
1. 得到p的最近邻近点
2. 计算p的曲面法线n
3.检查n是否始终朝向视点,否则翻转
视点默认为(0,0,0),可以用以下命令进行更改:
setViewPoint (float vpx, float vpy, float vpz);
若要计算单点法线,请使用:
computePointNormal (const pcl::PointCloud<PointInT> &cloud, const std::vector<int> &indices, Eigen::Vector4f &plane_parameters, float &curvature);
点云指输入包含点,索引代表的点云的k近邻,平面参数和曲率代表法向的正常估计,平面参数包含法向(nx、ny、nz),和第四个坐标是D中心点。输出曲面曲率估计为协方差矩阵特征值之间的关系(如上所示),为,
σ = λ 0 λ 0 + λ 1 + λ 2 \sigma=\frac{\lambda_0}{\lambda_0+\lambda_1+\lambda_2} σ=λ0+λ1+λ2λ0
使用OpenMP加速法向估计
对于速度敏感的用户,PCL提供了一个额外的表面法向估计,它使用使用OpenMP的多核/多线程范例来加速计算。类的名称是pcl::NormalEstimationOMP,它的API与单线程pcl:: normalestimate 100%兼容,这使它适合作为drop-in替换。在一个有8个内核的系统上,您应该可以获得6-8倍的计算速度。
注意
如果您的数据集是有序的(例如,使用TOF相机、立体声相机等获得的,也就是说,它有一个宽度和一个高度),要获得更快的结果,请查看使用积分图像的法向估计。