计算机应用   2016, Vol. 36 Issue (10): 2863-2869  DOI: 10.11772/j.issn.1001-9081.2016.10.2863
0

引用本文 

赵京东, 杨凤华. 激光散乱点云K最近邻搜索算法[J]. 计算机应用, 2016, 36(10): 2863-2869.DOI: 10.11772/j.issn.1001-9081.2016.10.2863.
ZHAO Jingdong, YANG Fenghua. K-nearest neighbor searching algorithm for laser scattered point cloud[J]. JOURNAL OF COMPUTER APPLICATIONS, 2016, 36(10): 2863-2869. DOI: 10.11772/j.issn.1001-9081.2016.10.2863.

基金项目

国家自然科学基金资助项目(61104136);曲阜师范大学实验室开放基金资助项目(sk201418)

通信作者

赵京东(1962—),男,山东莱州人,教授,主要研究方向:CAD、数字图像处理,E-mail:zhaojd@mail.qfnu.edu.cn

作者简介

杨凤华(1962—),女,山东新泰人,副教授,硕士,主要研究方向:非线性泛函分析

文章历史

收稿日期:2016-03-23
修回日期:2016-05-03
激光散乱点云K最近邻搜索算法
赵京东, 杨凤华    
曲阜师范大学 数学科学学院, 山东 曲阜 273165
摘要: 针对激光散乱点云的数据量大,且具有面型的特点,为降低存储器使用量,提高散乱点云的处理效率,提出了一种散乱点云K最近邻(KNN)搜索算法。首先,利用多级分块、动态链表的存储方式,只存储非空的子空间编号。对相邻子空间进行3进制编码,利用编码的对偶关系,建立相邻子空间之间的指针连接,构造出包含KNN搜索所需的各类信息的广义表,然后再搜索KNN。KNN搜索过程中,在计算被测点到候选点距离时,直接删除筛选立方体内切球之外的点,可将参入按距离排序的候选点数减少为现有算法的一半。依赖K值和不依赖K值的分块原则,均可计算不同的K邻域。实验结果表明,该算法不仅具有低的存储器使用量,而且具有较高的效率。
关键词: 散乱点云    子空间    K最近邻    广义表    
K-nearest neighbor searching algorithm for laser scattered point cloud
ZHAO Jingdong, YANG Fenghua     
School of Mathematical Sciences, Qufu Normal University, Qufu Shangdong 273165, China
Abstract: Aiming at the problem of large amount of data and characteristics of surface in laser scattered point cloud, a K-Nearest Neighbors (KNN) searching algorithm for laser scattered point cloud was put forward to reduce memory usage and improve processing efficiency. Firstly, only the non-empty subspace numbers were stored by multistage classification and dynamic linked list storage. Adjacent subspace was coded in ternary, and the pointer connection between adjacent subspaces was established by dual relationship of code, a generalized table that contained all kinds of required information for KNN searching was constructed, then KNN were searched. In the process of KNN searching, the candidate points outside inscribed sphere of filtration cube were directly deleted when calculating the distance from measured point to candidate points, thus reducing the candidate points that participate in the sort by distance to half. Both dividing principles, whether it relies on K value or not, can be used to calculate different K neighborhoods. Experimental results prove that the proposed algorithm not only has low memory usage, but also has high efficiency.
Key words: scattered point cloud    subspace    K-Nearest Neighbors (KNN)    generalized list    
0 引言

K最近邻(K-Nearest Neighbors,KNN)搜索是提高点云去噪、点云简化以及三维重建效率的关键一步。目前的K邻域搜索算法分为三类:Voronoi图法、树法和分块法。其中Voronoi图法[1-2] 时间复杂度高,算法效率普遍较低。树法[3-6]在K邻域搜索中,当邻近点和候选点处于树的不同层时会造成搜索范围扩大,搜索效率降低。最优近邻算法虽可减少K-D树[6]的搜索范围,但在确定最优近邻算法的参数和寻找最佳路径时又会增额外开销,且一组参数的确定只适合于一个特定数据集的K邻域搜索。Voronoi图法和树法的共同优点是不需要连续的存储空间分配,但搜索效率低。算法简单且高效的是分块法[7-13]。其中,周儒荣等[7]的算法中不可避免地使用了较多的浮点运算和排序运算。熊邦书等[8]在文献[7]的基础上作了改进,并推广了文献[9]中只能解决2维问题的搜索方法,但是如果点云数据分布不均匀,会使某个子空间内的点特别多,在K邻域搜索时也会进行大量的浮点运算[13]。马骊溟等[10]通过对点云数据坐标进行3维排序,然后根据坐标点的3维序号计算3个方向的交集,最后再根据传统方法搜索邻近点,避免了大量的浮点运算,但增加了一次对整个点云数据的3维排序过程。文献[11]中以候选点到所对应立方体小栅格的环六壁的距离为球体半径的取值范围,通过逐步增加球体半径扩大搜索范围来进行K邻域搜索,但是该算法不适合K值较大的情况[13]。文献[12]借鉴了文献[4]提出的“凡在同一区域的采样点,它们邻域的搜索范围有着一定程度的重叠”的思想,采用二次划分策略,利用多个子空间结构逐步减小搜索范围以提高搜索速度,不过该算法的稳定性不是很理想。杨军等[13]针对文献[8, 10-12]中的不足,将整个点云空间分块,并将点归入各个相应的子空间,通过选取符合条件的坐标点,降低了点云密度不均匀的影响,还采用了子空间扩张机制来避免文献[12]中可能出现的死循环。然而该算法需要对子空间内的所有点进行3维排序,在扩大子空间编号时还需要将相邻子空间中的数据点归并到有序序列之中,需要较长的排序时间;在对候选点进行距离计算后的排序过程中,同文献[8-12]一样,使过多的、不可能为K邻域的点参入了排序过程。该算法的最大问题在于使用了过多对效率敏感的控制参数,严重影响了算法的通用性。

前述的分块算法,虽然存在一些不足,但比Voronoi图法和树法的效率均高。原因是现有的分块算法都是采用了分类效率最高的基数分类法对点云数据进行分块,在K邻域搜索中,当需要扩大子空间时,又都是通过基数分类时所采用的顺序存储结构的线性表来直接获得相邻子空间的编号。然而顺序存储结构的线性表需要连续空间,当散乱点的数量很大、K值很小时,会使子空间的个数激增,造成存储空间的大量浪费。由于实际中的点云边缘存在较大起伏,凹陷部分所对应的子空间内根本没有任何数据点;尤其是表面型的激光三维点云,由于其中心存在较大的空洞,空洞部分的子空间也没有任何数据点,这些不包含任何数据点的子空间编号出现在线性表中,必将造成大量的存储空间浪费。如果采用链表结构,则会避免存储空间的浪费。链表存储结构在K邻域搜索中需要子空间扩大时,如何做到不需要遍历链表也能直接获得其周围相邻子空间是问题的关键。

要采用分块思想来提高K邻域的搜索效率,既要避免使用大的连续存储器,又要不遍历链表也能直接获取周围相邻子空间,必须采用新的分块方式,并改变子空间编号的存储结构。为此,本文提出多级分块方式,只存储非空子空间的链式广义表存储结构;对广义表中的相邻子空间进行3进制编码,利用编码的对偶关系,快速建立相邻子空间之间的连接,使子空间扩大时能够直接获得相邻子空间的编号,避免了遍历子空间链表;为提高K邻域的搜索效率,在计算被测点到候选点距离的过程中,直接删除筛选立方体内切球之外的点,将参入按距离排序的点数减少一半。

1 算法

为了便于算法描述,先给出几点说明:

1) 将一次分块所形成的大块称为一次栅格。

2) 将一次栅格再次分块后所形成的小块称为二次栅格,将n-1次栅格再次分块后所形成的小块称为n次栅格,n次栅格相当于引言中所提到的子空间或小栅格。考虑到子空间划分次数的增加会降低整体算法的效率,且实际应用中,一般经过2次划分即可得到足够大的子空间数量,因此,本文以n=2为例说明子空间的划分算法。

3) 设被规格化到最小值为0的三维数据点,存放在以Data为头指针的链表中,线性表Tub为最终收集到的二次栅格(K邻域搜索中也称为子空间或小栅格)的头指针表,其索引编号为二次栅格的主编号,将建立了相邻子空间指针连接后的Tub称为最终生成的广义表,用于分块后的K邻域搜索。

4) 由于顺序存储结构的线性表具有数组的特性,因此,本文将所有顺序存储结构的线性表均用数组来表示。

1.1 算法的主体步骤与基本数据结构

算法的主体步骤分为两步:

第1步 子空间划分。对于不同的K值,本步骤可只进行一次,亦可根据K值的大小重新分配。它既能避免使用大的连续存储器,又能减少存储器使用量,也是在K邻域搜索中直接获取相邻子空间编号的关键。将点云空间进行一次划分生成一次栅格,再倒序对每一个非空的一次栅格进行二次划分生成二次栅格;取二次栅格分块个数最多的编号为栅格主编号,以主编号为索引建立线性指针表,表中的每个指针指向同一类二次栅格链表的表头,倒序收集非空的二次栅格,并按序插入到该类的链表中,再按索引编号由小到大分别建立相邻二次栅格之间的指针链接,构造出点云数据广义表。

第2步 K邻域搜索。采用步进筛选半径r,按坐标预筛选出2r为边长的候选点立方体。对筛选立方体中的点进行距离计算的过程中,直接删除筛选立方体内切球之外的点,只对剩余的内切球中的候选点按距离排序。

算法中的主要数据类型:

1) 三维数据点的节点类型。

typedef struct tagPOINT

{ unsigned int x,y,z; //三维数据点的坐标

struct tagPOINT *Next; //后向连接指针

}POINT;

2)二次栅格的节点类型。

typedef struct tagCUBE

{ union {int Cubes;int Points;}Num; //二次栅格的个数,或该栅格中的点数

int sx,sy,sz; //二次栅格的编号

struct tagCUBE *Next;//后向连接指针

struct tagPOINT *Content;//二次栅格中数据点的头指针

struct tagCUBE * Neighbors[27];//相邻栅格的头指针

}CUBE;

1.2 一次栅格与二次栅格的边长

借鉴参考文献[8]的分块原则,将整个点云数据的3个坐标规格化为最小值0,并保留3个坐标的偏移值(x0,y0,z0),则点云空间的范围为[0,xmax]、[0,ymax]、[0,zmax]。由立方体及其内切球的体积可得:

$2\times \frac{4}{3}\pi {{(\frac{{{L}_{2}}}{2})}^{3}}/K={{x}_{\max }}{{y}_{\max }}{{z}_{\max }}/n$

其中:二次栅格边长${{L}_{2}}\approx \sqrt[3]{{{x}_{\max }}{{y}_{\max }}{{z}_{\max }}K/n}$。令

$\left\{ \begin{align} & {{L}_{2}}=\left\lfloor \beta \sqrt[3]{{{x}_{\max }}{{y}_{\max }}{{z}_{\max }}K/n} \right\rfloor \\ & {{L}_{1}}={{L}_{2}}\times Ratio \\ \end{align} \right.$ (1)

其中: β为调整系数[8],用来调节二次栅格的边长;$\left\lfloor * \right\rfloor $表示下取整;K是一个较小邻域点的个数;n是三维数据点总数目;L1为一次栅格的边长;Ratio是一次栅格与二次栅格边长的比率,它确定了一个一次栅格可分成的二次栅格的数量为Ratio3

${{C}_{1}}=\left\lceil {{x}_{\max }}/{{L}_{1}} \right\rceil $, ${{C}_{2}}=\left\lceil {{y}_{\max }}/{{L}_{1}} \right\rceil $, ${{C}_{2}}=\left\lceil {{z}_{\max }}/{{L}_{1}} \right\rceil $,则一次栅格的数量为C1×C2×C3。令Cmax=max(C1,C2,C3),则栅格主编号为Cmax所指向的编号,定义指针类型的线性表POINT*Gride1 [C1] [C2] [C3]、 POINT*Gride2 [Ratio] [Ratio][Ratio]和CUBE*Tub[Cmax×Ratio] (在此用数组的形式表示顺序存储结构的线性表),分别用于存放一次栅格数据链表的头指针、二次栅格数据链表的头指针和以栅格主编号为索引的二次栅格链表的头指针。若一个指针占用4个字节,点云数据采用二次划分所需的三个连续空间之和为:4(C1×C2×C3+Ratio3+Cmax×Ratio)字节。在同样的情况下,如果使用传统的分块方式,则连续存储器分配的数量为: 4(C1×C2×C3×Ratio3)字节。若取Ratio=5,则传统的分块方式所需连续存储器的数量将是采用二次划分数量的125倍。

1.3 创建有序二次栅格链表

将点云空间进行一次划分生成一次栅格,再倒序对每一个非空的一次栅格进行二次划分,并倒序收集非空的二次栅格,按序插入到以栅格主编号为索引的链表中,最后再建立所有二次栅格与其相邻的26栅格的指针链接,为K邻域搜索作准备,在整个分块过程只改变三维数据点的指针,不移动数据点在内存中的位置。

1.3.1 一次栅格操作

初始化Gride1为空(NULL),遍历点云数据链表Data中的每一个数据点P,采用基数分类法将P归并到Gride1[Px/L1,Py/L1,Pz/L1]为头指针的链表中。归并时,只调整P的指针,不改变P的位置,从而避免了点云数据在内存中的移动。归并结束后Data链表变为空。

1.3.2 二次栅格操作

二次栅格操作是对每一个非空的一次栅格进行二次划分并回收的操作。由于二次划分后接着回收,且回收操作能够自动地将非空的二次栅格清空,因此,只需要在定义Gride2时将其清空即可,不必在每次二次划分前都要对Gride2清空。

i∈[0,C1-1],j∈[0,C2-1],k∈[0,C3-1],倒序遍历 Gride1,将Gride1[i][j][k]为非空头指针的一次栅格数据点,采用基数分类法归并到头指针为Gride2[ii][jj][kk]的链表中,并统计归并到该链表中的数据点数m[ii][jj][kk]。归并结束后Gride1[i][j][k]链表变为空。其中:

$\left\{ \begin{align} & ii=({{P}_{x}}-i\times {{L}_{1}})/{{L}_{2}} \\ & jj=({{P}_{y}}-j\times {{L}_{1}})/{{L}_{2}} \\ & kk=({{P}_{z}}-k\times {{L}_{1}})/{{L}_{2}} \\ \end{align} \right.$ (2)

每一个非空的一次栅格被二次划分到Gride2后,接着进行二次栅格的回收,回收过程是:令ii∈[0,Ratio-1],jj∈[0,Ratio-1],kk∈[0,Ratio-1],倒序遍历每一个二次栅格数据的头指针Gride2[ii][jj][kk],若非空,创建一个类型为CUBE的新节点,新节点的二次栅格编号为:

$\left\{ \begin{align} & sx=i\times {{L}_{2}}+ii \\ & sy=j\times {{L}_{2}}+jj \\ & sz=k\times {{L}_{2}}+kk \\ \end{align} \right.$ (3)

数据点的个数Num.Pointsm[ii][jj][kk],数据链表的头指针Content为Gride2[ii][jj][kk]。将新节点按sy、sz二维升序方式插入到Tub[sx]所指的链表中,并将Tub[sx]所指的头节点中的二次栅格数量Num.Cubes加1,其中Tub[sx]是以sx为索引的二次栅格链表的头指针(此处设二次栅格分块的主编号为sx)。

二次栅格操作的C语言伪代码如下:

For(i=C1-1 to 0,j= C2-1 to 0,k= C3-1 to 0)

{ if(Gride1[i][j][k]≠NULL) //只对非空的一次栅格进行二次划分

{ while(Gride1[i][j][k]≠NULL) //一次栅格的二次划分

{ 取Gride1[i][j][k]中的一个数据点P,按式(2)计算ii,jj,kk;

P归并到以Gride2[ii][jj][kk]为头指针的链表中;

m[ii][jj][kk] ++;

Gride1[i][j][k]指向下一个;

}

For(ii=Ratio-1 to 0,jj=Ratio-1 to 0,kk=Ratio-1 to 0) //回收二次栅格

{ if(Gride2[ii][jj][kk]≠NULL) //只回收非空二次栅格

{ 创建CUBE类型的新节点P1;

按式(3)计算P1的编号sx,sy,sz;

P1的Content=Gride2[ii][jj][kk];

P1的Num.Points=m[ii][jj][kk];

Gride2[ii][jj][kk]=NULL; m[ii][jj][kk]=0;//清非空的二次栅格指针和点数

P1按sy,sz二维有序(从小到大)插入到Tub[sx]指向的链表中;

Tub[sx]头节点的NUM.Cubes++;}}}}

二次栅格回收后的Tub图 1所示,其中每一行的二次栅格是按栅格编号关键字sy、 sz由小到大排列的。

图 1 有序二次栅格Tub的示意图
1.3.3 相邻栅格链接

K邻域搜索中,为了能够直接获得相邻栅格的编号及存储位置,必须建立相邻栅格的指针连接,使广义表Tub包含K邻域搜索所需的各种信息。

三维空间中的每一个二次栅格最多有26个相邻栅格,将相邻栅格的指针存入该栅格的Neighbors[f]中,其中f是相邻栅格的3进制编码。由于二次栅格是有序排列的,因此可在较小的范围内找到其相邻栅格,再利用相邻栅格编码的对偶关系,可使查找减半。

1) 相邻栅格编码与编码的对偶性。

设当前栅格A的编号为(sx,sy,sz),则其相邻栅格B的编号为(sx+i,sy+j,sz+k),其中i,j,k∈{-1,0,1},若以A为参考点,对B进行3进制编码,则BA中的编码为:

$\begin{align} & {{f}_{B}}~=\left( sx+i-sx+1 \right){{3}^{2}}+\left( sy+j-sy+1 \right){{3}^{1}} \\ & +\left( sz+k-sz+1 \right){{3}^{0}}=13+9i+3j+k \\ \end{align}$ (4)

若以B为参考点,则AB中的编码为:

$\begin{align} & {{f}_{A}}=\left[ sx-\left( sx+i \right)+1 \right]{{3}^{2}}+\left[ sy-\left( sy+j \right)+1 \right]{{3}^{1}}+ \\ & \left[ sz-\left( sz+k \right)+1 \right]{{3}^{0}}=13-9i-3j-k \\ \end{align}$

由于fB+fA=26,所以

${{f}_{A}}=26-\text{ }{{f}_{B}}$ (5)

2) 相邻栅格的链接过程。

Tub[m]链表中,栅格的编号sx均等于m,且编号sy,sz是由小到大二维有序的。设任意一个栅格A的编号为(m,sy,sz),则与A相邻的栅格只能出现在Tub[m-1]、Tub[m]和Tub[m+1]链表中。其中Tub[m-1]链表包含A的相邻栅格编码为0~8,Tub[m]中包含A的相邻栅格编码为9~17(其中13表示A自身),Tub[m+1]中包含A的相邻栅格编码为18~26。若ATub[0],则A的相邻栅格编码无0~8,若A为Tub[m]中的第一个,则A的相邻栅格编码无9~12。若相邻栅格搜索从Tub[0]中的第一个栅格开始,按行、再按列向后推进,则ATub[m]的相邻栅格搜索过程只包括m行搜索和m+1行搜索。这是因为Am-1行中的相邻栅格已经在m-1行搜索时,通过相邻栅格编码的对偶关系建立了链接,如图 2中的虚线栅格对编号(1,8,2)的链接。同样道理,在当前行m中进行相邻栅格搜索时,也只需要搜索A之后的相邻栅格,如图 2中的实线链接,A之前的点划线链接已在(1,7,1)、(1,7,2)、(1,8,1)的相邻栅格搜索时被建立。

m行搜索:ATub[m],A之后的相邻栅格编码14~17只能出现在m行上,且相邻栅格与A的编号之差的绝对值均不超过1,如果有这样的二次栅格B,按式(4)计算fB,将B的地址写入ANeighbors[fB]中,同时将A的地址写入B的Neighbors[26-fB]中。对于图 2中的(1,8,2)栅格,m行搜索只需建立与(1,9,1)、(1,9,2) 栅格的实线链接。

提前退出m行搜索的条件是:

$\begin{align} & B.sy-A.sy>1或 \\ & \left( B.sy-A.sy==1\text{ }且\text{ }B.sz-A.sz>1 \right) \\ \end{align}$ (6)

② m+1行搜索:在Tub[m+1]中如果有二次栅格,依次取一个二次栅格B,若BA的编号之差的绝对值均不超过1,按式(4)计算fB,将B的地址写入ANeighbors[fB]中,同时将A的地址写入BNeighbors[26-fB]中。提前退出本行搜索的条件仍然为式(6)。如图 2中,对于(1,8,2)栅格在m+1行上的搜索只需要建立与(2,7,2)、(2,8,1)、(2,8,2)、(2,9,1)和(2,9,3)的实线链接。

图 2 编号为(1,8,2)的相邻栅格链接示意图

相邻栅格链接操作的C语言伪代码如下:

For(m=0 to Cmax×Cmax-1)

{ if(Tub[m]≠NULL)

{ 依次取ATub[m],直到Tub[m]的最后一个;

A进行m行搜索和m+1行搜索; }}

图 3 筛选立方体与内切球的差集
1.4 K邻域搜索

Pi的K邻域搜索是先确定合适的筛选立方体,再删除筛选立方体与筛选立方体内切球的差集。如果剩下的内切球中的候选点数等于K,则K邻域搜索成功;如果大于K,则按距离由小到大排序,排序后的前K个点就是PiK邻域。所谓筛选立方体就是以Pi为中心、2倍筛选半径为边长的小立方体。如图 3中的虚线框是筛选立方体,它与内切球的差集是图中的方形点。由于立方体的体积约等于其内切球体的2倍,所以删除方形点后只对圆形点按距离排序,可使参入按距离排序的点数减少一半。

1.4.1 最大筛选半径

对于任意一个二次栅格A,当PiAK邻域在A中搜索时,其最大筛选半径是PiA边界的最小距离;若Pi的K邻域搜索扩大到A的相邻子空间时,其最大筛选半径是PiA的相邻子空间外边界的最小距离,同样的方法,最大筛选半径可以继续向外扩大。如果A的某侧无相邻子空间,或是整个点云的边界,则该侧的外边界为+∞,如图 4给出了几个最近相邻子空间的最大筛选半径,其中的深色框为A,浅色框为非空相邻子空间。图 4(a)(b)的实线为相邻子空间的有限外边界,其他外边界均为+∞,R为最大筛选半径,最大筛选立方体是以2R为边长的小立方体,如图 4(a)(b)中的虚线框。如果A无任何相邻子空间(即A为孤立子空间)如图 4(d),其最大筛选半径为+∞,表示A中的数据点均为候选点。

图 4 相邻子空间的外边界示意图
1.4.2 K邻域搜索过程

对于任意一个子空间(二次栅格)的K邻域搜索过程是:先将子空间中的数据点复制到一个数组中,如果数组中的点数大于2K,则先在自身范围内进行K邻域搜索,搜索成功的点从数组中删除,不成功的点留作扩大子空间范围后再搜索,直到数组中的所有点被删除。

设数组A用于存放当前子空间的数据点,nA中的点数目;数组B用于存放当前搜索范围的数据点,B中的每个元素除了包含点的三维坐标外,还包含一个距离数据项,用于记录该点到Pi点的距离;数组C存放符合筛选条件的B序号,m统计C中的点数目。若PjB,RPi的当前最大搜索半径,Pi在当前搜索范围内的K邻域搜索过程如下:

1) 初始化为m =0,筛选半径r=R/2-Δr,其中Δr = R/6

2) 若m<aK(a是预筛选点数目的最小控制阈值,一般取2)且r<R,则rr,筛选B中符合条件的数据点到C,直到maK,即:将满足|Pix-Pjx|≤r且|Piy-Pjy|≤r且|Piz-Pjz|≤r条件的点PjB中的序号依次放入数组C中,并累计C中点的数目m

3) 计算B序号在C中的所有点到Pi的距离,并删除C中满足|PiPj|>rB序号,每删除一个m减一,即删除筛选立方体与内切球的差集(删除操作无需移动C中的B序号,只要用C中的最后一个B序号覆盖当前的B序号,C的指针后退一步,继续计算|PiPj|)。

4) 若mK,则K邻域搜索成功:如果m>K,按|PiPj|由小到大排序C中的B序号,否则不用排序,取C中的前KB序号为PiK邻域,删除A中的Pi(删除过程也是将A中的最后一个覆盖当前的Pi),A中的点数目n-1;

m<Kr<R,则rr,再次筛选B中数据点到C,直到mK,使K邻域搜索成功;或r达到R时,仍有m<K,则退出当前搜索范围内的K邻域搜索。

2 实验与分析

由于Voronoi图法的复杂度高且效率很低,目前已不被采用。而树法的搜索效率均低于分块法,在文献[12]中,已通过对脚、小孩、牙齿和球套四个实例分别用八叉树算法[4]、K-D树算法[6]和分块算法进行了测试比较,详细情况可见文献[12]中的表 3,其结论是:分块算法的效率均高于树法,因此树法也不作为本文效率的比较对象。在分块法中,本文选择了文献[8, 10, 13]作为效率比较对象,其中杨军算法[13]不仅稳定且效率高,为此将其作为本文的主要比较对象。在以下的所有实验中,对杨军算法所使用的初始搜索步长l、步长增量Δl、左右控制阈值αβ,均按文献[13]的要求,经过了多次筛选,并取效率最高的一组作为实验参数。所有参与测试的算法均用VC++ 6.0语言实现,并在启天M1750商用计算机(Pentium Dual-core CPU,主频3.06GHz,内存2GB)上进行,测试统计的时间是被测对象所有点的K邻域搜索时间。

2.1 实际数据验证

用激光扫描分别获取茶壶、兔子和马模型的3组点云数据(stl格式文件),经反相处理后的点云图像如图 5所示,点云的点数分别是13552、21975和103754。点云数据文件存放在当前计算机的硬盘上,测试过程中,首先将被测点云数据读入内存,并以链表的形式存储,然后使用内存中的点云数据,针对不同的K值,分别用文献[8]算法、文献[10]算法、杨军算法[13]和本文算法计算出所有点的K邻域。为便于效率比较,所有算法均采用的是依赖于K值的分块算法、K邻域的计算时间(不包括数据从硬盘中的读取时间),对比结果如表 1所示。因为文献[8, 10]算法和杨军算法的分块数相同,所以在后面统计时统称为传统分块算法。表 2~3即为传统分块算法与本文算法的实际分块个数和分块时间对比。由于传统分块算法的分块都是整个点云包围盒中的所有块,因此没有给出它们的分块图,只给出了本文算法的分块结果(如图 6所示)。本文所采用的分块算法需要遍历2次点云数据,还要建立相邻子空间之间的指针连接,其分块时间必定大于传统分块算法的分块时间,表 3中给出的是依赖K值的分块时间。表 4是本文算法中,采用依赖K值的分块原则的子空间剖分时间和分块后的K邻域搜索时间对照。由于本文的分块算法在K邻域搜索中能够直接获得相邻子空间,因而也可以先将点云进行细分,用细分后的子空间同样可以快速计算不同的K邻域。表 5给出了不依赖K值与依赖K值分块原则的搜索时间对照,其中不依赖K值搜索算法的所有子空间边长均是用式(1)取K=10时计算的。

表 1 四种算法K邻域搜索时间对比
图 5 经反相处理后的点云图像
图 6 本文算法的实际分块结果
表 2 传统分块算法与本文算法的实际分块个数对比
表 3 传统分块算法与本文算法的实际分块时间对比
表 4 本文算法的子空间剖分时间与剖分后的K邻域搜索时间
表 5 本文算法不依赖K值与依赖K值分块原则的搜索时间
2.2 实际数据分析

表 1的K邻域搜索时间统计来看,杨军算法是传统分块算法中效率最高的一个,对于较小的K邻域更为有效,当K值逐渐变大时,其效率逐渐接近于文献[8]算法和[10]算法,但均高于文献[8]算法和[10]算法。进一步观察表 1中最后一列的时间比率可以看出,本文算法的效率约是杨军算法的2.5倍,且点云数量越大(马的点云数量最大)、K邻域越小时效率越高。这是因为本文算法在扩大K邻域搜索范围时,能够方便地取得相邻子空间编号;在最终确定K邻域时,能够将参入按距离排序的点数减少一半,使算法的效率显著提高。

表 2的分块个数比较来看,本文算法的分块数量远小于传统分块算法,尤其对于边缘起伏较大的马模型,当K邻域为10时,传统分块算法的分块数量为10625块,是本文算法的4.26倍。这是因为本文算法所得的块只是那些真正包含数据点的块,如图 6所示,而传统分块算法所得的块是整个点云包围盒中的所有块,其中包含着大量的空块。由于实际的三维点云表面一般都是非规则的,尤其是激光三维点云,其内部又是空的,因此本文算法的分块数量必定远小于传统分块算法的分块数量。从表 2中不同模型的分块比率可以看出,马模型的比率明显大于茶壶的,其主要原因是马模型包围盒中有很大部分的无点区(主要集中在马的头下和腹下)。如果点云数据充满包围盒,且模型内部也含有大量的数据点,则本文算法的分块数量将等于传统分块算法的分块数量,此时,所对应的模型一定是长方体实体,而实际中完整的长方体点云少见。观察表 2中不同K值的分块比率可知,本文算法更适合于点云数量较大或分块数量较多的情况,当点云数据面型特征并不明显且边缘规整时(如茶壶模型),本文子空间划分方法的性能会略有下降,但并不明显。这是因为,本文所采用的分块算法需要遍历2次点云数据,还要建立相邻子空间之间的指针连接,其分块时间会有所增加(见表 3),但由于分块过程中只调整数据点的指针,不改变数据点的存储位置,从而避免了三维数据点在内存中的移动,节约了分块时间,且实际得到的分块个数大量减少(其中不包含任何空块),又在建立广义表时采用了各种高效率的处理方式,使得子空间的剖分时间同样很小,在整个K邻域搜索算法中只占很小的比例,与分块后的K邻域搜索时间相比可以忽略不计(如表 4的子空间剖分时间与剖分后的K邻域搜索时间比率不大于0.362%),不会造成整体算法的效率降低。对于一组点云,若采用同一分块来计算不同的K邻域,即采用不依赖K值的分块原则来计算K邻域,本文的分块算法同样有效,表 5的统计结果进一步证明了这一点。这是因为虽然这种不依赖K值的分块原则在K邻域计算时,会增加子空间的搜索范围,但由于所有的相邻子空间之间具有指针链接,能够直接获得相邻子空间编号,所以整体效率基本不变,仍可保持在98%以上。

3 结语

针对大规模激光三维点云K邻域搜索中,利用分块算法存在的连续存储器使用量大、存储空间浪费、对于大数据量模型可能导致存储空间分配失败,以及搜索效率低的问题,本文提出了多级分块方式和动态链式存储结构,既避免了使用大的连续存储器,又使子空间编号的存储数量大量减少,从而降低了存储器的使用量。通过以子空间主编号为索引建立子空间链表的头指针表、逐级倒序收集非空子空间,以较少的时间开销,避免了子空间编号的排序过程,同时也避免了三维数据点在内存中的移动过程,使整个子空间的组织占用类似于树的所需空间;利用相邻子空间编码的对偶关系,建立相邻子空间之间的连接,避免了扩大子空间编号时所需的子空间链表遍历,使获取相邻子空间编号的过程类似于普通分块方式的获取过程;在K邻域搜索中,采用按坐标筛选的方式筛选出候选点,避免了子空间内的3维排序;在计算当前点到候点距离的过程中,直接删除筛选立方体内切球之外的点,将参入按距离排序的候选点个数减少一半。实验结果表明,该算法具有较高的效率,子空间编号的数量也大量减少。由于本文剖分的子空间包含K邻域搜索所需的各类信息,在子空间之间建立了指针连接,且耗时小,其分块时间与分块后的K邻域搜索时间相比可以忽略不计,因而可以采用依赖K值或不依赖K值的分块原则来计算不同的K邻域,且具有同等的效率。

参考文献
[1] GOODSELL G. On finding p-th nearest neighbors of scattered points in two dimensions for small p[J]. Computer Aided Geometric Design, 2000, 17 (4) : 387-392. (0)
[2] DICKERSON M T, DRYSDALE R L S, SACK J R. Simple algorithms for enumerating inter point distance and finding K nearest neighbors[J]. International Journal of Computational Geometry and Applications, 1992, 2 (3) : 221-239. (0)
[3] CONNOR M, KUMAR P. Fast construction of K-nearest neighbor graphs for point clouds[J]. IEEE Transactions on Visualization & Computer Graphics, 2010, 16 (4) : 599-608. (0)
[4] SANKARANARAY J, SAMET H, VARSHNEY A. A fast all nearest neighbor algorithm for applications involving large point-clouds[J]. Computers & Graphs, 2007, 31 (2) : 157-174. (0)
[5] ARYA S, MOUNT D M, NETANYAHU N S, et al. An optimal algorithm for approximate nearest neighbor searching fixed dimensions[J]. Journal of the ACM, 1998, 45 (6) : 891-923. (0)
[6] MUJA M, LOWE D G. Scalable nearest neighbor algorithms for high dimensional data[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2014, 36 (11) : 2227-2240. (0)
[7] 周儒荣, 张丽艳, 苏旭, 等. 海量散乱点的曲面重建算法研究[J]. 软件学报, 2001, 12 (2) : 249-255. ( ZHAO R R, ZHANG L Y, SU X, et al. Algorithmic research on surface reconstruction from dense scattered points[J]. Journal of Software, 2001, 12 (2) : 249-255. ) (0)
[8] 熊邦书, 何明一, 余华璟. 三维散乱数据的K个最近邻域快速搜索算法[J]. 计算机辅助设计与图形学学报, 2004, 16 (7) : 909-912. ( XIONG B S, HE M Y, YU H J. Algorithm for finding K-nearest neighbors of scattered points in three dimensions[J]. Journal of Computer-Aided Design & Computer Graphics, 2004, 16 (7) : 909-912. ) (0)
[9] PIEGL L A, TILLER W. Algorithm for finding all K-nearest neighbors[J]. Computer-Aided Design, 2002, 34 (2) : 167-172. (0)
[10] 马骊溟, 徐毅, 李泽湘. 基于动态网格划分的散乱点K邻域快速搜索算法[J]. 计算机工程, 2008, 34 (8) : 10-11. ( MA L M, XU Y, LI Z X. Fast K-nearest neighbors searching algorithm for scattered points based on dynamic grid decomposition[J]. Computer Engineering, 2008, 34 (8) : 10-11. ) (0)
[11] 刘越华, 廖文和, 刘浩. 逆向工程中散乱点云的K邻域搜索算法研究[J]. 机械设计与制造, 2012, 1 (3) : 256-258. ( LIU Y H, LIAO W H, LIU H. Research of K-nearest neighbors search algorithm in reverse engineering[J]. Machinery Design & Manufacture, 2012, 1 (3) : 256-258. ) (0)
[12] 赵俭辉, 龙成江, 丁乙华, 等. 一种基于立方体小栅格的k邻域快速搜索算法[J]. 武汉大学学报(信息科学版), 2009, 34 (5) : 615-618. ( ZHAO J H, LONG C J, DING Y H, et al. A fast algorithm for finding k-nearest neighbors based on small cube grids[J]. Geomatics and Information Science of Wuhan University, 2009, 34 (5) : 615-618. ) (0)
[13] 杨军, 林岩龙, 王阳萍, 等. 大规模散乱点的K邻域快速搜索算法[J]. 中国图象图形学报, 2013, 18 (4) : 399-406. ( YANG J, LIN Y L, WANG Y P, et al. Fast algorithm for finding the K-nearest neighbors of a large-scale scattered point cloud[J]. Journal of Image and Graphics, 2013, 18 (4) : 399-406. ) (0)