表面点云法线

表面法线是几何表面的重要属性，在许多领域（例如计算机图形应用程序）中大量使用，以应用正确的光源以产生阴影和其他视觉效果。

1. 使用表面网格化技术从获取的点云数据集中获取基础表面，然后从网格中计算表面法线；
2. 使用近似值直接从点云数据集中推断表面法线。
   通常我们一般使用后者，即给定点云数据集，直接计算云中每个点的表面法线。

基础理论

尽管有许多不同的法线估计方法，我们先了解其中最简单也是最常见的一个，确定表面一点法线的问题近似于估计表面的一个相切面法线的问题，因此转换过来以后就变成一个最小二乘法平面拟合估计问题。

基本概念

1. $p_i$:几何坐标下的一个点
2. $<p_i,p_j>$: 两个几何点向量的点乘
3. $p_n$:一个几何点$p_i$的表面法向量评估，有坐标 { n x , n y , n z } \{n_x,n_y,n_z\} {nx​,ny​,nz​}的法向量
4. P n = { p 1 , p 2 , . . . . } P_n=\{p_1,p_2,....\} Pn​={p1​,p2​,....}:含有n个点的集合，使用P作为代表
5. $P^k$:点$p_j$($j\le k$)的集合，这个集合在$p_i$的k临近范围内。
6. r:确定$P^k$临近范围的半径
7. $\bar{p}$:一个点集的均值
   KaTeX parse error: Got function '\sum' with no arguments as argument to '\dot' at position 25: …rac{1}{k} \dot \̲s̲u̲m̲_{i=1}^{k} p_i
8. $\vec{z}$ 一个被给予的坐标系统的z轴。

法向量计算

• $p_q$一个平面上的一个点，他的周围存在k个临近点，我们称之为$p^k$。

• 通过这个k个临近点加上$q_q$，我们可以估量平面的几何特征。其中表面法向量是一个重要的特征，再许多领域多有重要的应用。

• 存在很多的表面法向量的计算方法。我们介绍一种简单一定点的方法。

• 设平面最终会被评估为一个点$x$和一个表面法向量$\vec{n}$

• 如何评估$x$和$\vec{n}$，我们设定一个距离$d_i$，设这个距离$d_i=(p_i-x)\vec{n}$,我们的目标通过最小二乘法来最小化这个$d_i$(实际上是PCA降维方法),先设置以下内容
  $$x=\bar{p}=\frac{1}{k}\cdot \sum _{i=1}^k{p}_i$$

• 然后我们通过以下公式获得C，然后获取特征向量与特征值。
  $$C=\frac{1}{k}\sum _{i=1}^k\zeta (p_i-\bar{p})\cdot (p_i-\bar{p})$$

• 这里的$\zeta$是一个定值，通常会取1

• 对以上$C$做特征分解，获得特征值，由于C是对称向量且是半正定向量d，因此他的特征值都是大于等于0的。

• 这里再说明下，k是点的数目，$\bar{p}$是周围点的质心，${\lambda }_j$是协方差矩阵的第j个特征值，${\vec{v}}_j$是协方差矩阵的第j个特征向量。通过以下代码可以计算质心和协方差矩阵(PCL代码)

// 定义一小块表面区域的协方差矩阵
Eigen::Matrix3f covariance_matrix;
// 定义一小块表面区域的质心坐标
Eigen::Vector4f xyz_centroid;
// 计算质心坐标
pcl::compute3DCentroid(*cloud, xyz_centroid);
// 计算3x3的协方差矩阵
pcl::computeCovarianceMatrix(*cloud , xyz_centroid, covariance_matrix);

• if $0\le {\lambda }_0\le {\lambda }_1\le {\lambda }_2$
• ${\lambda }_0$对应的特征向量$\vec{n}$(或者是$-\vec{n}$)就是法向量.
• 如果法向量$\vec{n}$是单位向量那么我们可以用($\varphi ,\theta$),来表示这个法向量
  $$\varphi =arctan\frac{n_z}{n_y},\theta =arctan\frac{\sqrt{(n_x^2+n_y^2)}}{n_x}$$

那么到现在为止，还有一个问题没有被确定，那就是法向量的方向(正负方向，还没办法确定下来)

法向量的方向问题

• 数学上，对这个法向量方向问题没有实际解的方法，但可以从其他方面考虑。 论文45页 考虑了两种情况，见图4.8与4.9
• 下面是图4.8是没有向量正负方向没有调整过的情况 ,从中间图以及左侧那个图我们可以发现中间的大量的法向量是不对的
  

1. 存在视点，非常容易被解决，只要确保以下公式就可以了
   $${\vec{n}}_i\cdot (v_p-p_i)>0$$
   下面的图4.9是经过视点调整的，明显好了很多。
   2. 没法获得视点的值，那就可能要麻烦一点了
   1. 首先介绍下什么是最小生成树，欧式最小生成树，是以欧式距离为权重的最小生成树。
   2. 我们认为在一个光滑表面并且几何封闭中两个点$p_i和p_j$，它们的法向量应该满足以下公式(至少接近)
      $${\vec{n}}_i\cdot {\vec{n}}_j\approx 1$$
   3. 以上结论在稠密采样数据集中基本都是正确的
   • 我们认为在$p_i\in P$的情况下，所有的$p_i$都可以被构建为一个图。以数据点为点，以两点间的距离为权重，可以构建一个欧拉距离的最小生成树，我们需要在$P$中找一个法线方向正确的点$p_i$,然后遍历每一个做小生成树中的临近点，确保完成以下公式
     $$if{\vec{n}}_i\cdot {\vec{n}}_j\le 0\Rightarrow {\vec{n}}_i=-{\vec{n}}_j$$

• 这边欧拉最小生成树的边的权重有时导致失败在传播方向上，作者提出了一种新的边权重方案
  $$e_g=1-|{\vec{n}}_i\cdot {\vec{n}}_j|$$
• 以上权重代价，可以更好的遍历最小生成树，如果一些近似平行的法向量，以上值就越是接近于0，

曲率评估

还记得前面介绍的PCA方法吗，以$C$的法向量为主我们定义${\lambda }_0=min({\lambda }_i)$ ,而这个${\lambda }_0$对应特征向量的就是法向量,我们使用协方差矩阵C的${\lambda }_i$作为一种表面变化的近似。
$${\sigma }_p=\frac{{\lambda }_0}{{\lambda }_0+{\lambda }_1+{\lambda }_2}$$
以上这个比率表示点p的k临近$P^k$的曲面变化。这个变化是抗缩放的。如果${\sigma }_p$是个小值，就表明$P^k$中的所有点在表面的切平面上。
${\sigma }_{p_i}$主要取决于临近点K和协方差矩阵C，并且很容收到噪声的影响。为了提高抗噪声的能力，作者提出了一种抗噪声的方法。
这种方法其实就是将点分成外层点和内层点，区分方法就是将$P^k$领域内的点求平均值，然后设定最大距离$d_{max}$，如果点到平均点的距离比$d_{max}$大，我们称之为外层点，反之我们称之为内层点。然后通过以下公式调整常量$\zeta$

这里的$\mu$是检索点p到他所有的k临近点$p_i^k\in P^k$的平均值，和$d_i$检索点p到它某一个k临近点的具体距离。

选择合适的比例

如前所述，需要根据该点的周围点邻域支持，也称为k-neighborhood（k邻域）来估算该点的表面法线。 最近邻估计问题的细节提出了正确的比例因子的问题：给定采样点云数据集，如何选择正确的k（通过pcl::Feature::setKSearch给出）或r（通过pcl::Feature::setRadiusSearch 给出）值来确定点的最近邻居？

这个问题非常重要，并且构成了点特征表示的自动估计（即在没有用户给定阈值的情况下）的限制因素。为了更好地说明此问题，下图显示了选择较小的比例（即较小的r或k）与较大的比例（即较大的r或k）的效果。图的左侧部分描绘了一个合理选择的比例因子，其中估计的表面法线大致垂直于两个平面，并且在整个桌子上可见小的边缘。但是，如果比例因子太大（右侧部分），则相邻对象的集合会覆盖来自相邻表面的较大点，则估计的点要素表示会失真，在两个平面边缘处旋转的曲面法线会被涂抹边缘和压制的精细细节。

无需赘述太多，只需假设现在必须根据应用所需的详细程度来选择确定点的邻域的比例即可。简而言之，如果杯子的手柄和圆柱部分之间的边缘处的曲率很重要，则比例因子必须足够小以捕获这些细节，否则要大。

默认视点坐标为(0,0,0)，可以使用以下代码修改：

setViewPoint(float vpx, float vpy, float vpx);

计算表面法向量内部伪代码：

// 遍历每个点云P中的点p
for each point p in cloud P 
    // 得到p点的最近邻
    1. get the nearest neighbors of p
    // 计算p点的表面法线n
    2. compute the surface normal n of p
    // 检查n的方向是否指向视点，如果不是则进行反转
    3. check if n is consistently oriented towards the viewpoint and flip otherwise

代码实现

#include <pcl/visualization/cloud_viewer.h>
#include <iostream>
#include <pcl/io/io.h>
#include <pcl/io/pcd_io.h>
#include <pcl/features/normal_3d.h>

/**
 * 评估法向量
 */
int main(){
    // load point cloud
    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);
    pcl::io::loadPCDFile("./data/target.pcd", *cloud);

    // estimate normals     
    //法向量评估器
    pcl::PointCloud<pcl::Normal>::Ptr normals(new pcl::PointCloud<pcl::Normal>);

    // Object for normal estimation.
    pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> normalEstimation;


    //normalEstimation.setIndices()
    normalEstimation.setInputCloud(cloud);

    // For every point, use all neighbors in a radius of 3cm.
    normalEstimation.setRadiusSearch(0.03);

    // A kd-tree is a data structure that makes searches efficient. More about it later.
    // The normal estimation object will use it to find nearest neighbors.
    pcl::search::KdTree<pcl::PointXYZ>::Ptr kdtree(new pcl::search::KdTree<pcl::PointXYZ>);
    normalEstimation.setSearchMethod(kdtree);
    // Calculate the normals.
    normalEstimation.compute(*normals);

    // visualize normals
    pcl::visualization::PCLVisualizer viewer("PCL Viewer");
    viewer.setBackgroundColor(0.0, 0.0, 0.5);
    viewer.addPointCloud<pcl::PointXYZ>(cloud, "cloud");

    // 参数int level=2 表示每n个点绘制一个法向量
    // 参数float scale=0.01 表示法向量长度缩放为0.01倍
    viewer.addPointCloudNormals<pcl::PointXYZ, pcl::Normal>(cloud, normals, 2, 0.01, "normals");
    //    viewer.addCoordinateSystem(1.0);

    while (!viewer.wasStopped()) {
        viewer.spinOnce();
    }
    return 0;
}

积分图估算点云法线

此方法进行法线估计只适用于有序点云,对于无序点云就只能采用其他方法。
代码实现：

#include <pcl/io/io.h>
#include <pcl/io/pcd_io.h>
#include <pcl/features/integral_image_normal.h>
#include <pcl/visualization/cloud_viewer.h>

int main()
{
    // load point cloud
    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>);
    pcl::io::loadPCDFile("./data/table_scene_mug_stereo_textured.pcd", *cloud);

    // estimate normals
    //待评估的法向量
    pcl::PointCloud<pcl::Normal>::Ptr normals(new pcl::PointCloud<pcl::Normal>);

    // 创建一个用于法线估计的对象并计算法线
    //采用积分图来评估法向量，必须是有序点云图
    pcl::IntegralImageNormalEstimation<pcl::PointXYZ, pcl::Normal> ne;

    // 有以下可选的估算方式：
    /**
    enum NormalEstimationMethod
    {
      COVARIANCE_MATRIX,
      AVERAGE_3D_GRADIENT,
      AVERAGE_DEPTH_CHANGE
    };
     COVARIANCE_MATRIX模式创建9个积分图像，以根据其局部邻域的协方差矩阵为特定点计算法线。
     AVERAGE_3D_GRADIENT模式创建6个积分图像以计算水平和垂直3D渐变的平滑版本，并使用这两个渐变之间的叉积计算法线。
     AVERAGE_DEPTH_CHANGE模式仅创建单个积分图像，并根据平均深度变化计算法线。
     */
    ne.setNormalEstimationMethod(ne.AVERAGE_3D_GRADIENT);
    ne.setMaxDepthChangeFactor(0.02f);
    ne.setNormalSmoothingSize(10.0f);
    ne.setInputCloud(cloud);
    ne.compute(*normals);

    // visualize normals
    pcl::visualization::PCLVisualizer viewer("PCL Viewer");
    viewer.setBackgroundColor(0.0, 0.0, 0.5);
    viewer.addPointCloudNormals<pcl::PointXYZ, pcl::Normal>(cloud, normals);

    while (!viewer.wasStopped()) {
        viewer.spinOnce();
    }
    return 0;
}