一、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 | normalEigen::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 数组设置为 NaNstd::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 ）特征描述》的学习，欢迎喜欢的朋友订阅。