论文:Persistent Point Feature Histograms for
3D Point Clouds
原理:
根据给定的特征点计算其临近点,在临近点范围内计算两点法向量之间的4个特征:
在pcl源码中,没有两点连线的模长,即没有f2特征,总共计算3个特征。
其他特征都与角度有关,如果3个特征都各自分成5份。那么一共可以分成5的3次方=125份,即由125个bins构成的直方图。
还有最重要的一点,两点要确定谁是source,谁是target,只有确定了source和target该4个特征才唯一确定,确定的原则是根据法向量与两点连线之间的角度:
上面确定好了source和target意味着任意两点只能计算唯一的4个特征,那么n个点总共可以计算
C
n
2
=
n
∗
(
n
−
1
)
2
C_{\mathrm{n}}^{2}=\frac{\mathrm{n} *(n-1)}{2}
Cn2=2n∗(n−1)个特征,也意味着复杂度为此值。
“pfh.h”文件
h文件一般是声明类,pfh是一种计算特征的类,因此继承自Feature类。
类首先会引用父类Feature中的成员变量,大致分为特征名字、待计算点的索引、kdtree搜索相关的变量、输入的点云集和法向量集。注意的是计算法向量的点云集要和surface_变量一致。
// 首先定义自己类的指针
using Ptr = shared_ptr<PFHEstimation<PointInT, PointNT, PointOutT> >;
using ConstPtr = shared_ptr<const PFHEstimation<PointInT, PointNT, PointOutT> >;
// 引用特征方法的名字及函数
using Feature<PointInT, PointOutT>::feature_name_;
using Feature<PointInT, PointOutT>::getClassName;
// 要计算的特征点的索引
using Feature<PointInT, PointOutT>::indices_;
// kdtre搜索
using Feature<PointInT, PointOutT>::k_;
using Feature<PointInT, PointOutT>::search_parameter_;
// 输入的点云、法向量集
using Feature<PointInT, PointOutT>::surface_;
using Feature<PointInT, PointOutT>::input_;
using FeatureFromNormals<PointInT, PointNT, PointOutT>::normals_;
// Feature<PointInT, PointOutT>::PointCloudOut 为 pcl::PointCloud<PointOutT>
using PointCloudOut = typename Feature<PointInT, PointOutT>::PointCloudOut;
using PointCloudIn = typename Feature<PointInT, PointOutT>::PointCloudIn;
该类中有一个变量用来存储计算的特征,使用map类型,key为参与计算特征的两点(ps, pt),value为计算的特征值,数据类型为Eigen::Vector4f(计算4个特征,但是用到了3个)。
std::map<std::pair<int, int>, Eigen::Vector4f, std::less<>, Eigen::aligned_allocator<std::pair<const std::pair<int, int>, Eigen::Vector4f> > > feature_map_;
上面的feature_map_变量可以存储计算的特征值,同时可以查询(ps, pt)是否已经计算过,减小计算量。
“pfh.cpp”文件
一般来说,cpp文件是用来实例化类模板的,但是改文件添加了计算4个特征的函数,是因为fpfh同样需要用到该计算特征的函数,只不过在fph中是链接的hpp文件,在fpfh中是使用pfh_tools.h来链接的。
下面的代码和论文中的公式一样,只不过4个特征值的对应顺序有差异。
///
// 计算pfh特征,传入两个点的坐标和法向量,输出的是4个特征值
bool
pcl::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)
{
Eigen::Vector4f dp2p1 = p2 - p1;
dp2p1[3] = 0.0f;
// f4是两相邻点的距离
f4 = dp2p1.norm ();
// 在进行数学运算前后,需要检查输出错误信息
if (f4 == 0.0f)
{
PCL_DEBUG ("[pcl::computePairFeatures] Euclidean distance between points is 0!\n");
f1 = f2 = f3 = f4 = 0.0f;
return (false);
}
// 计算f3,首先先确定谁是source,谁是target
Eigen::Vector4f n1_copy = n1,
n2_copy = n2;
n1_copy[3] = n2_copy[3] = 0.0f;
float angle1 = n1_copy.dot (dp2p1) / f4;
// Make sure the same point is selected as 1 and 2 for each pair
float angle2 = n2_copy.dot (dp2p1) / f4;
if (std::acos (std::fabs (angle1)) > std::acos (std::fabs (angle2)))
{
// switch p1 and p2
n1_copy = n2;
n2_copy = n1;
n1_copy[3] = n2_copy[3] = 0.0f;
dp2p1 *= (-1);
f3 = -angle2;
}
else
f3 = angle1;
// Create a Darboux frame coordinate system u-v-w
// u = n1; v = (p_idx - q_idx) x u / || (p_idx - q_idx) x u ||; w = u x v
// 齐次坐标相乘,第四维不需要管,只需要计算完叉乘后赋值0即可,因为叉乘是计算与现有向量都垂直的向量,因此前三项肯定满足
Eigen::Vector4f v = dp2p1.cross3 (n1_copy);
v[3] = 0.0f;
float v_norm = v.norm ();
// 计算完叉乘后,检查结果模长
if (v_norm == 0.0f)
{
PCL_DEBUG ("[pcl::computePairFeatures] Norm of Delta x U is 0!\n");
f1 = f2 = f3 = f4 = 0.0f;
return (false);
}
// Normalize v
v /= v_norm;
Eigen::Vector4f w = n1_copy.cross3 (v);
// Do not have to normalize w - it is a unit vector by construction
v[3] = 0.0f;
f2 = v.dot (n2_copy);
w[3] = 0.0f;
// Compute f1 = arctan (w * n2, u * n2) i.e. angle of n2 in the x=u, y=w coordinate system
f1 = std::atan2 (w.dot (n2_copy), n1_copy.dot (n2_copy)); // @todo optimize this
return (true);
}
“pfh.hpp”文件
hpp文件是代码的主体部分,主要包括构造函数、计算特征直方图函数computePairFeatures、computeFeature重载函数。
//
template <typename PointInT, typename PointNT, typename PointOutT> void
pcl::PFHEstimation<PointInT, PointNT, PointOutT>::computePointPFHSignature (
const pcl::PointCloud<PointInT> &cloud, const pcl::PointCloud<PointNT> &normals,
const std::vector<int> &indices, int nr_split, Eigen::VectorXf &pfh_histogram)
{
int h_index, h_p;
// Clear the resultant point histogram
pfh_histogram.setZero ();
// 一个统计点所占的百分比 一共有C(n, 2)排列组合个点,因为要从n个临近点中随机选取两个点来计算特征
// Factorization constant
float hist_incr = 100.0f / static_cast<float> (indices.size () * (indices.size () - 1) / 2);
std::pair<int, int> key;
bool key_found = false;
// Iterate over all the points in the neighborhood
for (std::size_t i_idx = 0; i_idx < indices.size (); ++i_idx)
{
for (std::size_t j_idx = 0; j_idx < i_idx; ++j_idx)
{
// If the 3D points are invalid, don't bother estimating, just continue
// 有限数
if (!isFinite (cloud.points[indices[i_idx]]) || !isFinite (cloud.points[indices[j_idx]]))
continue;
// 表示将查询过的点对存入缓存,在每次计算点对特征时,首先查询缓存中是否已经存在
if (use_cache_)
{
// In order to create the key, always use the smaller index as the first key pair member
int p1, p2;
// if (indices[i_idx] >= indices[j_idx])
// {
p1 = indices[i_idx];
p2 = indices[j_idx];
// }
// else
// {
// p1 = indices[j_idx];
// p2 = indices[i_idx];
// }
key = std::pair<int, int> (p1, p2);
// Check to see if we already estimated this pair in the global hashmap
// Eigen管理内存和C++11中的方法是不一样的,所以需要单独强调元素的内存分配和管理
// 使用Eigen创建stl容器等时,需要写明内存分配器
// pfh_tuple_是存放四个feature的数组Eigen::Vector4f
std::map<std::pair<int, int>, Eigen::Vector4f, std::less<>, Eigen::aligned_allocator<std::pair<const std::pair<int, int>, Eigen::Vector4f> > >::iterator fm_it = feature_map_.find (key);
// 先搜索是否已经计算过(ps, pt)该组的特征
if (fm_it != feature_map_.end ())
{
pfh_tuple_ = fm_it->second;
key_found = true;
}
else
{
// Compute the pair NNi to NNj
if (!computePairFeatures (cloud, normals, indices[i_idx], indices[j_idx],
pfh_tuple_[0], pfh_tuple_[1], pfh_tuple_[2], pfh_tuple_[3]))
continue;
key_found = false;
}
}
else
if (!computePairFeatures (cloud, normals, indices[i_idx], indices[j_idx],
pfh_tuple_[0], pfh_tuple_[1], pfh_tuple_[2], pfh_tuple_[3]))
continue;
// Normalize the f1, f2, f3 features and push them in the histogram
// std::floor向下取整数
// nr_split每个角度特征分类数目,分为两类,共有16
/**
* f1 = arctan (w * n2, u * n2)
* f2 = v · nt
* f3 = u · (pt− ps) / || pt− ps ||
* f4 = || pt− ps ||
*/
// d_pi_ = 1.0 / (2.0 * M_PI)
// 注意:这里在计算属于直方图的哪个bin时,用的是3个特征,没有使用模长特征
// 因此,后面在计算直方图时,bin的个数为n的3次方
f_index_[0] = static_cast<int> (std::floor (nr_split * ((pfh_tuple_[0] + M_PI) * d_pi_)));
if (f_index_[0] < 0) f_index_[0] = 0;
if (f_index_[0] >= nr_split) f_index_[0] = nr_split - 1;
f_index_[1] = static_cast<int> (std::floor (nr_split * ((pfh_tuple_[1] + 1.0) * 0.5)));
if (f_index_[1] < 0) f_index_[1] = 0;
if (f_index_[1] >= nr_split) f_index_[1] = nr_split - 1;
f_index_[2] = static_cast<int> (std::floor (nr_split * ((pfh_tuple_[2] + 1.0) * 0.5)));
if (f_index_[2] < 0) f_index_[2] = 0;
if (f_index_[2] >= nr_split) f_index_[2] = nr_split - 1;
// Copy into the histogram
h_index = 0;
h_p = 1;
for (const int &d : f_index_)
{
h_index += h_p * d;
h_p *= nr_split;
}
// 在h_index处直方图上添加一个点的百分比
pfh_histogram[h_index] += hist_incr;
if (use_cache_ && !key_found)
{
// Save the value in the hashmap
feature_map_[key] = pfh_tuple_;
// Use a maximum cache so that we don't go overboard on RAM usage
key_list_.push (key);
// Check to see if we need to remove an element due to exceeding max_size
if (key_list_.size () > max_cache_size_)
{
// Remove the oldest element.
feature_map_.erase (key_list_.front ());
key_list_.pop ();
}
}
}
}
}
/**调用的计算FPH特征的函数
* @输入:存储特征的点云对象,使用pcl::PointCloud<PointOutT>定义,PointOutT为直方图数据类型
* @返回:无
* @流程:
* 1. 根据输入的特征点索引indices_计算临近点
* 2. 根据临近点计算特征,计算直方图
*/
//
template <typename PointInT, typename PointNT, typename PointOutT> void
pcl::PFHEstimation<PointInT, PointNT, PointOutT>::computeFeature (PointCloudOut &output)
{
// Clear the feature map
// key是(ps, pt),value是4个特征的值向量
feature_map_.clear ();
std::queue<std::pair<int, int> > empty;
// key_list_存放的是已经计算过特征的(ps, pt),这里的作用是将key_list_清空
std::swap (key_list_, empty);
pfh_histogram_.setZero (nr_subdiv_ * nr_subdiv_ * nr_subdiv_);
// Allocate enough space to hold the results
// \note This resize is irrelevant for a radiusSearch ().
std::vector<int> nn_indices (k_);
std::vector<float> nn_dists (k_);
output.is_dense = true;
// Save a few cycles by not checking every point for NaN/Inf values if the cloud is set to dense
// 输入的点云是密集dense的,不包含NaN/Inf值
if (input_->is_dense)
{
// Iterating over the entire index vector
for (std::size_t idx = 0; idx < indices_->size (); ++idx)
{
if (this->searchForNeighbors ((*indices_)[idx], search_parameter_, nn_indices, nn_dists) == 0)
{
for (Eigen::Index d = 0; d < pfh_histogram_.size (); ++d)
output.points[idx].histogram[d] = std::numeric_limits<float>::quiet_NaN ();
output.is_dense = false;
continue;
}
// Estimate the PFH signature at each patch
computePointPFHSignature (*surface_, *normals_, nn_indices, nr_subdiv_, pfh_histogram_);
// Copy into the resultant cloud
for (Eigen::Index d = 0; d < pfh_histogram_.size (); ++d)
output.points[idx].histogram[d] = pfh_histogram_[d];
}
}
// 如果输入的点云不是稠密点云,意味着点云集中可能存在NaN值,首先调用!isFinite ((*input_)[(*indices_)[idx]])语句检查,剩余的语句和稠密的一样
else
{
// Iterating over the entire index vector
for (std::size_t idx = 0; idx < indices_->size (); ++idx)
{
// 根据输入的特征点索引indiecs_,计算临近点nn_indices索引数组
if (!isFinite ((*input_)[(*indices_)[idx]]) ||
this->searchForNeighbors ((*indices_)[idx], search_parameter_, nn_indices, nn_dists) == 0)
{
// 如果找不到临近点,那么给直方图数组histogram赋值NaN
for (Eigen::Index d = 0; d < pfh_histogram_.size (); ++d)
output.points[idx].histogram[d] = std::numeric_limits<float>::quiet_NaN ();
output.is_dense = false;
continue;
}
// Estimate the PFH signature at each patch
// 计算特征直方图
computePointPFHSignature (*surface_, *normals_, nn_indices, nr_subdiv_, pfh_histogram_);
// Copy into the resultant cloud
// 这里的pfh_histogram_的大小与参数nr_split有关,如果nr_split=2,那么大小为2的3次方=8,注意这里是计算的3个特征,而不是论文中计算的4个特征,少了一个模长特征
// points定义:std::vector<PointT, Eigen::aligned_allocator<PointT> > points;
// histogram是PointT中的成员
for (Eigen::Index d = 0; d < pfh_histogram_.size (); ++d)
output.points[idx].histogram[d] = pfh_histogram_[d];
}
}
}