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

1. 背景

3DSC 是一種描述三維點云局部形狀的特征描述子，受二維 Shape Context 的啟發。它用于捕捉點云某一點局部的幾何分布信息，對點云配準、識別等任務非常有效。

2. 基本思想

3DSC 描述子通過統計某個關鍵點周圍點云在空間中分布的統計直方圖，編碼該關鍵點的局部結構信息。它的核心是將關鍵點附近的鄰域點投影到球坐標系統中，然后統計這些點在球坐標的不同區域內的分布。

3. 具體原理步驟

• 關鍵點選擇
  3DSC一般計算在關鍵點上（比如角點、邊緣點）【關鍵點提取算法會在后面章節進行講解說明】，提高效率和魯棒性。

• 鄰域定義
  對關鍵點以一定半徑定義鄰域（radius），找到鄰域內的所有點。

• 局部坐標系建立
  建立關鍵點局部參考坐標系，通常利用該點的法線方向定義z軸，通過構造局部坐標系保證描述子具有旋轉不變性。

• 球坐標系劃分
  在局部坐標系中，以關鍵點為球心，將鄰域空間劃分成若干個體素（bins），通常按半徑（r）、極角（θ）、方位角（φ）劃分為3個維度的離散格子，如下圖所示。
  

• 統計鄰域點分布
  計算鄰域中點在每個球體素內的點數，統計形成一個三維直方圖。

• 歸一化
  對直方圖進行歸一化，得到局部幾何結構的概率分布，作為特征描述子。

• 特征匹配
  計算兩個3DSC描述子之間的距離（例如χ2距離、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 ）特征描述》的學習，歡迎喜歡的朋友訂閱。