PCL点云库入门（第22讲）——PCL库点云特征之RIFT 局部特征描述Rotation Invariant Feature Transform descriptors（RIFT）

一、算法原理

1、RIFT 特征简介

RIFT（轮换不变特征变换）最初由 Chao Dong 等人在 2009 年提出，用于多模态图像匹配。与 SIFT 的“方向归一化”思路不同，RIFT 通过 相对梯度方向直方图 建立旋转不变性，无需估计关键点主方向。随后该思想被推广到 3 D 点云（PCL 中的 RIFTEstimation）。下面以二维灰度图像为例讲述原理，再说明其扩展到 3 D 的要点。

2、原理讲述

1. 局部坐标系与旋转不变性

• 关键点中心 p：从检测器（MSER、DoG、Harris…）获得。

• 局部邻域半径 R：通常与尺度关联，确保包含足够上下文。

• 旋转不变策略：

  1. 在极坐标系内，只记录梯度方向与径向方向之间的夹角

     $${\theta }_q=\angle (\nabla I(q),q-p)$$

  2. 当整个补丁旋转时，径向方向随之旋转，而二者夹角保持不变 ⇒ 描述子天然轮换不变。

  3. 无需像 SIFT 那样估计主方向并旋转补丁，因此 对噪声与遮挡更稳健。

2. 描述子结构

每个环建立一个方向直方图 → 把所有直方图串联成最终向量，并可再 L2 归一化 提升鲁棒性。

3. 计算步骤（二维图像）

1. 梯度计算
   Sobel / Scharr 得到每像素梯度 $\nabla I=(I_x,I_y)$，存储幅值 $|\nabla I|$ 与方向 $\varphi$。

2. 环分区 & 方向量化

   • 将邻域以中心为圆心划分为 $d_r$ 个同面积环。

   • 对于环内每像素 q，求

     $${\theta }_q={\varphi }_q-\arctan 2(y_q-y_p,x_q-x_p)$$

     并映射到 $[0,2\pi )$ 后落入某方向 bin。

3. 加权投票

   • 幅值权重：投射到对应方向 bin。
   • 距离高斯权：离中心越近权重越大（可选）。
   • 双线性插值：减小量化误差。

4. 特征向量拼接 + 归一化

   • 拼接所有环的直方图 → 得到 RIFT 向量。
   • L1/L2 归一化可提升光照鲁棒性。

3、 扩展到三维点云

PCL 默认参数：nr_distance_bins = 4, nr_gradient_bins = 8 → 32 维向量。

1.3D RIFT 特征提取算法流程

步骤概览：

1. 计算强度梯度图（Gradient）
2. 确定每个关键点的邻域
3. 构建方向直方图
4. 构建多环形结构（radial bins）
5. 计算 RIFT 描述子

Step 1: 强度梯度计算（使用法线 + intensity）

• 使用 pcl::IntensityGradientEstimation 或用户自定义方式获取每个点的强度梯度 $\nabla I$。
• 梯度可以从反射强度、灰度或其他标量字段中计算。

得到每个点的梯度向量：

$$\nabla I=\left(\frac{\partial I}{\partial x},\frac{\partial I}{\partial y},\frac{\partial I}{\partial z}\right)$$

Step 2: 定义关键点邻域

• 对每个关键点 $p$，在其周围定义一个球形邻域半径 $r$。
• 所有在这个邻域内的点 $q\in \mathcal{N}(p)$ 都会参与特征描述子计算。

Step 3: 计算方向直方图（gradient orientation histogram）

对每个邻域内点，计算其相对位置向量 $\vec{v}=q-p$，然后计算梯度方向与向量方向之间的夹角。

梯度方向编码成若干个 角度 bin（常见如 8 个方向 bin）：

• 将方向空间（0 到 2π）划分成 $d_{\theta }$ 个角度 bin
• 每个点的梯度方向 ${\theta }_q$ 落入某个 bin，并根据梯度幅值对该 bin 累计投票

Step 4: 构建多环区域（radial bins）

• 将邻域按距离关键点的距离划分为若干个 环形区域（radial bins）

• 常用的划分为 $d_r$ 个半径分段，例如：

  • 第1环：0 ~ 1/3*r
  • 第2环：1/3r ~ 2/3r
  • 第3环：2/3*r ~ r

Step 5: 构造 RIFT 描述子

• 每个 radial bin 中统计一组 gradient orientation histogram（方向直方图）
• 整体 RIFT 描述子是一个维度为 $d_r\times d_{\theta }$ 的向量

例如，设定：

• nr_distance_bins = 4
• nr_gradient_bins = 8

则 RIFT 描述子长度为 $4\times 8=32$ 维

4、RIFT 描述子特点总结

RIFT 描述通常用于：

• 点云配准前的特征匹配
• Loop closure detection
• 场景识别等

二、主要成员变量和成员函数

1.类模板定义：

template <typename PointInT, typename GradientT, typename PointOutT>
class RIFTEstimation : public Feature<PointInT, PointOutT>

• PointInT：输入点类型（通常带有强度/灰度通道）
• GradientT：梯度类型（一般为 pcl::GradientXY 或用户自定义）
• PointOutT：输出特征点类型（如 pcl::Histogram<N>）

2.成员变量：

PointCloudGradientConstPtr gradient_;

• 输入梯度点云的常量指针。
• 每个点包含梯度方向信息（如强度梯度 vector）。

int nr_distance_bins_;

• RIFT 描述子的“半径方向”的分箱数量。
• 默认值为 4，表示描述子将邻域划分为 4 个同心圆环。

int nr_gradient_bins_;

• RIFT 描述子的“梯度方向”维度的分箱数量。
• 默认值为 8，表示将方向划分为 8 个扇形方向（45° 角度分辨率）。

3.成员函数：

1.公有成员函数：

RIFTEstimation () : gradient_ (), nr_distance_bins_ (4), nr_gradient_bins_ (8)

• 构造函数，设置默认距离和梯度分箱数。
• 设置特征名称为 "RIFTEstimation"。

2.输入设置相关：

inline void setInputGradient (const PointCloudGradientConstPtr &gradient)

• 设置输入的梯度图（一般由 pcl::IntensityGradientEstimation 预处理生成）。

inline PointCloudGradientConstPtr getInputGradient () const

• 返回输入的梯度图指针。

inline void setNrDistanceBins (int nr_distance_bins)
inline int getNrDistanceBins () const

• 设置 / 获取 RIFT 描述子的径向（环形）分区数。

inline void setNrGradientBins (int nr_gradient_bins)
inline int getNrGradientBins () const

• 设置 / 获取 RIFT 描述子的方向（角度）分区数。

3.特征计算相关：

void computeRIFT (const PointCloudIn &cloud, const PointCloudGradient &gradient,
                  int p_idx, float radius,
                  const std::vector<int> &indices,
                  const std::vector<float> &squared_distances,
                  Eigen::MatrixXf &rift_descriptor);

• 核心函数：计算某一点的 RIFT 特征描述子。

• 参数含义：

  • cloud：输入点云。
  • gradient：梯度图。
  • p_idx：查询点索引。
  • radius：邻域半径。
  • indices：查询点邻域点索引。
  • squared_distances：邻域内各点到中心点的距离平方。
  • rift_descriptor：输出的 RIFT 描述子矩阵（大小为 nr_distance_bins_ * nr_gradient_bins_）。

4.保护函数：

void computeFeature (PointCloudOut &output) override;

• 重载 Feature 接口中的 computeFeature() 函数。

• 该函数：

  • 遍历每个 interest point；
  • 查找邻域；
  • 调用 computeRIFT 填充输出特征；
  • 输出为 PointOutT 格式（通常是 pcl::Histogram<N>）。

三、关键部分代码注释

1.计算单点RIFT特征

template <typename PointInT, typename GradientT, typename PointOutT> 
void pcl::RIFTEstimation<PointInT, GradientT, PointOutT>::computeRIFT (
      const PointCloudIn &cloud, const PointCloudGradient &gradient, 
      int p_idx, float radius, const std::vector<int> &indices, 
      const std::vector<float> &sqr_distances, Eigen::MatrixXf &rift_descriptor)
{
  // 若邻域索引为空，则报错退出
  if (indices.empty ())
  {
    PCL_ERROR ("[pcl::RIFTEstimation] Null indices points passed!\n");
    return;
  }

  // 获取 RIFT 描述子的行数和列数，对应距离和梯度方向的 bin 数
  int nr_distance_bins = static_cast<int> (rift_descriptor.rows ());
  int nr_gradient_bins = static_cast<int> (rift_descriptor.cols ());

  // 获取中心点 p_idx 的 3D 坐标向量（Eigen::Vector3fMap 类型）
  pcl::Vector3fMapConst p0 = cloud.points[p_idx].getVector3fMap ();

  // 将 RIFT 描述子初始化为零矩阵
  rift_descriptor.setZero ();

  // 遍历每个邻域点
  for (std::size_t idx = 0; idx < indices.size (); ++idx)
  {
    // 获取当前点的 3D 坐标
    pcl::Vector3fMapConst point = cloud.points[indices[idx]].getVector3fMap ();

    // 获取梯度向量（转换为 Eigen::Vector3f）
    Eigen::Map<const Eigen::Vector3f> gradient_vector (& (gradient.points[indices[idx]].gradient[0]));

    // 计算梯度幅值
    float gradient_magnitude = gradient_vector.norm ();

    // 计算该点梯度与中心点连线方向之间的夹角（弧度）
    float gradient_angle_from_center = std::acos (
        gradient_vector.dot ((point - p0).normalized ()) / gradient_magnitude);
    
    // 检查角度值是否为有效浮点数
    if (!std::isfinite (gradient_angle_from_center))
      gradient_angle_from_center = 0.0;

    // 引入一个非常小的数防止除以零
    const float eps = std::numeric_limits<float>::epsilon ();

    // 将距离和角度归一化并映射到 bin 编号
    float d = static_cast<float> (nr_distance_bins) * std::sqrt (sqr_distances[idx]) / (radius + eps);
    float g = static_cast<float> (nr_gradient_bins) * gradient_angle_from_center / (static_cast<float> (M_PI) + eps);

    // 计算距离维度和角度维度上影响的 bin 范围，准备进行插值加权
    int d_idx_min = (std::max)(static_cast<int> (std::ceil (d - 1)), 0);
    int d_idx_max = (std::min)(static_cast<int> (std::floor (d + 1)), nr_distance_bins - 1);
    int g_idx_min = static_cast<int> (std::ceil (g - 1));
    int g_idx_max = static_cast<int> (std::floor (g + 1));

    // 遍历角度 bin 范围
    for (int g_idx = g_idx_min; g_idx <= g_idx_max; ++g_idx)  
    {
      // 因为角度是周期性的，所以需要对角度 bin 索引进行 wrap-around 处理
      int g_idx_wrapped = ((g_idx + nr_gradient_bins) % nr_gradient_bins); 

      // 遍历距离 bin 范围
      for (int d_idx = d_idx_min; d_idx <= d_idx_max; ++d_idx)
      {
        // 使用双线性插值（线性权重）更新对应 bin 值
        float w = (1.0f - std::abs (d - static_cast<float> (d_idx))) *
                  (1.0f - std::abs (g - static_cast<float> (g_idx)));

        // 将加权后的梯度幅值加入到对应的 RIFT 描述子 bin 中
        rift_descriptor (d_idx, g_idx_wrapped) += w * gradient_magnitude;
      }
    }
  }

  // 对描述子进行 L2 归一化（使得向量模长为 1），增强鲁棒性
  rift_descriptor.normalize ();
}

简要总结核心逻辑：

1. 输入数据：中心点、其邻域点索引、半径、邻域点到中心点的平方距离、以及对应的梯度。

2. 核心处理：

   • 根据每个邻域点的位置和梯度方向，计算其“距离 bin”与“方向 bin”。
   • 对每个有效 bin 进行插值加权累加。

3. 输出结果：

   • 得到一个二维的直方图 rift_descriptor，其每个元素为梯度幅值的加权和。
   • 最后归一化该描述子。

2.计算所有点RIFT特征

template <typename PointInT, typename GradientT, typename PointOutT> void
pcl::RIFTEstimation<PointInT, GradientT, PointOutT>::computeFeature (PointCloudOut &output)
{
  // 检查是否设置了搜索半径，如果为 0 则报错并清空输出
  if (search_radius_ == 0.0)
  {
    PCL_ERROR ("[pcl::%s::computeFeature] The search radius must be set before computing the feature!\n",
               getClassName ().c_str ());
    output.width = output.height = 0;
    output.points.clear ();
    return;
  }

  // 检查梯度方向的 bin 数是否有效（必须 > 0）
  if (nr_gradient_bins_ <= 0)
  {
    PCL_ERROR ("[pcl::%s::computeFeature] The number of gradient bins must be greater than zero!\n",
               getClassName ().c_str ());
    output.width = output.height = 0;
    output.points.clear ();
    return;
  }

  // 检查距离的 bin 数是否有效（必须 > 0）
  if (nr_distance_bins_ <= 0)
  {
    PCL_ERROR ("[pcl::%s::computeFeature] The number of distance bins must be greater than zero!\n",
               getClassName ().c_str ());
    output.width = output.height = 0;
    output.points.clear ();
    return;
  }

  // 检查输入是否设置了梯度图
  if (!gradient_)
  {
    PCL_ERROR ("[pcl::%s::computeFeature] No input gradient was given!\n", getClassName ().c_str ());
    output.width = output.height = 0;
    output.points.clear ();
    return;
  }

  // 检查梯度图中的点数与输入点云是否一致
  if (gradient_->points.size () != surface_->points.size ())
  {
    PCL_ERROR ("[pcl::%s::computeFeature] ", getClassName ().c_str ());
    PCL_ERROR ("The number of points in the input dataset differs from the number of points in the gradient!\n");
    output.width = output.height = 0;
    output.points.clear ();
    return;
  }

  // 创建一个用于存储 RIFT 描述子的矩阵（行：距离 bin，列：梯度 bin）
  Eigen::MatrixXf rift_descriptor (nr_distance_bins_, nr_gradient_bins_);
  
  // 存储邻域内的点索引
  std::vector<int> nn_indices;
  // 存储邻域内点与中心点之间的平方距离
  std::vector<float> nn_dist_sqr;

  // 遍历每个有效索引点（通常是采样下来的点）
  for (std::size_t idx = 0; idx < indices_->size (); ++idx)
  {
    // 在 search_radius_ 范围内查找当前点的邻居点
    tree_->radiusSearch ((*indices_)[idx], search_radius_, nn_indices, nn_dist_sqr);

    // 对当前点计算 RIFT 描述子
    computeRIFT (*surface_, *gradient_, (*indices_)[idx], static_cast<float> (search_radius_), nn_indices, nn_dist_sqr, rift_descriptor);

    // 将计算出的描述子内容复制到输出点云对应点的 histogram 数组中
    std::size_t bin = 0;
    for (Eigen::Index g_bin = 0; g_bin < rift_descriptor.cols (); ++g_bin)
      for (Eigen::Index d_bin = 0; d_bin < rift_descriptor.rows (); ++d_bin)
        output.points[idx].histogram[bin++] = rift_descriptor (d_bin, g_bin);
  }
}

总结

该函数的核心流程如下：

1. 参数检查：确保搜索半径、距离与梯度方向 bin 数有效，梯度图存在且维度匹配。
2. 邻域搜索：对每个采样点，使用 KdTree 进行半径邻域搜索。
3. 特征提取：调用 computeRIFT 函数计算当前点的 RIFT 描述子。
4. 写入输出：将描述子矩阵扁平化填入到输出点云的 histogram 数组中。

四、使用代码示例

/*****************************************************************//**
 * \file   PCLLearnRIFTFeathermain.cpp
 * \brief  
 * 
 * \details    
 * \author ZhengShi.Yang
 * \date   July 2025
 * \version   1.0
 * \dependencies  
 *********************************************************************/

#include <pcl/point_types.h>
#include <pcl/register_point_struct.h>
#include <pcl/io/pcd_io.h>
#include <pcl/common/common.h>
#include <pcl/features/normal_3d.h>
#include <pcl/features/rift.h>
#include <pcl/features/intensity_gradient.h>
#include <pcl/search/kdtree.h>
#include <pcl/kdtree/kdtree_flann.h>
#include <pcl/features/narf_descriptor.h>

#include <iostream>
#include <Eigen/Dense>

using namespace std;


int TestRIFT(int argc, char** argv)
{
	// 1. 读取点云数据
	std::string file = "E:/db_data/PCLLearnData/cloud_bin_1.pcd";
	// 加载输入点云（要求具有强度字段）
	pcl::PointCloud<pcl::PointXYZI>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZI>());
	if (pcl::io::loadPCDFile(file, *cloud) < 0)
	{
		PCL_ERROR("Cannot load PCD file!\n");
		return -1;
	}

	// 创建 KdTree 搜索结构（用于邻域查询）
	pcl::search::KdTree<pcl::PointXYZI>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZI>());
	// 计算强度梯度
	pcl::PointCloud<pcl::IntensityGradient>::Ptr gradients(new pcl::PointCloud<pcl::IntensityGradient>());
	pcl::IntensityGradientEstimation<pcl::PointXYZI, pcl::Normal, pcl::IntensityGradient> gradient_est;
	gradient_est.setInputCloud(cloud);

	//// 法线计算器（用于梯度估计）
	pcl::PointCloud<pcl::Normal>::Ptr normals(new pcl::PointCloud<pcl::Normal>());
	pcl::NormalEstimation<pcl::PointXYZI, pcl::Normal> normal_est;
	normal_est.setInputCloud(cloud);
	normal_est.setSearchMethod(tree);
	normal_est.setRadiusSearch(0.05); // 根据场景调整

	normal_est.compute(*normals);
	gradient_est.setInputNormals(normals);
	gradient_est.setSearchMethod(tree);
	gradient_est.setRadiusSearch(0.05);
	gradient_est.compute(*gradients);

	// 设置 RIFT 特征估计器
	pcl::RIFTEstimation<pcl::PointXYZI, pcl::IntensityGradient, pcl::Histogram<32>> rift_est;
	rift_est.setInputCloud(cloud);
	rift_est.setSearchSurface(cloud); // 可选：完整点云用于搜索
	rift_est.setInputGradient(gradients);
	rift_est.setSearchMethod(tree);
	rift_est.setRadiusSearch(0.05);  // 搜索半径
	rift_est.setNrDistanceBins(4);     // 空间划分距离 bins
	rift_est.setNrGradientBins(8);     // 梯度方向 bins

	// 输出特征点云
    pcl::PointCloud<pcl::Histogram<32>>::Ptr rift_descriptors(new pcl::PointCloud<pcl::Histogram<32>>());
	rift_est.compute(*rift_descriptors);

	std::cout << "RIFT descriptor computed for " << rift_descriptors->size() << " points." << std::endl;

	return 0;
}



int main(int argc,char *argv[])
{

	TestRIFT(argc, argv);
  std::cout<<"Hello PCL!"<<std::endl;
  std::system("pause");
  return 0;
}

注意事项

1. 点类型必须包含强度：RIFT 是基于强度梯度的描述子，因此必须是 PointXYZI。
2. 法线影响梯度质量：确保法线估计合理，radius 参数应匹配数据密度。
3. 匹配点对：RIFT 通常配合对应点对匹配方案使用，例如 FLANN 最近邻搜索。

结果展示：

至此完成第二十二节PCL库点云特征之RIFI特征描述，下一节我们将进入《PCL库中点云特征之NARF (Normal Aligned Radial Features)特征描述》的学习，欢迎喜欢的朋友订阅。