【PCL自学：ocTree】八叉树（octree）的原理及应用案例(点云压缩，搜索，空间变化)

一、什么是八叉树ocTree？

1.八叉树原理

  上世纪80年代，八叉树结构被提出来，用来表示空间中的区域划分，简单来说，空间可以被分为8个象限，想象一下假设空间中存在一个笛卡尔坐标系，则该坐标系将空间分为了8个象限，每个象限又可以按照这种方式再建立一个笛卡尔坐标系，再划分为8个象限，以此类推，空间中任何一块区域都可以被n个笛卡尔坐标划分。

  以上图为例：
  如果所生成的八叉树的节点可分为三类：

灰节点： 它对应的立方体部分的为V所占据；
   白节点： 它对应的立方体没有V的内容；
   黑节点： 它对应的立方体全部为V所占据。

后两类又称为叶节点。立方体的八叉树的逻辑结构是这样的：将立方体内的空间划分，并记录成树结构，其上的节点要么是叶节点，要么是有八个子节点的灰节点。
  因此任意空间都可以被划分为树状结构，而如果要找到某子空间，只需遍历八叉树即可找到，且效率较高。
  关于八叉树的存储结构，有兴趣深入研究的同学可以参考这篇文章：
  https://blog.csdn.net/qq_37855507/article/details/90957798

二、八叉树应用案例

1.点云压缩

  点云由巨大的数据集组成，这些数据集描述了与距离、颜色、法线等附加信息相关的三维点。此外，它们可以以很高的速度创建，因此会占用大量内存资源。一旦点云必须在速率受限的通信信道上存储或传输，就必须想办法对这类数据进行压缩。点云库提供点云压缩功能。它允许编码所有类型的点云，包括无组织的点云，其特征是不存在点索引、分辨率、密度和/或点顺序。此外，底层的八叉树数据结构能够有效地合并来自多个数据源的点云数据。

  下面，我们将解释如何有效地压缩单点云和点云流。在给出的示例中，我们使用OpenNIGrabber捕获点云，并使用PCL点云压缩技术进行压缩。请看如下示例代码：

#include <pcl/point_cloud.h>
#include <pcl/point_types.h>
#include <pcl/io/openni_grabber.h>
#include <pcl/visualization/cloud_viewer.h>

#include <pcl/compression/octree_pointcloud_compression.h>

#include <stdio.h>
#include <sstream>
#include <stdlib.h>

#ifdef WIN32
# define sleep(x) Sleep((x)*1000)
#endif

class SimpleOpenNIViewer
{
public:
 SimpleOpenNIViewer () :
    viewer (" Point Cloud Compression Example")
  {
  }
  
//在OpenNIGrabber捕获循环执行的回调函数中(以下部分代码)，首先将捕获的点云压缩到stringstream缓冲区中。接下来是一个解压步骤，该步骤将压缩的二进制数据解码为一个新的点云对象。然后将解码后的点云发送到点云查看器。
  void
  cloud_cb_ (const pcl::PointCloud<pcl::PointXYZRGBA>::ConstPtr &cloud)
  {
    if (!viewer.wasStopped ())
    {
      // 存储被压缩后点云的流
      std::stringstream compressedData;
      // 创建输出的解压缩后点云指针
      pcl::PointCloud<pcl::PointXYZRGBA>::Ptr cloudOut (new pcl::PointCloud<pcl::PointXYZRGBA> ());
      // 压缩点云
      PointCloudEncoder->encodePointCloud (cloud, compressedData);

      // 解压缩点云
      PointCloudDecoder->decodePointCloud (compressedData, cloudOut);


      // 显示解压缩后的点云
      viewer.showCloud (cloudOut);
    }
  }

  void
  run ()
  {
	// 布尔值：是否显示统计信息（是）
    bool showStatistics = true;

    // 配置文件列表查看这个头文件: /io/include/pcl/compression/compression_profiles.h
    // 此处是枚举类型，包括压缩配置(分辨率，在/离线，颜色信息)
    pcl::io::compression_Profiles_e compressionProfile = pcl::io::MED_RES_ONLINE_COMPRESSION_WITH_COLOR;
    // 实例化八叉树点云压缩用于编码和解码
    PointCloudEncoder = new pcl::io::OctreePointCloudCompression<pcl::PointXYZRGBA> (compressionProfile, showStatistics);
    PointCloudDecoder = new pcl::io::OctreePointCloudCompression<pcl::PointXYZRGBA> ();

    // 创建OpenNI设备的抓取器
    pcl::Grabber* interface = new pcl::OpenNIGrabber ();

    // 定义回调函数，用于回调点云采集设备
    std::function<void(const pcl::PointCloud<pcl::PointXYZRGBA>::ConstPtr&)> f =
      [this] (const pcl::PointCloud<pcl::PointXYZRGBA>::ConstPtr& cloud) { cloud_cb_ (cloud); };

    // 注册链接回调函数.此例子中是带有颜色信息的点云。
    boost::signals2::connection c = interface->registerCallback (f);

    // 开始接收点云，
    interface->start ();

    while (!viewer.wasStopped ())
    {
      sleep (1);
    }

    interface->stop ();

    // 删除点云压缩实例
    delete (PointCloudEncoder);
    delete (PointCloudDecoder);

  }

  pcl::visualization::CloudViewer viewer;

  pcl::io::OctreePointCloudCompression<pcl::PointXYZRGBA>* PointCloudEncoder;
  pcl::io::OctreePointCloudCompression<pcl::PointXYZRGBA>* PointCloudDecoder;

};

int
main ()
{
  SimpleOpenNIViewer v;
  v.run ();

  return (0);
}

  在run()函数中，我们创建了用于编码和解码的OctreePointCloudCompression类的实例。他们可以将压缩文件作为配置压缩算法的参数。所提供的压缩文件为openNI设备捕获的点云预定义了公共参数集。在本例中，我们使用MED_RES_ONLINE_COMPRESSION_WITH_COLOR它的意思是采用坐标编码精度为5立方毫米的轮廓线，并能进行颜色编码。它是进一步优化的快速在线压缩参数。压缩配置文件的完整列表，包括它们的配置，可以在文件/io/include/pcl/compression/compression_profiles.h找到。
这里列举一些参数以供参考：
第一个参数LOW,MED,HIGH代表的是压缩程度。
第二个ONLINE/OFFLINE代表是否在线压缩。
最后WITH.WITHOUT代表是否压缩颜色信息。

LOW_RES_ONLINE_COMPRESSION_WITHOUT_COLOR 	1立方厘米分辨率(八叉树体素体积)，无颜色，快速在线编码
LOW_RES_ONLINE_COMPRESSION_WITH_COLOR 		1立方厘米分辨率，彩色，快速在线编码
MED_RES_ONLINE_COMPRESSION_WITHOUT_COLOR 	5立方毫米分辨率，无色彩，快速在线编码
MED_RES_ONLINE_COMPRESSION_WITH_COLOR 	
HIGH_RES_ONLINE_COMPRESSION_WITHOUT_COLOR 	1立方毫米分辨率，无色彩，快速在线编码
HIGH_RES_ONLINE_COMPRESSION_WITH_COLOR 		以下以此类推
LOW_RES_OFFLINE_COMPRESSION_WITHOUT_COLOR 	
LOW_RES_OFFLINE_COMPRESSION_WITH_COLOR 	
MED_RES_OFFLINE_COMPRESSION_WITHOUT_COLOR 	
MED_RES_OFFLINE_COMPRESSION_WITH_COLOR 	
HIGH_RES_OFFLINE_COMPRESSION_WITHOUT_COLOR 	
HIGH_RES_OFFLINE_COMPRESSION_WITH_COLOR 	启用高级参数化的手动配置

  当然也可以选择手动更改压缩配置，为了能够完全访问所有与压缩相关的参数，OctreePointCloudCompression类的构造函数可以初始化附加的压缩参数。请注意，为了启用高级参数化，compressionProfile_arg参数需要设置为MANUAL_CONFIGURATION。

pcl::io::OctreePointCloudCompression< PointT, LeafT, BranchT, OctreeT >::OctreePointCloudCompression	(	

compression_Profiles_e 	compressionProfile_arg = MED_RES_ONLINE_COMPRESSION_WITH_COLOR,//定义压缩配置文件
bool 	showStatistics_arg = false,//获取压缩统计信息
const double 	pointResolution_arg = 0.001,//点坐标精度
const double 	octreeResolution_arg = 0.01,//最低八叉树级别的八叉树分辨率
bool 	doVoxelGridDownDownSampling_arg = false,//是否进行体素网格降采样
const unsigned int 	iFrameRate_arg = 30,//编码率
bool 	doColorEncoding_arg = true,//颜色编码
const unsigned char 	colorBitResolution_arg = 6 //颜色位分辨率
)		

  高级参数说明如下:
  compressionProfile_arg：启用高级参数时，需要将该参数设置为MANUAL_CONFIGURATION。showStatistics_arg:将压缩相关的统计数据打印到命令窗口。
  pointResolution_arg：定义点坐标的编码精度。该参数应设置为低于传感器噪声的值。
  octreeResolution_arg：该参数定义部署的八叉树的体素大小。较低的体素分辨率可以使压缩速度更快，但压缩性能下降。这在高帧/更新速率和压缩效率之间实现了平衡。
  doVoxelGridDownDownSampling_arg：如果激活，则只对层次八叉树数据结构进行编码。解码器在体素中心生成点。这样，点云在压缩时得到下采样，同时获得较高的压缩性能。
  iFrameRate_arg：点云压缩方案对点云进行区别编码。通过这种方式，传入点云和之前编码的点云之间的差异被编码，以实现最大的压缩性能。iFrameRate_arg允许指定流中进入点云不进行差异编码的帧速率(类似于视频编码中的I/ p帧)。
  doColorEncoding_arg:该选项启用颜色组件编码。
  colorBitResolution_arg:这个参数定义了每个颜色组件需要编码的位数。

2.用八叉树进行空间划分和搜索操作

  八叉树是一种基于树的数据结构，用于管理稀疏的三维数据。每个内部节点正好有8个子节点。在本节中，我们将学习如何在点云数据中使用八叉树进行空间划分和邻居搜索。特别是，我们解释了如何在体素搜索、K近邻搜索和半径搜索中执行邻点搜索。
  请看如下示例代码：

 #include <pcl/point_cloud.h>
 #include <pcl/octree/octree_search.h>
 
 #include <iostream>
 #include <vector>
 #include <ctime>
 
 int
 main ()
 {
   srand ((unsigned int) time (NULL));
	// [1]创建点云指针
   pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);
 
   // [2] 构造1000个随机的点云数据
   cloud->width = 1000;
   cloud->height = 1;
   cloud->points.resize (cloud->width * cloud->height);
 
   for (std::size_t i = 0; i < cloud->size (); ++i)
   {
     (*cloud)[i].x = 1024.0f * rand () / (RAND_MAX + 1.0f);
     (*cloud)[i].y = 1024.0f * rand () / (RAND_MAX + 1.0f);
     (*cloud)[i].z = 1024.0f * rand () / (RAND_MAX + 1.0f);
   }
 
   float resolution = 128.0f;
 	// [3]八叉树点云搜索实例
   pcl::octree::OctreePointCloudSearch<pcl::PointXYZ> octree (resolution);
 	// [4]将点云设置成八叉树结构
   octree.setInputCloud (cloud);
   octree.addPointsFromInputCloud ();
 	// [5]随机定义一个要查找的点
   pcl::PointXYZ searchPoint;
 
   searchPoint.x = 1024.0f * rand () / (RAND_MAX + 1.0f);
   searchPoint.y = 1024.0f * rand () / (RAND_MAX + 1.0f);
   searchPoint.z = 1024.0f * rand () / (RAND_MAX + 1.0f);
 
  //**************************体素搜索（搜索相同体素内的点）************************
 	// [6]创建点云索引容器
   std::vector<int> pointIdxVec;
	// [7]开始进行体素内近邻点搜索
  if (octree.voxelSearch (searchPoint, pointIdxVec))
  {
  // [8] 将查找点所在体素内的邻点全都输出打印到命令窗口
     std::cout << "Neighbors within voxel search at (" << searchPoint.x 
      << " " << searchPoint.y 
      << " " << searchPoint.z << ")" 
      << std::endl;
               
     for (std::size_t i = 0; i < pointIdxVec.size (); ++i)
    std::cout << "    " << (*cloud)[pointIdxVec[i]].x 
        << " " << (*cloud)[pointIdxVec[i]].y 
        << " " << (*cloud)[pointIdxVec[i]].z << std::endl;
   }
 
     //*************************K近邻搜索（搜索K个最近点）************************
   int K = 10;
 	// 创建点索引容器，用于按距离保存搜索到的点
   std::vector<int> pointIdxNKNSearch;
   	// 创建点距离容器 
   std::vector<float> pointNKNSquaredDistance;
 
   std::cout << "K nearest neighbor search at (" << searchPoint.x 
             << " " << searchPoint.y 
             << " " << searchPoint.z
             << ") with K=" << K << std::endl;
 	// 开始最近K邻搜索
   if (octree.nearestKSearch (searchPoint, K, pointIdxNKNSearch, pointNKNSquaredDistance) > 0)
   {
     for (std::size_t i = 0; i < pointIdxNKNSearch.size (); ++i)
       std::cout << "    "  <<   (*cloud)[ pointIdxNKNSearch[i] ].x 
                 << " " << (*cloud)[ pointIdxNKNSearch[i] ].y 
                 << " " << (*cloud)[ pointIdxNKNSearch[i] ].z 
                 << " (squared distance: " << pointNKNSquaredDistance[i] << ")" << std::endl;
   }
 
   //*************************K近邻搜索（搜索半径范围内的点）************************
 	// 创建半径范围内点索引的容器
   std::vector<int> pointIdxRadiusSearch;
   	// 创建搜索到的点的距离容器
   std::vector<float> pointRadiusSquaredDistance;
 
   float radius = 256.0f * rand () / (RAND_MAX + 1.0f);
 
   std::cout << "Neighbors within radius search at (" << searchPoint.x 
       << " " << searchPoint.y 
       << " " << searchPoint.z
       << ") with radius=" << radius << std::endl;
 
 // 开始半径搜索
   if (octree.radiusSearch (searchPoint, radius, pointIdxRadiusSearch, pointRadiusSquaredDistance) > 0)
   {
     for (std::size_t i = 0; i < pointIdxRadiusSearch.size (); ++i)
       std::cout << "    "  <<   (*cloud)[ pointIdxRadiusSearch[i] ].x 
                 << " " << (*cloud)[ pointIdxRadiusSearch[i] ].y 
                 << " " << (*cloud)[ pointIdxRadiusSearch[i] ].z 
                 << " (squared distance: " << pointRadiusSquaredDistance[i] << ")" << std::endl;
   }
 
}

  附加内容:
PCL八叉树组件提供了几种八叉树类型。它们的主要区别在于它们的单个叶节点特征。

OctreePointCloudPointVector(等于OctreePointCloud):这个八叉树可以在每个叶节点上保存一个点索引列表。
OctreePointCloudSinglePoint:这个八叉树类在每个叶节点上只有一个单点索引。仅存储分配给叶节点的最近的点索引。
OctreePointCloudOccupancy:这个八叉树在它的叶子节点上不存储任何点信息。它可以用于空间占用检查。
OctreePointCloudDensity:这个八叉树计算每个叶节点体素内的点的数量。它允许空间密度查询。
  如果需要高速创建八叉树，请查看八叉树双缓冲实现(Octree2BufBase类)。这个类同时在内存中保存两个并行的八叉树结构。除了搜索操作之外，这还支持空间变化检测。此外，高级内存管理减少了八叉树构建过程中的内存分配和回收操作。双缓冲八叉树实现可以通过模板参数OctreeT赋值给所有的OctreePointCloud类。
  所有八叉树都支持八叉树结构和八叉树数据内容的序列化和非序列化。

3.无序点云数据的空间变化检测

  八叉树是一种基于树的数据结构，用于组织稀疏的三维数据。在本节中，我们将学习如何使用八叉树实现来检测多个无序点云之间的空间变化，这些点云可以在大小、分辨率、密度和点顺序上变化。通过递归地比较八叉树的树结构，可以识别由体素配置差异表示的空间变化。此外，我们解释了如何使用pcl八叉树双缓冲技术，使我们能够在一段时间内有效地处理多个点云。

 #include <pcl/point_cloud.h>
 #include <pcl/octree/octree_pointcloud_changedetector.h>
 
 #include <iostream>
 #include <vector>
 #include <ctime>
 
 int
 main ()
{
  srand ((unsigned int) time (NULL));

  // 设置八叉树分辨率 即体素边长
  float resolution = 32.0f;

  // 实例化基于八叉树的点云变化检测类
  pcl::octree::OctreePointCloudChangeDetector<pcl::PointXYZ> octree (resolution);

	// 创建点云A的指针
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloudA (new pcl::PointCloud<pcl::PointXYZ> );

  // 为点云A生成一些随机点云数据
  cloudA->width = 128;
  cloudA->height = 1;
  cloudA->points.resize (cloudA->width * cloudA->height);
	
  for (std::size_t i = 0; i < cloudA->size (); ++i)
  {
    (*cloudA)[i].x = 64.0f * rand () / (RAND_MAX + 1.0f);
    (*cloudA)[i].y = 64.0f * rand () / (RAND_MAX + 1.0f);
    (*cloudA)[i].z = 64.0f * rand () / (RAND_MAX + 1.0f);
  }

  // 将点云A变成八叉树状结构
  octree.setInputCloud (cloudA);
  octree.addPointsFromInputCloud ();
  
// 点云A是我们的参考点云，用八叉树结构描述它的空间分布。
// OctreePointCloudChangeDetector类继承自Octree2BufBase类，
// 后者允许同时在内存中保存和管理两棵八叉树。
// 此外，它实现了一个内存池，可以重用已经分配的节点对象，因此在生成多个点云的八叉树时减少了昂贵的内存分配和回收操作。
//通过调用octree. switchbuffers()，我们重置了八叉树类，同时在内存中保留了之前的八叉树结构。

  // 八叉树缓冲区:重置八叉树，但在内存中保留以前的树结构。
  octree.switchBuffers ();
 // 创建点云B的指针
  pcl::PointCloud<pcl::PointXYZ>::Ptr cloudB (new pcl::PointCloud<pcl::PointXYZ> );
   
  // 为点云B生成一些随机点云数据
  cloudB->width = 128;
  cloudB->height = 1;
  cloudB->points.resize (cloudB->width * cloudB->height);

  for (std::size_t i = 0; i < cloudB->size (); ++i)
  {
    (*cloudB)[i].x = 64.0f * rand () / (RAND_MAX + 1.0f);
    (*cloudB)[i].y = 64.0f * rand () / (RAND_MAX + 1.0f);
    (*cloudB)[i].z = 64.0f * rand () / (RAND_MAX + 1.0f);
  }

  // 将点云B设置成八叉树结构
  octree.setInputCloud (cloudB);
  octree.addPointsFromInputCloud ();

// 为了检索存储在当前八叉树结构体素(基于clouddB)中的点，这些点在之前的八叉树结构(基于clouddA)中不存在，我们可以调用getPointIndicesFromNewVoxels方法，它返回结果点索引的向量。

	// 创建输出点云的索引容器
  std::vector<int> newPointIdxVector;

  // 从八叉树体素中获取点索引向量，这在之前的缓冲区中不存在
    octree.getPointIndicesFromNewVoxels (newPointIdxVector);

  // Output points
  std::cout << "Output from getPointIndicesFromNewVoxels:" << std::endl;
  for (std::size_t i = 0; i < newPointIdxVector.size (); ++i)
    std::cout << i << "# Index:" << newPointIdxVector[i]
              << "  Point:" << (*cloudB)[newPointIdxVector[i]].x << " "
              << (*cloudB)[newPointIdxVector[i]].y << " "
              << (*cloudB)[newPointIdxVector[i]].z << std::endl;

}

  由以上代码可知，点云A生成八叉树结构后，使用
pcl::octree::OctreePointCloudChangeDetectorpcl::PointXYZ octree (resolution);
和 octree.switchBuffers ();来分析无组织的点云的空间变化,将又一次生成的点云B放入之前点云A的树状结构中，并检索出点云B中没在之前点云A创建的树状结构中的点。
   通俗来说可以这样理解，假设A是一片人型点云，由人型点云创建了一个八叉树状结构，此时八叉树涵盖了A中所有点的位置，当第二片人形点云B出现，B和A的唯一区别就是B的左臂向上抬，此时利用上述方式，将B也放入之前A的树状结构，但是B的左臂上的点肯定不在之前A的树状结构中，因此会产生额外的树叶，这些额外的树叶内的点就是变化的点，即左臂的点。
例如下图中红色区域。

  有人会问，为什么黑色的影子部分的点明明消失了，为何不变红？因为点虽然消失了，但是左图中点创建的树状结构依然存在，因此不被判定为变化的点。
  需要注意的是，使用OcTree进行无组织的点的差分分析确实能起到较好结果，但是并不建议用在有序点的差分上，对于通常基于固定宽度和高度的单一2D深度/视差图像的有组织点云，直接对相应的2D深度图像进行差分分析可能会更快。

---

【博主简介】
  斯坦福的兔子，男，天津大学机械工程工学硕士。毕业至今从事光学三维成像及点云处理相关工作。因工作中使用的三维处理库为公司内部库，不具有普遍适用性，遂自学开源PCL库及其相关数学知识以备使用。谨此将自学过程与君共享。
博主才疏学浅，尚不具有指导能力，如有问题还请各位在评论处留言供大家共同讨论。
若前辈们有工作机会介绍欢迎私信。