[PCL] 点云特征估计

3D点云特征描述与提取是点云信息处理中最基础也是最关键的一部分，点云的识别、分割、重采样、配准、曲面重建等处理大部分算法，都严重依赖特征描述与提取的结果。从尺度上来分，一般分为局部特征的描述和全局特征的描述，例如局部的法线等几何形状特征的描述，全局的拓朴特征的描述，都属于3D点云特征描述与提取的范畴。

和数字图像处理中的特征点概念类似的，3D系统概念中所定义的点，使用它们的笛卡尔坐标x、y、z来表示。假设坐标系统的原点不随时间变化，可能有两个点p1和p2，分别在t1和t2处得到，具有相同的坐标。然而，比较这些点是一个ill-posed的问题，因为即使它们相对于某些距离度量(如欧几里得度量)是相等的，它们也可能是在完全不同的表面上采样的，因此，当它们与周围的其他点放在一起时，代表完全不同的信息。这是因为不能保证世界空间上在t1和t2之间没有变化。一些采集设备可能会为采样点提供额外的信息，如强度，甚至是颜色，但这并不能完全解决问题。

由于各种原因需要比较点的应用程序需要更好的特性和度量，以便能够区分几何表面。因此，三维点作为具有笛卡尔坐标的奇异实体的概念，取而代之的是一个新的概念，即局部描述符。

通过包含周围的邻居，可以在特征公式中推断表面几何形状。理想情况下，对于位于相同或相似表面上的点，合成的特征将非常相似(相对于某些度量)，而对于位于不同表面上的点，合成的特征将非常不同，如下图所示。一个好的点特征表示法与一个坏的点特征表示法不同，它能够在以下情况下捕获相同的局部表面特征，即下面几个条件通过能否获得相同的局部表面特征值判定点特征表示方式的优劣：

1. 刚体转换，即数据中的3D旋转和3D平移不应影响最终的特征向量F估计，即特征向量具有旋转平移不变性；
2. 改变采样密度，原则上，采样密度或多或少的局部表面贴片应该具有相同的特征向量特征；
3. 噪声，点特征表示必须在其特征向量中保留相同或非常相似的值，即使数据中存在轻微的噪声。

通常，PCL中特征向量利用快速kd-tree查询，使用近似法来计算查询点的最近邻元素，通常有两种查询类型：

• 确定查询点的k(用户给定参数)邻域(也称为K-搜索)；
• 确定半径为r的球体(也称为半径搜索)内查询点的所有邻域。

输入点云数据

pcl中的所有类都继承自pcl::PCLBase类，pcl::feature类通过两种方式接收输入数据：

整个点云数据集：setInputCloud (PointCloudConstPtr &)，试图估计给定输入点云中的每个点的特征。

点云数据集的子集：setInputCloud (PointCloudConstPtr &) 和 setIndices (IndicesConstPtr &)结合，尝试估计给定输入点云中在给定索引列表中有索引的每个点上的特征。默认情况下，如果没有给定一组索引，则会考虑云中的所有点。

此外，可以通过调用setSearchSurface(PointCloudConstPtr &)指定要使用的点邻域集。这个调用是可选的，当搜索平面没有给定时，默认使用输入点云数据集。

因为setInputCloud()接口总是必需的，所以我们可以使用最多四种组合来创建点云输入（setInputCloud(), setIndices(), setSearchSurface()）。

图示两个点云P和Q的特征点情况。

1. setIndices() = false, setSearchSurface() = false。因为我们不期望根据给定的一组索引或搜索表面来维护不同的实现副本，所以此时，PCL创建一组内部索引(std::vector<int>)，指向整个数据集(indexes=1..N，其中N为云中的点数)。在上图中，这对应于最左边的情况，首先，我们估计p_1的最近邻域，然后是p_2的最近邻域，以此类推，直到我们耗尽P中的所有点。

2. setIndices() = true, setSearchSurface() = false。在上图中，这对应于第二种情况。在这里，我们假设p_2的索引不是给定索引向量的一部分，所以在p2处没有邻居或特征。

3. setIndices() = false, setSearchSurface() = true。在上图中，这对应于第三种情况。如果Q={q_1, q_2}是另一个给定的点云作为输入，与P不同，P是Q的搜索面，那么q_1和q_2的邻居将从P中寻找和计算。

4. setIndices() = true, setSearchSurface() = true。在上图中，这对应于最后一种情况。邻域点将会在两个点云中进行计算，且考虑点的索引是否位于索引向量中。

当我们有一个非常密集的输入数据集，但我们不想估计所有点的特征，而是只希望估计一些关键点的特征，或希望处理下采样点云(例如使用pcl::VoxelGrid过滤器获取的下采样点云)，此时，setSearchSurface()应该使用。

例程

一旦确定邻域以后，查询点的邻域点可以用来估计一个局部特征描述子，它用查询点周围领域点描述采样面的几何特征，描述几何表面图形的一个重要属性，首先是推断它在坐标系中的方位，也就是估计他的法线，表面法线是表面的一个重要的属性，在许多领域都有重要的应用，如果用光源来生成符合视觉效果的渲染等。

下面的代码用于估计输入数据集中所有点的一组表面法线：normal_estimate.cpp

#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 ()
}

下面的代码片段将为输入数据集中的一个子集估计一组表面法线：

#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 a set of indices to be used. For simplicity, we're going to be using the first 10% of the points in cloud
  std::vector<int> indices (std::floor (cloud->points.size () / 10));
  for (std::size_t i = 0; i < indices.size (); ++i) indices[i] = i;

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

  // Pass the indices
  boost::shared_ptr<std::vector<int> > indicesptr (new std::vector<int> (indices));
  ne.setIndices (indicesptr);

  // 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 indicesptr->size ()
}

下面的代码片段将估计输入数据集中所有点的一组表面法线，但是将使用另一个数据集估计它们的最近邻。如前所述，当输入是曲面的下采样数据集时，这是一个很好的用途。

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

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

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

  ... create a downsampled version of it ...

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

  // Pass the original data (before downsampling) as the search surface
  ne.setSearchSurface (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 surface dataset.
  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_downsampled->points.size ()
}