pcl里面的点特征直方图(PFH)

表面法线和曲率可以好的代表一个点的几何特征。然而它们算得很快，而且算法简单，但是它们不能捕获细节，它们只是点的近邻的几何特征的近似估计。作为一个直接的结论，大多数的场景往往会包括很多有着相似特征的点，这会减少它们所带来的消息量。

这次我们将引进一个叫做PFH(point feature histgrams)的3D特征描述器，同时还将显示出它的一些理论优势，并讨论它的一些实现细节。

PFH的目标是通过使用一个点周围的多维直方图的平均曲率来编码一个点的k个最近邻的几何属性。这个高维空间提供了一些有用的特征代表，同时关于6位位姿不变的，同时和可以很好的应对不同的采样密度和近邻的噪声水平。

一个点特征直方图是与点的最近邻和法线有关的，简单的说，它企图通过考虑各个预测法线方向之间的影响来捕获最好的采样平面。最终的高维空间因此和每个点的表面法线的方向有关。

下面的图展示了一个PFH在计算Pq这个点时候的影响区域图，Pq这个点用红色标记并放在半径为r的圆圈的中央，然后它的所有近邻(离点的距离小于半径r)都通过一个网格来间接的相邻着。下面的PFH描述器作为了一个直方图计算了所有匹配的点之间的关系，有一个O(k方)的复杂度。

计算pi和pj之间的绝对差，和与它们相关联的法线ns和nj，我们定义了一个固定的坐标系在某个点上。

我发现官网里面的上面的式子与下面的图对应不上，以图为标准，式子里面的那个v向量的方向和大小与图里面的都对不上。使用上面的uvw的坐标系，ns法线和ni法线的差可以化解为下面的3个角度的差异(因为他们的模都是一样的)：

d是两点距离的平方，是每对点的4要素，这样就把原来的12个值(2x每个点的{(x,y,z)，(nx,ny,nz)})减低为4个值。

使用PFH得到一个点的四要素可以用下面的方法:

computePairFeatures (const Eigen::Vector4f &p1, const Eigen::Vector4f &n1,
const Eigen::Vector4f &p2, const Eigen::Vector4f &n2,
float &f1, float &f2, float &f3, float &f4);

去生成最终的对每个查询点的PFH代表，四要素这个集合将放到直方图里面。分发的过程是把所有的特征(上面的3个角度，一个距离的平方)平均分割成b小份，然后数每个子区间上的个数。因为有3个特征是法线间的角度，所以他们可以很好的归一化并把它们放到三角函数的圆(高中知识)里面。一个分发的案例是把每个特征区间进行平均分配，然后我们就有了一个有b的四次方分发的柱状图。在这个图里面有一个区间的值不为0代表着，它这个特征有一定的值。下面的这幅图展示了一个点云里面不同点的PFH图。

有的情况下，d不是一个好的特征对于2.5D的数据集，比如机器人来说，距离会随着视点发生变化。所以，有的时候省略d是比较好的。

PFH特征在PCL里面是pcl_features库里面

默认的PFH用5个区间来分发，并不包括距离，因此有125比特的数组生成5的3次方，并把它们放在pcl::PFHSignature125这个点类型里面。

下面是一个代码段。

#include <pcl/point_types.h>
#include <pcl/features/pfh.h>
{
pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
pcl::PointCloud<pcl::Normal>::Ptr normals (new pcl::PointCloud<pcl::Normal> ());
... read, pass in or create a point cloud with normals ...
... (note: you can create a single PointCloud<PointNormal> if you want) ...
// Create the PFH estimation class, and pass the input dataset+normals to it
pcl::PFHEstimation<pcl::PointXYZ, pcl::Normal, pcl::PFHSignature125> pfh;
pfh.setInputCloud (cloud);
pfh.setInputNormals (normals);
// alternatively, if cloud is of tpe PointNormal, do pfh.setInputNormals (cloud);
// Create an empty kdtree representation, and pass it to the PFH 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> ());
//pcl::KdTreeFLANN<pcl::PointXYZ>::Ptr tree (new pcl::KdTreeFLANN<pcl::PointXYZ> ()); -- older call for PCL 1.5-
pfh.setSearchMethod (tree);
// Output datasets
pcl::PointCloud<pcl::PFHSignature125>::Ptr pfhs (new pcl::PointCloud<pcl::PFHSignature125> ());
// Use all neighbors in a sphere of radius 5cm
// IMPORTANT: the radius used here has to be larger than the radius used to estimate the surface normals!!!
pfh.setRadiusSearch (0.05);
// Compute the features
pfh.compute (*pfhs);
// pfhs->points.size () should have the same size as the input cloud->points.size ()*
}

PFHEstimation类的内部做了几件事情:

1.得到p的最近邻

2.给每个邻近的配对，计算3个角度

3.把所有的结果分发到直方图里面。

只计算一个单一的PFH特征代表:

computePointPFHSignature (const pcl::PointCloud<PointInT> &cloud,
const pcl::PointCloud<PointNT> &normals,
const std::vector<int> &indices,
int nr_split,
Eigen::VectorXf &pfh_histogram);

nr_spilit是被分发的某个区间，即代表了某个角度特征，pfh_histgram是结果的直方图。

因为有的时候法线的值可能是NaN(不是数字)或者无效值，所以可以用下面的代码进行优化，下面这段代码放在compute()这个函数的上面

for (int i = 0; i < normals->points.size(); i++)
{
if (!pcl::isFinite<pcl::Normal>(normals->points[i]))
{
PCL_WARN("normals[%d] is not finite\n", i);
}
}