1 引言

随着计算机软、硬件的更新和发展,以及网络技术的不断成熟,数字博物馆技术^[1]逐渐成为研究热点。数字博物馆具有一些特有功能,如存储数字化和资源共享等。数字博物馆技术成为文化遗产和自然遗产虚实展示的重要技术手段。通过硬件设备和网络快速获得博物馆展品已经成为研究热点。其中,精确有效地对点云数据进行简化和压缩成为虚实展示的关键技术^[2]。通过三维(3D)扫描仪^[3]获得的数据非常密集,数据量庞大。随着3D扫描设备精度的提高,研究人员所获得的数据量增加,数据冗余现象越来越严重。考虑到数字博物馆的存储技术数字化以及资源共享的功能,如果不对原始数据进行简化、直接建模,大量的时间和资源将会被消耗掉,同时网络传输的负担会不断增加,数字博物馆的用户体验也会变差。

常用的三维模型简化方法可以分为两类:基于网格的简化^[4]和基于点云的简化^[5-6]。基于网格的简化需要构造多边形网格,该方法的缺点是需要网格拓扑相互一致,并且其网格生成复杂。基于点云的简化方法不仅可避免存储点的连接信息,减少内存消耗,加快简化速度,还能够保证重建表面的质量。因此,基于点云的简化方法更具优势,已成为研究的热点。

三维点云数据的简化可以分为两大类:基于模型拓扑结构的点云简化和基于特征保留的点云简化。Weir等^[7]依据体包围盒对点云进行简化,将大包围盒不断拆分成多个大小近似一致的小包围盒,并用小包围盒最靠近中心的点代表整个包围盒中所有点。该方法简单高效,但难以保证简化模型的精度。Shen等^[8]提出一种特征保留的三维点云简化算法:1) 根据点云密度进行下采样;2) 在点云下采样后计算平均曲率,得到点云的特征;3) 使用区域增长聚类方法根据曲率阈值得到简化点云。该方法不但需要消耗大量时间,而且无法获取更多的全局特征。李仁忠等^[9]提出了一种散乱点云的均匀精简算法,通过创建一个k邻域三维体素栅格,同时结合包围盒法对输入的点云数据进行k邻域距离计算和法线估计,以此来确定每个小立方栅格的重心,并用其近似表示小立方栅格内所有数据点。该方法虽然可以达到精简点云的目的,但是在曲率变化明显的区域损失了过多的特征点。Han等^[10]提出一种基于法向量的边缘保留点云简化算法,但该方法在平坦区域容易产生孔洞。张雨禾等^[11]提出一种基于局部特征检测算子和潜在曲面平均弯曲度的散乱点云简化算法,在较好地保留原始点云轮廓特征的同时可以简化点云,降低了点云简化误差,但该方法的检测特征点计算量大,比较耗时。Mahdaoui等^[12]提出一种基于密度函数和香农熵估计的点云迭代简化方法,该算法适用于不同密度的不同点云,简化点云的重建效果接近原始点云模型。袁小翠等^[13]提出了一种特征点保留的三维点云简化算法,该方法不但可以减少点云,而且能够较好地保持原始曲面的局部和全局特征;但是,该方法在对一些非封闭曲面的特征进行简化时容易失去边界特征。总体来说,现有基于点云的简化方法适用于曲率变化不大、外貌轮廓不复杂的点云模型,并且现有点云简化方法所得的结果无法保留其全局或局部特征信息。

为此,本文提出一种特征保留的三维点云简化方法。通过提取边缘点尽可能地保留点云的全局特征;使用期望最大化算法根据点的局部分布对点云进行聚类,计算平均曲率,确定局部特征点;对非特征点进行边缘感知的点云简化。

2 三维点云轮廓提取

点云简化问题可以概括为:假设原始点云模型为M,求M的一个子集P,在满足简化比的条件下,使得子集P在删除冗余点的同时能够尽可能地保留全局特征和局部特征。简化比可表示为

η=P^M^,(1)

式中: P^为简化点云模型中的点云个数; M^为原始点云模型中的点云个数。

三维点云模型简化的目的在于简化点云冗余数据的同时保留模型的几何特征。一般情况下,边缘点云比非边缘点云具有更多的特征,因此应在点云简化过程中保留这些边缘点云处的尖锐点。为此,可通过构造八叉树建立每个点的空间拓扑关系,从而找出每个点的邻近点,逐步检测并提取边缘点^[14]。

2.1 k邻点的建立与法向量的估计

点云通常是海量且无序的,故需要构造最小空间包围盒来划分空间,并使用八叉树搜索k个相邻点^[15]。空间包围盒的大小为

Csize=max(Δx,Δy,Δz),(2)

式中:Δx=x_max-x_min,Δy=y_max-y_min,Δz=z_max-z_min,其中x_max和x_min,y_max和y_min,z_max和z_min分别是点云x,y,z的最大值和最小值。空间包围盒是八叉树的根节点,将空间包围盒平均分成8个小立方体,8个立方体是根节点的叶节点。根据此规则划分每个立方体,直到达到指定的分割级别,如图1所示。八叉树的结构如图2所示。

图 1. 八叉树的空间分割

Fig. 1. Space segmentation of octree

下载图片 查看所有图片

图 2. 八叉树结构

Fig. 2. Octree structure

下载图片 查看所有图片

因此,八叉树具有以下特征:如果八叉树非空,那么树中的任何一个节点都没有子节点或有8个子节点。如果空间包围盒被划分为n层,则可以构造一个具有n+1层的八叉树,每个立方体对应一个唯一节点,每个节点的空间位置编码Q为

Q=qn-1…qm…q1q0,(3)

式中:q_m是八进制数,m∈{0,1,…,n-1},表示n层中的任意一层。如果q_m表示该节点在其兄弟节点之间的序号,则q_m+1表示q_m节点的父节点在其同胞兄弟节点间的序号。这样,从q₀到q_n-1完整地表示出八叉树中的每个叶子节点到树根的路径,如图2中黑色节点之间表示的路径,其黑色叶子节点的编码为{12,67}。

搜索节点p的k邻点(近邻点分为点近邻、边近邻和面近邻),即搜索p的叶节点最近的k个节点和p的相邻叶节点,也就是搜索立方体中的节点(节点p所在的立方体)及其26个相邻立方体,故使用八叉树结构可极大地提高相邻点的搜索速度。

由于某些点云可能不包含法向量信息,故需要计算每个点的法向量。计算点云法向量常用的方法有基于Delaunay方法和基于回归方法。本文使用经典的主成分分析(PCA)回归方法来计算点云的法向量。1) 通过最小二乘法将任意点及其相邻点拟合到局部平面;2) 用拟合平面代替原来的局部曲面;3) 点云中某点的法向量可以根据其拟合平面的法向量计算而得。在上述步骤中,点云中的点p及其k个相邻点的拟合平面F为

F(n,d)=argmin(n,di)∑i=1knp--di)2,p-=1k∑i∈Nppi,(4)

式中:n是平面F的法向量;d_i是点p邻域点到坐标原点的距离; p-是质心;N_p是点p的k个相邻点的集合;i为点序号;p_i是点p邻域中的任意邻近点。平面F包含点p,法向量n满足|n|²=1,因此该问题简化为一个半正定协方差矩阵U的特征值分解,U的表达式为

U=1k∑i=1k(pi-p-)(pi-p-)T,(5)

其最小特征值对应的特征向量即为点p的法向量估计。

2.2 边缘点检测

边缘点检测的步骤为:1) 利用点和法向量的拓扑关系构造最小二乘平面;2) 将拟合平面附近的每个点投影;3) 根据投影的均匀性提取边缘点。均匀性投影的判断标准可以是角度偏差和角度间隔,但二者的计算复杂度高。本文抛弃角度比较的思路,采用坐标值比较。

假设p是点云中的某点,使用八叉树搜索其k近邻,构造最小二乘平面并将邻近点投影到平面,创建三个平面xpy、xpz和ypz,分别平行于平面xoy、xoz和yoz。通过点p,将平面xpy、xpz和ypz作为参考平面,当参考平面两侧的投影数量之差与参考平面两侧投影数量之和的比大于指定阈值时,点p被判断为边缘点。以点p和平面xpy为例,如图3所示,平面xpy上方有11个点,平面下方有4个点,因此差值与和的比值为7/15。

图 3. 投影点附近点示意图

Fig. 3. Schematic of neighbor points near projection point

下载图片 查看所有图片

3 局部分布的聚类

聚类方法旨在确定小均匀区域中的点集并使点集在空间上相似,在不显著丢失几何特征信息的情况下去除冗余点。聚类过程使用考虑点集的局部分布的期望最大化(EM)算法^[16]。

EM算法^[17]是一种基于概率的聚类算法,类似于k均值方法。 k均值方法是将n个对象分成k个簇(簇内相似度较高,簇与簇之间相似度较低)。EM算法扩展了这种基本方法,允许更均匀的聚类生成。通常,EM算法在每次迭代中具有两个步骤——期望和最大化。在期望步骤中,假设观察到潜在变量,利用算法估计可能性的期望,即对每一个数据点计算其所属聚类的概率,作为权重。在最大化步骤中,计算最大似然值,最大化先前估计的预期可能性,即利用上述计算的权重来估计每个聚类的均值和方差,进而求取聚类的总体概率或极大似然。新的参数值用于后续的交互,直到算法收敛。

给定变量数据集{x₁,x₂,…,x_n}和迭代混合模型的初始参数值(w_0,c,μ_0,c),其中,c是k个初始簇的一个。基本的迭代分组算法步骤如下:

1) 在第j次迭代中计算每个簇中每个点的成员概率。

2) 更新混合模型参数。

3) 验证是否达到收敛,如果不是,则继续步骤1中的第j+1次迭代。

在第j次迭代中,对于任意一点x_i,簇c中的成员概率为

p(c|xi)=wj,cpj(xicpj(xi),(6)

式中:p_j表示第j次迭代中样本所属簇的概率。

如上所述,新参数值为

wj+1,c=1n∑i=1np(c|xi),(7)μj+1,c=∑i=1nxi·p(c|xi)∑i=1np(c|xi),(8)∑0,c=∑i=1np(c|xi)(xi-μj+1,c)(xi-μj+1,c)T∑i=1np(c|xi)。(9)

对于收敛验证,使用在第j和j+1次迭代中混合模型的对数似然之差。在第j次迭代中的对数似然度为

Ej=∑i=1nlg[pj(xi)]=∑i=1nlg∑c=1kwj+1,c·pj(xic。(10)

如果|E_j-E_j+1|≤ε,则达到收敛,否则计算新的迭代,其中ε表示给定迭代误差阈值。需要注意的是:对于簇的规模来说,随着簇规模的增加,时间以非线性的方式增加。在整个过程中,簇计算是最耗时的阶段,简化步骤是较快的阶段。原始点云和EM聚类簇分别如图4、5所示。

图 4. 原始点云

Fig. 4. Original point cloud

下载图片 查看所有图片

图 5. EM聚类簇

Fig. 5. EM cluster

下载图片 查看所有图片

4 曲率计算

点云的曲率^[18-19]可以反映局部表面的变化程度,它在几何变换(例如平移和旋转)下保持不变,因此也被称为不变量。从几何意义上来讲,曲率是切向量变化的速率。若S(r,t)表示一个常规的连续的三维参数曲面,S上在S(r₀,t₀)处的主曲率k₁(r₀,t₀)和k₂(r₀,t₀)被定义为该点曲率的最大值和最小值。根据欧拉定理,曲面S(r,t)切平面方向上法线曲率为

kn(θ)=k1cos2θ+k2sin2θ,(11)

式中:θ是第1主方向与切线方向上的夹角。高斯曲率K和平均曲率H是由曲面的主曲率唯一确定的,即

K=k1·k2H=(k1+k2)2。(12)

高斯曲率的主要限制是当任何主曲率为零时其等于零。平均曲率的问题在于极小曲面总是具有H=0。另一个有利的度量是绝对曲率A,其表达式为

A=|k1+k2|。(13)

因此,无论曲率类型如何,人们都可以使用最小二乘回归算法进行局部曲面拟合,在给定点拟合曲面,计算曲率K、H和A。曲率可以用于描述局部几何变化。因此,为了避免简化尖锐的特征点,在估计点云的曲率时,需要标记估计曲率超过固定阈值的曲率。阈值集通过计算该点所属簇内局部平均曲率的均值来获得。具体来说,估计给定聚类中所有点的平均曲率A_i,计算平均曲率的均值 Ai¯,并满足以下规则:

1) 如果A_i> Ai¯,表示在点i处发生表面的显著区域,将此点标记为特征点。

2) 如果A_i< Ai¯,表示在点i处是没有尖锐特征的平滑区域,此点可以在下一阶段减少。

图 6. 曲率估计选择的特征点和非特征点。(a)原始点云;(b)曲率估计点云

Fig. 6. Feature points and non-feature points selected by curvature estimation. (a) Original point cloud; (b) curvature estimation point cloud

下载图片 查看所有图片

图6显示了原始模型及其曲率估计图(轮廓点和特征点),在简化点云模型过程中,保存这些点可以保留边界和几何特征。因此,下一阶段只使用非特征点进行简化。

5 基于Hausdorff距离的点云简化算法

5.1 有向Hausdorff距离

在数学中,Hausdorff距离^[20-21]测量两个点集之间的距离,可以描述两组点集之间的相似程度。假设有两组集合X={x₁,x₂,…,x_n},Y={y₁,y₂,…,y_m},则这两个点集之间的Hausdorff距离定义为

H(X,Y)=maxl(X,Y),l(Y,X),(14)l(X,Y)=maxxi∈xminyi∈yx-y,(15)

图 7. 有向Hausdorff距离

Fig. 7. Directed Hausdorff distance

下载图片 查看所有图片

式中:H(X,Y)称为双向Hausdorff距离,l(X,Y)和l(Y,X)分别称为从X集合到Y集合和从Y集合到X集合的有向Hausdorff距离。使用有向Hausdorff距离(DHD)进行点云简化,如图7所示。假定原始点云{x₁,x₂,x₃,y₁,y₂},并且N={y₁,y₂}表示已经采样的特征点,使用DHD从M中选择一个新的特征点加入N。此时,(15)式则可以理解为:1) 分别计算从x₁、x₂和x₃到y₁的欧几里德距离,并保存最近的值d₁₁;2) 与第1步相同,分别计算从x₁、x₂和x₃到y₂的距离,并保存最近的值d₃₂;3) 比较并选择d₁₁和d₃₂之间的最大值。显然,d₃₂>d₁₁,因此,将x₃作为新的采样点,并且将N更新为N={y₁,y₂,x₃}。

DHD方法是一种从随机种子点开始的增量区域生长方法。本质上来讲,DHD是一种基于距离的自适应加权方法。与其他简化方法不同的是,在DHD中,当前采样点的位置对如何从原始点云中选择新特征点起着至关重要的作用。因为点云的边界点与内点的距离相对较远,所以DHD方法对轮廓和角点的采样是非常有效的。

5.2 点云简化

综上所述,点云简化算法步骤如下:

1) 通过构造八叉树搜索每个点的k近邻点并计算每个点的法向量。

2) 检测边缘点并保留。

3) 使用期望最大化算法,根据点的局部分布对点云进行聚类。

4) 确定高曲率的点(不会被删除的特征点)。

5) 边缘感知的DHD点云简化。

6) 合并步骤2、步骤4和步骤5获得的采样点,删除重复点。

6 实验结果和分析

选取斯坦福兔子、椅子和兵马俑(TW)点云模型作为实验对象。斯坦福兔子模型包含弯曲度较大的曲面,椅子模型包含尖锐特征和较平坦区域,兵马俑模型包含曲面和较复杂的特征。为了验证本文点云简化算法的有效性,从简化效果、重建效果和简化速度三方面进行评估,并与曲率自适应方法和K均值聚类简化方法分别进行对比实验。

6.1 简化效果

不同简化率η下的简化效果如图8所示,可以看出,本文算法在有效提取尖锐特征点(簇中高平均曲率的点)的基础上进一步检测并保留了边缘点。两者结合,在简化点云的过程中更好地保留了不同点云模型的轮廓,使三维点云模型的体貌展现得更为真实全面。

6.2 重建效果

在简化阶段之后,将特征点集和简化数据集连接起来并去除重复点,形成简化后的点云。图9显示不同点云模型在不同简化率下的重建效果(为了便于展示,对重建模型进行了三角剖分、封装和渲染)。可以看出,随着简化率的减少,简化后的封装模型会失去一些纹理信息和局部细节特征,但封装简化后的模型基本保留了原始模型的几何特征和轮廓外貌(在简化率20%以上)。重建效果进一步证明,在不同模型中,基于边缘感知的DHD点云简化有效地保留了边界和角点等几何特征。

图 8. 不同点云模型的简化效果。(a)斯坦福兔子;(b)椅子;(c)兵马俑

Fig. 8. Simplification effect of each point cloud model. (a) Bunny; (b) chair; (c) TW

下载图片 查看所有图片

6.3 简化速度

整个点云简化过程分为边缘点检测、聚类曲率估计和点云简化三个阶段。采用整个过程的平均时间作为简化速度的评判标准,评估数据如表1所示,E表示边缘点检测,C表示聚类曲率估计,S表示点云简化。由表1可知,在不同简化率下,当簇数不变时,边缘点检测和聚类曲率估计阶段的时间是一致的,其简化速度的差异主要集中在点云简化阶段(基本上是一个数量级的)。在相同简化率下,随着簇数的增加,计算成本逐渐增加,简化时间以非线性的方式增加,简化速度也逐渐变慢。随着簇数的增加和简化率的降低,误差逐渐减少,这是可预测的,因为该算法提取出更多的特征点,使拟合曲面误差变小,但时间成本增加。

图 9. 不同简化率下的重建效果。(a)斯坦福兔子;(b)椅子;(c)兵马俑

Fig. 9. Reconstruction effect at different simplification rates. (a) Bunny; (b) chair; (c) TW

下载图片 查看所有图片

6.4 对比实验及误差分析

本文提出一种特征保留的点云简化方法,通过结合EM聚类和曲率估计方法提取局部细节特征。因此,对比实验中,在同等软硬件环境下同曲率自适应简化算法(CAA)和聚类简化方法(CSM)进行对比分析,如表2所示。在相同点云模型和简化比下,不同方法的简化时间、最大误差和平均误差的差异进一步验证了本文所提算法的有效性。

7 结论

为了满足数字博物馆存储数字化和资源共享等功能的要求,需要精确和有效地对三维点云数据进行简化,故提出一种结合边缘点保留和EM聚类(允许生成更均匀的聚类)的三维点云简化算法。1) 通过构造八叉树搜索每个点的k近邻点,并计算每个点的法向量,以此检测边缘点并保留;2) 使用期望最大化算法根据点的局部分布对点云进行聚类,并确定高曲率的点;3) 使用边缘感知的DHD方法进行点云简化,合并上述点云并删除重复点。该方法可以在曲率较大的地方尽可能保留稠密的点云,在曲率较小的地方仅保留稀疏的点云,在保留原始模型尖锐特征的同时尽可能展示其整体轮廓。重建的简化点云基本保留了原始模型的几何特征和轮廓外貌,实现了较低的几何简化误差。

本文采用基于EM聚类的特征提取方法,因此存在聚类簇数量的估计问题。研究一种生成更均匀的聚类簇的方法和解决聚类簇数量的自动估计方法将是未来研究的重点。