简述

点云法向量估计这个问题，相信很多人在点云处理，曲面重建的过程中遇到过。表面法线是几何体面的重要属性。而点云数据集在真实物体的表面表现为一组定点样本。对点云数据集的每个点的法线估计，可以看作是对表面法线的近似推断。在开源库提供我们调用便利的同时，了解其实现原理也有利于我们对问题的深刻认识！格物要致知：）

---

原理

确定表面一点法线的问题近似于估计表面的一个相切面法线的问题，因此转换过来以后就变成一个最小二乘法平面拟合估计问题。

• 平面方程
  

  为平面上点处法向量的方向余弦，为原点到平面的距离。
  
  求平面方程即转化为求四个参数。

• 求解过程

1. 待拟合平面点集
待拟合的平面方程：
任意点到平面的距离：

2. 要获得最佳拟合平面，则需要满足：

因此，转化为求解极值的问题，

3. 分别对求偏导


将带入任意点到平面的距离公式：

继续求偏导

令
则：

同理：

将上述三式统一：

易得：

即转化到了求解矩阵A的特征值与特征向量的问题，矩阵即为n个点的协方差矩阵。即为该矩阵的一个特征向量。

4. 求最小特征向量
如上所示，求得的特征向量可能不止一个，那么如何来选取特征向量，使得求得法向量为最佳拟合平面的法向量呢？
由(内积形式)，
,


由
因此，最小特征值对应的特征向量即为法向量

---

程序应用

• PCL中的NormalEstimation

    #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 ...

  // Create the normal estimation class, and pass the input dataset to it
  pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne;
  ne.setInputCloud (cloud);

  // Create an empty kdtree representation, and pass it to the normal estimation object.
  // Its content will be filled inside the object, based on the given input dataset (as no other search surface is given).
  pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ> ());
  ne.setSearchMethod (tree);

  // Output datasets
  pcl::PointCloud<pcl::Normal>::Ptr cloud_normals (new pcl::PointCloud<pcl::Normal>);

  // Use all neighbors in a sphere of radius 3cm
  ne.setRadiusSearch (0.03);

  // Compute the features
  ne.compute (*cloud_normals);

  // cloud_normals->points.size () should have the same size as the input cloud->points.size ()*
}

• OpenMP加速法线估计
  PCL提供了表面法线估计的加速实现，基于OpenMP使用多核/多线程来加速计算。 该类的名称是pcl :: NormalEstimationOMP，其API与单线程pcl :: NormalEstimation 100％兼容。 在具有8个内核的系统上，一般计算时间可以加快6-8倍。

include <pcl/point_types.h>
#include <pcl/features/normal_3d_omp.h>

{
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);

  ... read, pass in or create a point cloud ...

  // Create the normal estimation class, and pass the input dataset to it
  pcl::NormalEstimationOMP<pcl::PointXYZ, pcl::Normal> ne;
  ne.setNumberOfThreads(12);  // 手动设置线程数，否则提示错误
  ne.setInputCloud (cloud);

  // Create an empty kdtree representation, and pass it to the normal estimation object.
  // Its content will be filled inside the object, based on the given input dataset (as no other search surface is given).
  pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ> ());
  ne.setSearchMethod (tree);

  // Output datasets
  pcl::PointCloud<pcl::Normal>::Ptr cloud_normals (new pcl::PointCloud<pcl::Normal>);

  // Use all neighbors in a sphere of radius 3cm
  ne.setRadiusSearch (0.03);

  // Compute the features
  ne.compute (*cloud_normals);

  // cloud_normals->points.size () should have the same size as the input cloud->points.size ()*
}