PCL点云库入门（第18讲）——PCL库点云特征之3DSC特征描述3D shape context descriptor

一、3DSC（3D Shape Context）特征算法原理

1. 背景

3DSC 是一种描述三维点云局部形状的特征描述子，受二维 Shape Context 的启发。它用于捕捉点云某一点局部的几何分布信息，对点云配准、识别等任务非常有效。

2. 基本思想

3DSC 描述子通过统计某个关键点周围点云在空间中分布的统计直方图，编码该关键点的局部结构信息。它的核心是将关键点附近的邻域点投影到球坐标系统中，然后统计这些点在球坐标的不同区域内的分布。

3. 具体原理步骤

• 关键点选择
  3DSC一般计算在关键点上（比如角点、边缘点）【关键点提取算法会在后面章节进行讲解说明】，提高效率和鲁棒性。

• 邻域定义
  对关键点以一定半径定义邻域（radius），找到邻域内的所有点。

• 局部坐标系建立
  建立关键点局部参考坐标系，通常利用该点的法线方向定义z轴，通过构造局部坐标系保证描述子具有旋转不变性。

• 球坐标系划分
  在局部坐标系中，以关键点为球心，将邻域空间划分成若干个体素（bins），通常按半径（r）、极角（θ）、方位角（φ）划分为3个维度的离散格子，如下图所示。
  

• 统计邻域点分布
  计算邻域中点在每个球体素内的点数，统计形成一个三维直方图。

• 归一化
  对直方图进行归一化，得到局部几何结构的概率分布，作为特征描述子。

• 特征匹配
  计算两个3DSC描述子之间的距离（例如χ²距离、Hellinger距离等）完成匹配。

4. 3DSC特点

• 保持对旋转和尺度的鲁棒性（通过局部坐标系和归一化处理）
• 能够有效捕捉局部三维结构信息
• 特征维度较高（通常上千维），计算和存储较大，适合关键点描述

---

二、主要成员函数和变量

1、成员变量详解

变量名	类型	说明
azimuth_bins_	std::size_t	方位角方向划分的 bin 数（默认 12）
elevation_bins_	std::size_t	仰角方向划分的 bin 数（默认 11）
radius_bins_	std::size_t	半径方向划分的 bin 数（默认 15）
min_radius_	double	最小半径，表示从点开始计算特征时的起始距离（默认 0.1）
point_density_radius_	double	用于计算点密度的搜索半径（默认 0.2）
radii_interval_	std::vector<float>	半径分桶区间，用于特征球体的划分
theta_divisions_	std::vector<float>	方位角分桶区间
phi_divisions_	std::vector<float>	仰角分桶区间
volume_lut_	std::vector<float>	每个 bin 对应的小体积 LUT（look-up table）
descriptor_length_	std::size_t	每个点的描述子维度 = azimuth_bins × elevation_bins × radius_bins
rng_	std::mt19937	随机数生成器，用于计算局部坐标系中 X 轴的选择（不确定性）
rng_dist_	std::uniform_real_distribution<float>	0~1 的均匀分布，用于上述随机采样

---

2、 核心成员函数讲解

构造函数

ShapeContext3DEstimation(bool random = false)

• 构造函数初始化默认参数和随机数生成器。
• 如果 random = true，则使用当前时间作为种子，否则使用固定值 12345。

---

initCompute()

bool initCompute () override;

• 初始化函数：

  • 计算并填充 radii_interval_, theta_divisions_, phi_divisions_
  • 计算每个体素的体积并填入 volume_lut_
  • 是描述子计算的前置步骤。

---

computePoint(...)

bool computePoint (std::size_t index, const pcl::PointCloud<PointNT>& normals, float rf[9], std::vector<float>& desc);

• 对单个点进行特征计算。

• 输入：

  • 点的索引
  • 法线集合
  • 点的局部参考坐标系 rf[9]

• 输出：

  • 描述子 desc（长度为 descriptor_length_）

---

computeFeature(...)

void computeFeature (PointCloudOut &output) override;

• 继承自 Feature<PointInT, PointOutT> 的虚函数。
• 用于计算整云的 3DSC 特征。
• 遍历所有索引点，调用 computePoint() 逐个生成描述子。

---

setMinimalRadius() / getMinimalRadius()

void setMinimalRadius(double radius);
double getMinimalRadius();

• 设置 / 获取最小半径 rmin（决定特征球从哪里开始采样）

---

setPointDensityRadius() / getPointDensityRadius()

void setPointDensityRadius(double radius);
double getPointDensityRadius();

• 设置 / 获取计算点密度时的搜索半径。

---

getAzimuthBins(), getElevationBins(), getRadiusBins()

• 获取当前 azimuth/elevation/radius 的 bin 数量。

---

rnd()

inline float rnd() { return rng_dist_(rng_); }

• 内部随机数生成函数（0~1 均匀分布），用于随机选择局部坐标轴方向。

---

3、 特征维度说明

描述子维度为：

descriptor_length_ = azimuth_bins_ × elevation_bins_ × radius_bins_;

默认情况下为：

12 × 11 × 15 = 1980维（对应 pcl::ShapeContext1980）

---

4、总结：整体流程（computeFeature）

1. 初始化（initCompute()） → 构建球形体素划分

2. 遍历输入点集（computeFeature()）

3. 对每个点：

   • 建立局部坐标系（通常由法线和随机 X 方向构建）
   • 邻域点转换到该局部坐标系
   • 根据方位角/仰角/半径落入 bin 中计数
   • 根据 bin 的体积归一化 → 得到直方图

4. 存入输出点云 PointCloudOut（即每个点一个1980维向量）

---

三、主要实现代码注解

1、计算单点 3D 形状上下文描述子（Shape Context 3D Descriptor）

// 模板函数：计算指定点的 Shape Context 描述子
template <typename PointInT, typename PointNT, typename PointOutT> 
bool pcl::ShapeContext3DEstimation<PointInT, PointNT, PointOutT>::computePoint (
    std::size_t index,                          // 当前计算的点在索引列表中的索引
    const pcl::PointCloud<PointNT> &normals,   // 所有法线
    float rf[9],                                // 输出：参考坐标系（Reference Frame，3x3 矩阵展开为数组）
    std::vector<float> &desc)                   // 输出：形状上下文描述子
{
  // rf 中的三个向量组成了局部参考系 RF：x_axis | y_axis | normal
  Eigen::Map<Eigen::Vector3f> x_axis (rf);
  Eigen::Map<Eigen::Vector3f> y_axis (rf + 3);
  Eigen::Map<Eigen::Vector3f> normal (rf + 6);

  // 在指定搜索半径内查找邻域点
  pcl::Indices nn_indices;
  std::vector<float> nn_dists;
  const std::size_t neighb_cnt = searchForNeighbors ((*indices_)[index], search_radius_, nn_indices, nn_dists);
  if (neighb_cnt == 0)
  {
    // 若无邻居，返回 NaN 描述子和空的参考系
    std::fill (desc.begin (), desc.end (), std::numeric_limits<float>::quiet_NaN ());
    std::fill (rf, rf + 9, 0.f);
    return (false);
  }

  // 找到最近邻点索引（用来获取该点的法线）
  const auto minDistanceIt = std::min_element(nn_dists.begin (), nn_dists.end ());
  const auto minIndex = nn_indices[std::distance (nn_dists.begin (), minDistanceIt)];

  // 获取当前点的位置
  Vector3fMapConst origin = (*input_)[(*indices_)[index]].getVector3fMap ();
  // 使用最近邻的法线作为当前点的法线
  normal = normals[minIndex].getNormalVector3fMap ();

  // 初始化 x_axis 为一个随机向量，用于之后与法线正交化
  x_axis[0] = rnd ();
  x_axis[1] = rnd ();
  x_axis[2] = rnd ();
  if (!pcl::utils::equal (normal[2], 0.0f))
    x_axis[2] = - (normal[0]*x_axis[0] + normal[1]*x_axis[1]) / normal[2];
  else if (!pcl::utils::equal (normal[1], 0.0f))
    x_axis[1] = - (normal[0]*x_axis[0] + normal[2]*x_axis[2]) / normal[1];
  else if (!pcl::utils::equal (normal[0], 0.0f))
    x_axis[0] = - (normal[1]*x_axis[1] + normal[2]*x_axis[2]) / normal[0];

  x_axis.normalize ();

  // 断言 x_axis 与 normal 正交
  assert (pcl::utils::equal (x_axis.dot(normal), 0.0f, 1E-6f));

  // y_axis = normal × x_axis，构成正交参考系
  y_axis.matrix () = normal.cross (x_axis);

  // 遍历邻域内每个邻居点
  for (std::size_t ne = 0; ne < neighb_cnt; ne++)
  {
    if (pcl::utils::equal (nn_dists[ne], 0.0f))
		  continue;

    // 获取邻居坐标
    Eigen::Vector3f neighbour = (*surface_)[nn_indices[ne]].getVector3fMap ();

    // --- 计算邻居点的极坐标 ---
    float r = std::sqrt (nn_dists[ne]); // 与中心点距离

    // 将点投影到切平面（normal 所在平面）
    Eigen::Vector3f proj;
    pcl::geometry::project (neighbour, origin, normal, proj);
    proj -= origin;
    proj.normalize ();

    // phi：投影在切平面后，与 x_axis 形成的角度（0 ~ 360°）
    Eigen::Vector3f cross = x_axis.cross (proj);
    float phi = pcl::rad2deg (std::atan2 (cross.norm (), x_axis.dot (proj)));
    phi = cross.dot (normal) < 0.f ? (360.0f - phi) : phi;

    // theta：当前邻居点与法线的夹角（0 ~ 180°）
    Eigen::Vector3f no = neighbour - origin;
    no.normalize ();
    float theta = normal.dot (no);
    theta = pcl::rad2deg (std::acos (std::min (1.0f, std::max (-1.0f, theta))));

    // --- 计算邻居所在的直方图 Bin(j,k,l) ---
    const auto rad_min = std::lower_bound(std::next (radii_interval_.cbegin ()), radii_interval_.cend (), r);
    const auto theta_min = std::lower_bound(std::next (theta_divisions_.cbegin ()), theta_divisions_.cend (), theta);
    const auto phi_min = std::lower_bound(std::next (phi_divisions_.cbegin ()), phi_divisions_.cend (), phi);

    const auto j = std::distance(radii_interval_.cbegin (), std::prev(rad_min));
    const auto k = std::distance(theta_divisions_.cbegin (), std::prev(theta_min));
    const auto l = std::distance(phi_divisions_.cbegin (), std::prev(phi_min));

    // --- 计算当前邻居点的局部点密度 ---
    pcl::Indices neighbour_indices;
    std::vector<float> neighbour_distances;
    int point_density = searchForNeighbors (*surface_, nn_indices[ne], point_density_radius_, neighbour_indices, neighbour_distances);
    if (point_density == 0)
      continue;

    // 权重 w = 体素归一化体积 LUT / 邻域点数
    float w = (1.0f / static_cast<float> (point_density)) *
              volume_lut_[(l*elevation_bins_*radius_bins_) +  (k*radius_bins_) + j];

    assert (w >= 0.0);
    if (w == std::numeric_limits<float>::infinity ())
      PCL_ERROR ("Shape Context Error INF!\n");
    if (std::isnan(w))
      PCL_ERROR ("Shape Context Error IND!\n");

    // 将权重累加到对应的直方图 bin 中
    desc[(l*elevation_bins_*radius_bins_) + (k*radius_bins_) + j] += w;

    assert (desc[(l*elevation_bins_*radius_bins_) + (k*radius_bins_) + j] >= 0);
  }

  // 注意：Shape Context 3D 的参考系不具备重复性，因此输出设为 0，提示用户
  std::fill (rf, rf + 9, 0);
  return (true);
}

此函数的关键步骤为：

1. 计算局部参考坐标系（使用最近邻的法线，并构造 x/y 轴）；
2. 搜索邻域点，并将邻居映射到形状上下文的三维极坐标空间；
3. 根据点在 (r, θ, φ) 空间中的位置，统计加权直方图；
4. 输出形状上下文描述子 desc（一维向量，长度为 radius_bins_ * elevation_bins_ * azimuth_bins_）；
5. 将 RF 置零表示其不可重复。

2、计算所有点 3D 形状上下文描述子（3DSC）

下面是你提供的 pcl::ShapeContext3DEstimation::computeFeature 函数的 逐行中文注释版，便于理解其作用和流程。

---

template <typename PointInT, typename PointNT, typename PointOutT> 
void pcl::ShapeContext3DEstimation<PointInT, PointNT, PointOutT>::computeFeature(PointCloudOut &output)
{
  // 确保描述子的长度为 1980（由半径/角度划分决定）
  assert (descriptor_length_ == 1980);

  // 假设输出点云初始是 dense（没有无效点）
  output.is_dense = true;

  // 遍历每一个需要计算描述子的点（由 indices_ 指定）
  for (std::size_t point_index = 0; point_index < indices_->size(); point_index++)
  {
    // 如果当前点不是有限（如有 NaN 或 Inf），填充 NaN 描述子并跳过
    if (!isFinite((*input_)[(*indices_)[point_index]]))
    {
      // 将 descriptor 数组设置为 NaN
      std::fill(output[point_index].descriptor, 
                output[point_index].descriptor + descriptor_length_,
                std::numeric_limits<float>::quiet_NaN());
      
      // 将局部参考帧 rf 设置为 0（表示无效）
      std::fill(output[point_index].rf, output[point_index].rf + 9, 0);

      // 标记整个输出点云为非 dense（包含无效点）
      output.is_dense = false;
      continue;
    }

    // 初始化当前点的描述子向量，长度为 descriptor_length_（1980）
    std::vector<float> descriptor(descriptor_length_);

    // 调用 computePoint 函数计算当前点的 3D Shape Context 描述子和局部参考帧 rf
    // 如果失败（如搜索不到邻域点），将该点视作无效点
    if (!computePoint(point_index, *normals_, output[point_index].rf, descriptor))
      output.is_dense = false;

    // 将计算出的 descriptor 拷贝到 output 中对应点的 descriptor 成员里
    std::copy(descriptor.begin(), descriptor.end(), output[point_index].descriptor);
  }
}

---

说明：

• descriptor_length_ == 1980 是由 Shape Context 直方图的划分决定的：
  通常为：radius_bins × elevation_bins × azimuth_bins，如 5 × 6 × 66 = 1980。

• rf 是局部参考帧（Local Reference Frame），由三个正交轴组成（x, y, normal），用于保证描述子的方向不变性。

• 如果点本身是无效的，或在 computePoint() 中无法找到有效邻域点，会将该点描述子设置为 NaN。

• output.is_dense 是一个标志位，表示结果是否全部为有效点（无 NaN）。

---

四、使用代码示例

/*****************************************************************//**
* \file   PCLFeature3DSCmain.cpp
* \brief
*
* \author YZS
* \date   Jun 2025
*********************************************************************/

#include <iostream>

#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/features/normal_3d.h>
#include <pcl/features/3dsc.h>
#include <pcl/search/kdtree.h>



int PCL3DSC(int argc, char** argv)
{
    // 加载点云
    std::string fileName = "E:/PCLlearnData/18/fragment.pcd";
    pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>());
    if (pcl::io::loadPCDFile<pcl::PointXYZ>(fileName, *cloud) == -1)
    {
        std::cerr << "无法读取点云文件"<< fileName << std::endl;
        return -1;
    }

    std::cout << "点云加载成功，点数: " << cloud->size() << std::endl;

    // 计算法线
    pcl::PointCloud<pcl::Normal>::Ptr normals(new pcl::PointCloud<pcl::Normal>());
    pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne;
    ne.setInputCloud(cloud);
    pcl::search::KdTree<pcl::PointXYZ>::Ptr tree(new pcl::search::KdTree<pcl::PointXYZ>());
    ne.setSearchMethod(tree);
    ne.setRadiusSearch(0.05); // 法线估计半径
    ne.compute(*normals);

    std::cout << "法线估计完成。" << std::endl;

    // 创建ShapeContext3D特征估计器
    pcl::ShapeContext3DEstimation<pcl::PointXYZ, pcl::Normal, pcl::ShapeContext1980> sc3d;
    sc3d.setInputCloud(cloud);  // 输入点云
    sc3d.setInputNormals(normals); // 输入法线
    sc3d.setSearchMethod(tree);
    sc3d.setRadiusSearch(0.2);  // 特征提取的球形邻域半径

    // 输出描述子
    pcl::PointCloud<pcl::ShapeContext1980>::Ptr descriptors(new pcl::PointCloud<pcl::ShapeContext1980>());
    sc3d.compute(*descriptors);

    std::cout << "3DSC特征计算完成，特征维度: " << pcl::ShapeContext1980::descriptorSize() << std::endl;
    std::cout << "输出特征个数: " << descriptors->size() << std::endl;

    // 打印前3个点的描述子片段
    for (size_t i = 0; i < std::min<size_t>(3, descriptors->size()); ++i)
    {
        std::cout << "点 " << i << " 描述子 (前50维): ";
        for (int j = 0; j < 50; ++j)
            std::cout << descriptors->points[i].descriptor[j] << " ";
        std::cout << std::endl;
    }

    return 0;
}

int main(int argc, char* argv[])
{
    PCL3DSC(argc, argv);
	std::cout << "Hello World!" << std::endl;
	std::system("pause");
	return 0;
}

结果展示：

---

至此完成第十八节PCL库点云特征之3DSC特征描述，下一节我们将进入《PCL库中点云特征之GASD（Globally Aligned Spatial Distribution ）特征描述》的学习，欢迎喜欢的朋友订阅。