1 引言

随着三维点云数据在地形测绘、文物修复、工程测量、医疗卫生等众多领域被广泛应用^[1],围绕其开展的研究成为了当下的热点。目前点云数据获取技术已较为成熟,与之相比,点云的数据处理技术还相对滞后。制约其发展的原因是:点云数据的海量性、高空间分辨率以及散乱性^[2]严重影响了其在实际工程应用中的可操作性。为了实现对点云数据的高效应用,关键是建立合理有效的空间索引方式^[3]。

目前,基于内存来构建空间索引是实现海量点云组织和管理的主要方法。依据分割方式,空间索引方式可分为规则分割法、对象分割法和组合索引技术三类^[4]。规则分割法是按照某种规则将空间分割成均匀的单元,并将空间中的实体对应分配到这些单元中。比较有代表性的规则分割方法包括规则格网方法、四叉树索引、八叉树索引、K维(KD)树索引、二叉空间分割树(BSP树)和R树^[5]等。该方法适用于点云的分布均匀且稀疏的情况,针对不均匀点云而言,很难找到合适的分割单元,从而导致空间索引效率低下;对象分割法利用一定形状的几何体以包围盒的方法对空间进行剖分。常用的包围体有:包围球、轴向包围盒、方向包围盒和凸包^[6]。但是目前没有一种包围盒可以在任意场景下达到最佳的分割效果。

无论是规则分割还是对象分割技术都存在局限性,为了克服上述弊端,引入了组合型的索引技术。通过将不同的索引结构有机结合起来,优势互补,达到很好的索引效果。赵江洪等^[6]提出了一种网格和KD树思想相结合的网格组合空间索引结构,该方法不仅解决了单一分辨率数据的冗余而且具有较好的查询效率;余飞等^[7]提出了一种基于四叉树和KD树的两级索引相互集成的数据组织管理方法应用于海量点云处理过程中,有效提高了点云数据可视化管理和处理的便捷程度;龚俊等^[8]提出了一种八叉树和R树集成的新三维索引方法,解决了高效迅速生成点云空间索引的难题;赵尔平等^[9]在八叉树编码的基础上,提出邻点差值渐进压缩和基于裁剪重叠区域进行冗余处理的R树空间索引方法。王永志等^[10]在上述方法的基础上进行了创新提出了一种八叉树与三维R^*树集成的组合索引结构——3DOR^*树,相较于文献[ 8]和文献[ 9]的方法具有高效的空间储存和查询的优势,并且可以应用于海量点云的处理。

鉴于前述算法中构建的索引结构对海量点云处理能力较差、管理效率低下的问题,本文提出了一种结合KD树空间切分思想的类八叉树混合点云索引结构,以此来实现对海量点云空间的组织。另外,为了对描述对象的数据进行有效管理,为散乱点云建立一个系统化的邻域结构,本文对传统的Morton编码方式进行了优化,并提出了一种新的k近邻搜索方式。从而在上述索引的基础上实现高效查询、检索等多种操作。实验结果表明,本文方法不仅使数据组织的整体结构趋于合理,也提高了海量点云数据的存取和检索效率,而且节省了存储空间,实现了对海量点云空间的合理管理

2 类八叉树索引结构

空间索引是依据空间对象的位置或按某种排列顺序来表征对象之间关系的一种数据结构,其中,KD树和线性八叉树被广泛应用于三维点云数据组织中^[11-12]。

KD树在涉及高维空间查找邻域等方面具有自身独特优势^[13],在无数据拓扑关系的情况下可以快速实现邻域查找,但其建立数据关系需要占据大量的内存空间。另外,对于海量点云进行建树,会由于深度过大而增加数据查找和回溯时间,从而影响搜索效率。

线性八叉树^[14]可以利用较小存储空间对海量数据进行高效操作与处理。但是当最小分割粒度过大时,后续的定位及查询效率会降低;粒度过小,八叉树的深度增加,同样也会导致效率下降^[15]。此外,相当数量、不存储任何数据的“空白节点”会影响树的平衡性。

综上所述,KD树在邻域搜索中具有优势,而八叉树可以利用较少的存储空间来对海量点云数据进行处理,但是这两种方法在点云存储和处理中都各有利弊。通过观察可知,它们能以合理的混合方式实现优势互补,考虑到组合索引方式的优越性^[16],本文提出了一种融合KD树和八叉树的组合索引结构。

2.1 索引构建思想

基于模型点云建立多层索引结构,在上层采用KD树的思想分割管理全局的点云数据,在下层利用八叉树组织和存储局部的点云数据,每个划分子空间的信息都保存在八叉树末端关联的叶结点中,并利用优化后的编码方式来构建索引结构。全局KD树和局部八叉树的多层索引结构如图1所示,其中虚线框内表示KD树叶结点与八叉树的对应关系。

图2为二维坐标下类八叉树分割空间示意图,其中红色点代表KD树的分割域(中值点),黄色虚线框内表示某KD树子空间,并对此子空间进行八叉树分割。

2.2 索引构建流程

建立类八叉树多层索引结构的具体算法步骤如算法1所示。

图 1. 类八叉树多层索引结构

Fig. 1. Octree-like multi-layer index structure

下载图片 查看所有图片

图 2. 类八叉树平面分割示意图

Fig. 2. Schematic of the octree-like plane division

下载图片 查看所有图片

Algorithm 1: Construct an octree-like index structure

Input: Point set P.

Output: Set of octree leaf nodes.

Initialization:KD tree segmentation threshold T_kd,Octree segmentation threshold T_oc,KD tree recursively decomposes the number of layers m.

Algorithm steps:

Step1:Build a KD tree. The original data point set P is used as the root node of the KD tree.

Step2:If the number of data points in the current KD tree node is larger than that in the KD tree segmentation parameter T_kd, the KD tree is decomposed in this space. Among them, T_kd=N/2^m,N represents the total number of point cloud data, m is the number of recursive decomposition layers of KD tree. According to the difference of original data, it is obtained by many tests. In this experiment, m = 5.

According to the cube{x_min, y_min, z_min, x_max, y_max, z_max} in this space, the partition domain is calculated, and the dimensions of the maximum coordinate variance in three dimensions of x, y, and z are calculated.

s_d=max sx,sy,sz

Where sx=1n∑i=1nxi-x-2,sy=1n∑i=1nyi-y-2,sz=1n∑i=1nzi-z-2, the median point p_m(x_m, y_m, z_m) on the dimension d is used as the segmentation point of the space, the over segmentation point p_m is parallel to the plane of the d axis, and the space is divided into two subspaces, the space cube is decomposed into two sub cubes, and the point set in the two subspaces is used as the two sub nodes, and the encircling box space is saved.

Step3:Each subspace is decomposed based on the KD tree, until the number of points in all spaces is less than equal to T_kd, and the n point set Q of the n leaf nodes of the KD tree and the bounding box space B for each point set, that is, the split point cloud data.

Step4:An octree is established for the split point cloud data. The following operation for each point set Q_i is as follows:

① The point set Q_i is used as the root node of octree, and its Morton code and bounding box space are calculated and preserved.

② If the number of data points in the octree node is larger than the octree partition parameter T_oc, the octree is decomposed into 8 sub cubes, and the point set in each sub cube is divided into 8 sub nodes, and the Morton code and packet enclosure space are calculated and saved. In this experiment, T_oc is defined as 0.3% of the original data points.

③ Gradually decomposing the KD tree recursively until all the points in the space are less than or equal to T_oc, save the leaf node of octree, that is, the segmented point cloud data.

Finally, n octree is obtained, and the segmented point cloud data are stored at each leaf node of each octree.

图3为构建类八叉树索引结构的流程图。

图 3. 建立类八叉树索引结构流程图

Fig. 3. Flowchart for establishing the octree-like index structure

下载图片 查看所有图片

3 基于优化编码的k近邻搜索方法

编码是为了对描述对象的数据进行有效的管理,为散乱点云建立一个系统化的邻域结构,以实现多种不同的操作,如查询、检索等^[17]。但是利用类八叉树思想切分空间会造成子空间边长各异且层数不一,另外传统KD树无法进行编码。基于此,本文对Morton编码方式进行了优化,提出了一种新的k近邻搜索方式来实现对点云数据的高效检索和管理。

3.1 Morton编码

Morton编码方式可以对点云数据进行有效编码,实现对数据的线性有序排列,同时保留空间位置关系。

针对空间中的某个正立方体而言,利用八叉树方法进行切分后,其结点位置对应的Morton编码M可表示为

M=mn-1…mk…m1m0,(1)

式中:m_k为八进制数,表示该节点的序号,k∈{0,1,…,n-1};m_k+1是m_k父节点的序号。这样,从m₀到m_n+1完整地表示了八叉树中每个叶子结点到根的路径,编码示意如图4所示。

图 4. 线性编码示意图

Fig. 4. Schematic of linear coding

下载图片 查看所有图片

3.2 基于优化编码的k近邻搜索方法

在类八叉树分割后的空间中,首先利用优化后的方法将数据点进行编码,然后利用编码搜索此点所在空间周围相邻的空间,并在搜索所得的相邻空间中选取k个点作为k近邻点。

搜索算法的具体思想如下:

1)对比点p与每个KD树的子空间包围盒的空间位置来判断其所处的子空间S_kd。点与KD空间的关系示意图如图5所示。

图 5. 点与KD空间的关系

Fig. 5. Relationship between point and KD space

下载图片 查看所有图片

2)对任意数据点p进行Morton编码。

若点所在S_kd空间包围盒为(x_min, y_min, z_min, x_max, y_max, z_max),其八叉树的最大分割层数为n,计算点p的索引值(i, j, k)为

i=floorx(xmax-xmin)/2nj=floory(ymax-ymin)/2nk=floorz(zmax-zmin)/2n,(2)

式中:floor()为向下取整操作。将索引值(i, j, k)转换为二进制表示方式,

i=i020+i121+…+il2l+…+in-12n-1j=j020+j121+…+jl2l+…+jn-12n-1k=k020+k121+…+kl2l+…+kn-12n-1,(3)

式中:i_l, j_l, k_l∈{0,1},1<l<n。

设点p对应的线性八叉树结点编码为M=m₁m₂…m_l(l≤n),则点p所在最大分割层数子空间的编码表示为

ml=4il+2jl+kl(1≤l≤n)。(4)

3)搜索点p相邻的八叉树空间,构成k近邻搜索的候选点集R。当搜索点的k近邻时,在数据点所在的子立方体及其周围相邻的26个子立方体中搜索,若p所在空间子立方体的索引值为(i, j, k),则其周围26个子立方体的索引值可由(i±δ, j±δ, k±δ)表示,δ∈{0,1}。保存这些相邻空间内的点,构成候选点集R。

4)计算点p与R中所有点的欧氏距离,取其中距离最短的k个点作为点p的k个近邻点:

di,j=(xi-xj)2+(yi-yj)2+(zi-zj)2。(5)

在上述编码过程中,利用类八叉树划分后的点云子空间搜索相邻空间,只需在局部进行k近邻搜索,大大提高了搜索速度。但空间分割停止的依据为是否低于阈值,因此并不是均匀分割。在查找某点相邻立方体空间时,其分割层次是未知的,所以不能有效地查找其周围的26个立方体空间。

图 6. 八叉树空间分割示意图。(a)均匀分割;(b)非均匀分割

Fig. 6. Spatial segmentation schematic of octree. (a) Uniform segmentation; (b) nonuniform segmentation

下载图片 查看所有图片

因此,假设点p对应的空间编码为M=m_n-1…m_k…m₁m₀,l(p)表示其编码长度,在点p的k近邻搜索过程中,本文算法确定搜索范围的策略如下。

对每个八叉树叶结点空间S_i:

1)若空间S_i的Morton码Q_s=Q_p,则对空间S_i进行搜索,将S_i加入候选空间集合R';

2)若空间S_i的编码长度l(S)<l(p)且sub[Q_p,l(p)]=Q_s,其中sub[Q_p,l(p)]表示Q_p的前l(p)位,则对空间S_i进行搜索,将S_i加入候选空间集合R'。

删去候选空间集合R'中重复的空间,得到空间集合R,即为点p所对应空间的相邻空间集合,取此空间集合内的所有点进行k近邻搜索,大大提高了搜索效率。

4 分析与讨论

为验证所提类八叉树混合索引结构有效性,本文从构造索引消耗时间、占用内存以及邻域搜索时间三个方面对KD树、八叉树、3DOR^*树、四叉+KD组合树以及本文类八叉树索引的基础性能进行测试;同时,通过下采样实验来展示上述方法对点云数据的管理效果;最后,通过实验探索了类八叉树索引结构的时间与空间效率与KD树分割阈值的关系。

4.1 实验数据及环境说明

本文采用了5组10⁵~10⁷数量级的点云数据进行实验。为了验证该方法的普适性,其中两组数量级较小,剩余三组均来源于Robotic 3D Scan Repository的海量点云数据。本文实验环境为Intel Core 4核i7-6700HQ 2.60 GHz,4 GB RAM内存,GeForce-GTX960M显卡,显存为4 GB。

4.2 构造索引速度与内存占比测试

索引构造时间对比如表1所示,其中,八叉树建树速度最快,KD树最慢,其耗时远大于其他方法,

而且该差异随着点云数量级的攀升而愈发明显。当点云数量超过400万时,KD树会由于内存溢出而无法建树。另外三种组合索引相比,本文方法建树时间略低于3DOR^*树方法,与最快的八叉树方法相差不多。

5种索引建树的内存占用对比如表2所示。本文方法占用内存最少,其他4种建树占用的内存较大,且差异随数量级的增加而变大。另外,对分布不均匀的点云而言,由于所构建的树深度过大而占用更多的内存,甚至导致堆栈溢出。

4.3 邻域搜索速度测试

索引结构的空间管理效率可以通过邻域搜索效率来进行表征,实验结果如表3所示。其中,组合索引在邻域搜索方面有着明显的优势,相比于单一八叉树搜索速度提高了1~2个数量级。另外,在点云数量较少且分布较为均匀的情况,KD树索引进行邻域搜索的时间最快,但对海量点云而言,仍会有堆栈溢出的困扰。将三种组合索引相比,本文方法比其他两种方法邻域搜索时间更加迅速,具有更好的检索性能。

4.4 点云采样效果对比

本实验对含88万点的Turbine Blade模型进行采样,采样为5万的精简点云。上述索引结构采样与随机采样效果如图7所示。其中,索引结构采样效果较随机采样更加规整自然;与单一索引结构相比,组合索引方式下采样所得点云更准确地保留了边界特征且分布均匀规律。而在三种组合索引方法中,本文方法处理之后的点云模型,纹理特征较其他两种方法得到了更好的保留。

图 7. Turbine Blade不同空间索引结构下采样效果局部对比。(a)随机采样;(b) KD树采样;(c)八叉树索引采样;(d)四叉+KD组合树采样;(e) 3DOR^*树采样;(f)类八叉树采样

Fig. 7. Local comparison of sampling effects under different spatial structures of Turbine Blade. (a) Random sampling; (b) KD-tree sampling; (c) octree index sampling; (d) quad+KD-tree sampling; (e) 3DOR^*-tree sampling; (f) octree-like sampling

下载图片 查看所有图片

另外,为了进一步验证对点云模型(外部轮廓和内部特征)的管理效果,对局部点云采样效果图进行比较分析,效果如图8所示。随机采样所得局部点集散乱且原始特征丢失最多;基于KD树索引的采样点分布较为均匀,但过于稀疏不能细致刻画特征;基于八叉树索引的采样效果与KD树相差无几但相对稠密;本文采样方法所得点云分布均匀且局部带状特征保留较好。而其他的组合索引方法采样得到的点云,如四叉+KD组合树采样方法局部点云较为稀疏,但是分布不太均匀;3DOR^*树采样方法处理后的局部点云仍较为稠密,且分布不够规整。

图 8. Turbine Blade不同空间索引结构下采样效果局部对比。(a)随机采样;(b) KD树采样;(c)八叉树索引采样;(d)四叉+KD组合树采样;(e) 3DOR^*树采样;(f)类八叉树采样

Fig. 8. Local comparison of sampling effects under different spatial structures of Turbine Blade. (a) Random sampling; (b) KD-tree sampling; (c) octree index sampling; (d) quad+KD-tree sampling; (e) 3DOR^*-tree sampling; (f) octree-like sampling

下载图片 查看所有图片

为了验证该方法的普适性和鲁棒性。本文采用另外54万的点云模型Haloxylon进行实验,该模型不同于上述较为规整的模型,其整体轮廓粗糙,局部特征多拟合为曲线。通过上述实验发现,组合索引明显优于单一索引结构,因此这里只利用随机方法以及三种组合索引方法进行了实验,采样后的整体与局部效果如图9、10所示。

从图9可以看出,就整体效果图而言,随机采样效果较差,其他三种采样结果相差无几。但是对图10中局部曲线特征进行分析可知,相比图10(a)~(c)的采样效果,图10(d)中类八叉树采样能够更好地保留模型中的曲线边界,而且就局部采样结果而言更加规整,点云分布更加均匀合理。

图 9. Haloxylon不同空间索引结构下采样效果局部对比。(a)随机采样;(b) KD树采样;(c)八叉树索引采样;(d)四叉+KD组合树采样;(e) 3DOR^*树采样;(f)类八叉树采样

Fig. 9. Local comparison of sampling effects under different spatial structures in Haloxylon. (a) Random sampling; (b) quad+KD-tree sampling; (c) 3DOR^*-tree sampling; (d) octree-like sampling

下载图片 查看所有图片

图 10. Haloxylon不同空间索引结构下采样效果局部对比。(a)随机采样;(b) KD树采样;(c)八叉树索引采样;(d)四叉+KD组合树采样;(e) 3DOR^*树采样;(f)类八叉树采样

Fig. 10. Local comparison of sampling effects under different spatial structures of Haloxylon. (a) Random sampling; (b) quad+KD-tree sampling; (c) 3DOR^*-tree sampling; (d) octree-like sampling

下载图片 查看所有图片

4.5 基于优化编码的k近邻搜索方法

本文方法对海量点云进行组织与管理的过程中,KD树叶子结点数对性能有影响。阈值过大会降低八叉树的搜索效率;反之,当阈值较小时增加了边界数据点的数量,会影响采样的准确度。在不同的KD阈值下,通过对构建索引结构的时间和内存占比进行对比,来分析能实现最佳搜索的分割参数。

由图11、图12可知,对数量级较小的点云而言,不同KD阈值对构建索引时间及空间占用的影响不大。对数量级较大的海量点云而言,当KD阈值较小时,建树所消耗的时间有明显的减少,但到达一定值后,减少趋势逐渐平缓;内存空间占用增加也有着类似的趋势。

因此,在建立类八叉树索引结构的过程中,时间和空间效率的提升可以转化为一个最优化问题,即寻找某一合适的KD阈值,使得建立索引结构的耗时较短且空间占用较小。如何寻求最优的KD阈值,需要针对不同点云数据的实际情况进行分析和实验。

图 11. 不同KD阈值下建立索引时间对比

Fig. 11. Comparison of indexing time under different KD thresholds

下载图片 查看所有图片

图 12. 不同KD阈值下建立索引占用内存对比

Fig. 12. Comparison of index memory usage under different KD thresholds

下载图片 查看所有图片

5 结论

本文提出了一种对海量点云进行快速组织与高效管理的类八叉树空间索引结构。该索引结构上层采用KD树的思想进行全局分割,在此基础上利用八叉树组织和存储分割后的点云数据,每个划分子空间的信息都保存在八叉树末端关联的叶结点中。另外,对传统八叉树编码方式进行了改进,使得在八叉树分割空间并不均匀的情况下,能够有效地对邻域进行搜索,大大提高搜索效率。该技术不仅有助于海量的点云数据存储和管理,也为基于点云数据进行的分析应用提供了技术支持。

本文使用几个数量级各异的数据集对类八叉树索引结构的搜索效率进行测试,从时间、内存、邻域搜索速度三个方面对KD树、八叉树和本文类八叉树索引进行对比分析,验证了类八叉树索引结构是一种有机嵌套的新型索引方式,综合了传统八叉树和KD树在时间和空间方面的优点,不仅能够提高传统KD树构造索引的速度和空间效率,而且能有效地改善传统八叉树对空间占用过大、邻域搜索耗时过长的问题。本文通过对海量点云进行基于类八叉树切分空间后的精简采样,显示出类八叉树索引结构对点云空间管理的合理性,通过对比实验可知,类八叉树索引方法可以实现点云空间的合理切分和管理。

但是本文提出的索引结构仍存在诸多不足:1)针对数量级更大、分布更加不均匀的点云数据集而言,在建立索引时对时间和内存的消耗依旧相当庞大,甚至会发生溢出等问题;2)该索引结构的效率受KD树阈值的影响较大,该阈值的选择需要通过最优问题的求解,但在本实验中仅依赖于经验和多次实验获得,不够智能化;3)该索引结构仅针对点云数据的空间位置进行操作,但对点云数据的曲率、法向量等信息无法进行存储和操作。