1 引言

随着三维扫描技术的不断发展,三维测量技术已广泛应用于人机交互、虚拟现实(VR)与增强现实(AR) 游戏、机器人等领域。现有的三维扫描仪虽然可以轻松获取海量点云数据,但是其中包含大量冗余数据点。过密的数据不仅会制约计算机处理的速度,而且会影响曲面的光顺性^[1-3]。因此在保留模型特征的前提下,有必要对点云进行精简。

国内外许多学者对点云精简开展了相关研究。Martin等^[4]提出均匀网格法,将点云分配至均匀网格中,以网格中的中值点代替该网格内所有点,该算法未考虑特征点,易丢失模型的细节特征。Lee等^[5]提出一种非均匀网格精简法,根据均匀网格内点云法矢偏差值对网格进一步细分,将点云数据分为特征区域和非特征区域,特征区域保留更多点,非特征区域保留少数点。Chen等^[6]提出基于三角网格的点云精简算法,对扫描的点云数据进行三角化处理,利用向量加权算法对冗余的三角网络进行判定并删除。Pauly等^[7]根据曲面变化程度分别提出了顶点聚类法和迭代精简法,顶点聚类法运行速度快,但平均误差较大,迭代法误差小,但采样率不好控制。Weir等^[8]提出包围盒法,将点云模型包含在一个立体包围盒中,对包围盒不断均匀细分,取细小包围盒的中心点代替该包围盒的所有数据,该算法原理简单、易于实现,对均匀分布点云数据的精简效果良好,但包围盒的尺寸大小由人为设定,中心点由计算而非原始点云数据获得,且人为区分模型的特征和非特征区,易造成细节特征的丢失。周煜等^[9]将八叉树空间分割方法和平均曲率法相结合应用于点云精简,将点云模型压入空间六面体包围盒内,设定平均曲率阈值,根据包围盒内点云数据平均曲率大小判定是否进一步细分,当平均曲率大于阈值时,细分包围盒,反之则保留立方体内曲率值最接近平均曲率的点,并删除其他数据点。张雨禾等^[10]提出一种基于泊松分布的点云简化算法,将单位距离上法向变化作为局部检测算子,根据曲面弯曲程度选择不同的简化策略。刘迎等^[11]根据所有点云的微分信息,将点云数据分为平面数据和非平面数据,对平面数据进行边界提取,对非平面数据根据曲面弯曲度不同进行不同程度的精简。文献[ 1]提出了一种保特征的点云精简算法,根据点云的微分信息进行K均值聚类,并对点云空间进行全局聚类,将各个类划分为平坦区域和非平坦区域,针对非平坦区域再细分多个子类,最后以聚类中心点代替该类数据点。

当前主流精简算法是通过点云微分信息将点云数据划分为特征区和非特征区,或平面区和非平面区^[12-14],针对不同区域采取不同精简策略,以在实现精简的同时较好地保留模型细节特征。但对于高密度的点云可能无法准确计算其微分信息,且点云微分信息值与K邻域的K值有关,不同K值时会得到不同特征值^[15-16],继而无法准确识别特征区域,导致精简效果不佳;当模型中含有噪声时,其计算结果的正确性更低;高密度的点云数据导致算法耗时长、计算成本大。为此本文提出一种基于动态栅格划分的点云精简算法,其特点是无需计算所有点云的微分信息,就可区分特征区与非特征区,再根据不同区域采用不同的精简策略,能够在精简的同时较好地保留模型的细节特征。具体算法如下:对点云模型进行初步空间栅格划分,根据栅格内点云数据平整度判断是否进一步细分,在点云特征丰富区域,由于特征凹凸的存在使得该区域点云数据起伏较大,平整度相对较差。经动态划分后,特征丰富区域的栅格间距较小,而平坦区域的栅格间距较大,提取特征区域内点云点并保留,根据栅格尺寸大小对非特征区域进行不同程度的精简。

2 动态空间划分

如图1所示,所提出的动态空间划分法非常类似于空间八叉树剖分,划分栅格的编码完全可使用八叉树编码^[17-18],与八叉树划分不同之处在于栅格动态细分条件。本文栅格动态划分分为初划分与细划分,首先将点云模型进行初步空间栅格划分,通过随机采样一致性(RANSAC)算法计算各栅格内的平整度值,根据平整度情况决定是否应进一步细分,经多次动态细分后,点云模型的平坦区域将被划分到大栅格内,而模型的特征丰富区域则被细小栅格包围。

图 1. 动态划分示意图

Fig. 1. Schematic of dynamic division

下载图片 查看所有图片

2.1 空间栅格初划分

设点云数据X、Y、Z三个方向的边值分别为X_max、X_min、Y_max、Y_min、Z_max、Z_min,则可利用一个空间长方体盒将点云模型包围,该长方体盒的8个顶点坐标分别为(X_min,Y_min,Z_min)、(X_min,Y_max,Z_min)、(X_max,Y_min,Z_min)、(X_max,Y_max,Z_min)、(X_min,Y_min,Z_max)、(X_min,Y_max,Z_max)、(X_max,Y_min,Z_max)、(X_max,Y_max,Z_max)。

设定初始栅格划分间距D_x、D_y、D_z的大小:

Dx=2NXsdxDy=2NYsdyDz=2NZsdz,(1)

式中d_x、d_y、d_z为点云数据X、Y、Z三个方向的最小间距,其大小根据点云原始数据的密度选取。理论上将点云模型以d_x、d_y、d_z大小分割,每一个空间栅格内只包含一个点云数据,实际操作中可能无法保证每个空间栅格内只含一个点云数据,因此,为增加算法的可操作性及稳定性,将栅格细分的最小尺度均放大X_s、Y_s、Z_s倍,放大因子X_s、Y_s、Z_s取值参考文献[ 19],即2≤X_s、Y_s、Z_s≤5,本文方法中取X_s=Y_s=Z_s=3。N表示初始间隔最多可细分的次数。

点云模块被初划分为N_x×N_y×N_z个,其中N_x、N_y、N_z分别表示为

{Nx=「(Xmax-Xmin)/Dx⌉Ny=「(Ymax-Ymin)/Dy⌉Nz=「(Zmax-Zmin)/Dz⌉,(2)

式中「·⌉表示向上取整。

N次细分后,最小栅格的间距为minD_x=X_sd_x、minD_y=Y_sd_y、minD_z=Z_sd_z。

2.2 栅格动态细分

将平整度定义为栅格内点云到拟合平面的最大距离,其大小决定是否对该栅格进一步细分。由于拟合平面的状况将直接影响到栅格划分情况,为了提高算法的抗噪能力,减少粗大误差数据点对平面拟合的影响,首先利用RANSAC法剔除异常点(明显远离拟合平面的点)^[20],再利用最小二乘法进行平面拟合,这样可保证绝大多数点云参与拟合,确保平整度计算的正确性。

具体细分方法如下。

步骤1 RANSAC平面拟合。任选一个含有点云数据的初始栅格,使用RANSAC和整体最小二乘法求得该栅格内点云的拟合平面。

1) 在栅格内任选非共线的3点,可计算该三点的拟合平面方程为z=ax+by+c,其中a、b、c为待定系数。

2) 计算栅格内其他各点到上述平面的距离d_i=axi+byi+c-zi/a2+b2+1,并计算距离标准差σ=∑i=1k(di-d-)2k-1,其中 d-为平均距离。

3) 取阈值t=2σ,当d_i>t时,表示该点为异常点,反之,则认为是有效点,保留并统计个数M。

4) 一般重复以上步骤3~4次,即可以获得足够多的有效数据,至此RANSAC算法结束,随后可以利用最小二乘法对有效数据进行平面拟合,求解拟合的平面方程z=ax+by+c的参数。

步骤2 平整度计算及细分判断。计算栅格内各点到最终拟合平面的距离 dpi=(ax_i+by_i+c-z_i)/a2+b2+1,定义距离差d_s=max( dpi)-min( dpi)为平整度,记平整度逻辑函数d_flag=d_s>d_threshold,记栅格划分次数的逻辑值N_flag=n<N,根据d_flag、N_flag的逻辑值进行相应操作,具体伪代码说明如下:

if d_flag=true {∥此条件说明栅格内的数据起伏变化较大,栅格内含丰富特征信息。

if N_flag=false {∥栅格已被细分到最小间距,此栅格已经为最小栅格并对该栅格编码,但栅格内含丰富特征信息,计算细小栅格内点云数据的法向量。

} else {∥说明此栅格并非最小间距栅格,栅格内的数据起伏变化较大,含丰富特征信息,对栅格进一步细分为8个子栅格并分别编码,具体编码方式参考八叉树编码,累加细分次数n++,返回步骤 1重新计算各个子栅格平整度,再判断。}}

else{∥此条件说明栅格内数据平整度较好,栅格数据基本处于同一平面内,停止细分,并对该栅格编码。}

重复操作步骤 1、步骤 2,直到所有初始栅格全部处理完毕。

3 特征点提取保护

经动态栅格划分后,点云的特征丰富区域被划分至细小栅格内,平坦区域划分至相对较大间距的栅格内。点云精简的原则是尽可能保证模型的细微特征,为此对细小栅格内点云特征点进行提取及保护。

模型的特征点多在凹凸处,即在该区域内相邻点云的法向夹角和曲面曲率较大,文献[ 7]以曲面变化度H=λ₁/(λ₁+λ₂+λ₃)大小衡量是否为特征点,即H值越大则该点为特征点的可能性就越大。其中λ₁、λ₂、λ₃为半正定矩阵C的特征值,

C=p1-p-︙pk-p-Tp1-p-︙pk-p-,(3)

式中{p₁,p₂,…,p_k}为点p_i的邻域点集, p-=1k∑j=1kp_j为该点集的中心点。因矩阵C为三阶对称的半正定矩阵,必定存在特征值λ₁、λ₂、λ₃(λ₁≥λ₂≥λ₃>0)。e₁、e₂、e₃为矩阵C的特征向量,特征向量e₃可近似认为是该点p的法向量n_vector=e₃。

由(3)式可知k取值的大小会影响λ₁、λ₂、λ₃值,且各点之间的H值变化较小,将较大程度降低特征识别的准确性。文献[ 13]、[21]通过使用数据点p_i与所有邻域点的法向夹角和识别特征点,即定义法向夹角和 Wθij=1k∑j=1kθ_ij,其中θ_ij=arccos npi·npjnpi×npj, Wθij越大表示该点p_i为特征点的可能性就越大。由于法向量值也与k取值有关,因此单以法向夹角识别特征点的方法也存在一定缺陷。文献[ 10]、[22]验证了由于点云数据的不均匀性,使得部分非特征点与远距离点的夹角过大,导致特征点误判的现象。对此问题文献[ 10]、[22]引入两点欧氏距离的特征调节系数,即以两点间单位距离上的法向夹角值β_ij=θij‖pi-pj‖代替夹角信息,通过降低远距离点对特征识别的贡献权重,提高识别的准确性和稳定性。

为提高特征点识别的准确性和稳定性,在综合上述文献思路的基础上,本文方法将曲面变化度和邻域法向量夹角信息相结合共同确定特征点,即

Ci=λHHi+(1-λH)1k∑j=1kθij×g(‖pi-pj‖),(4)

式中g(‖p_i-p_j‖)为高斯函数用于调节各点法矢量夹角的权值,其本质思想与文献[ 10]、[22]的欧氏距离权重类似,即距离点p_i越远的点对法向量夹角和的贡献就越小。C_i大于设定阈值C_Threshold表示该点p_i为特征点,反之为非特征点。使用本文方法对bunny模型进行特征点提取,具体结果如图2所示。

图 2. (a)原始点云和(b)特征点

Fig. 2. (a) Original point cloud and (b) feature points

下载图片 查看所有图片

4 精简策略

本文方法点云精简基本思想为针对不同间距栅格区域采取不同的精简策略。在特征丰富区域内强制保留特征点并随机采样其余点,平坦区域根据栅格间距的大小和栅格内点云的密度信息选择不同采样率的随机采样法。由于经动态栅格划分后,细小栅格的位置区域是丰富特征区域,因此特征点的识别无需计算所有点的微分信息即可确定,这样能快速完成点云精简同时保留模型的细微特征。

假设原始点云点数为N₀,提取的特征点数为N_fp,细分次数N,各间距栅格数目 Ngridi,点云精简率为δ,精简后点云数为(1-δ)N₀,则各间距栅格随机采样率δ_i满足

∑i=1Nδi=(1-δ)N0-εNfpNgridi。(5)

间距大的栅格区域由于平坦性较好,可采用较低的随机采样率,而间距相对较小的栅格采用较大采样率采样方法,其中ε∈(0,1]。当精简后点数小于或约等于检测到的特征点数时,根据精简率值保留部分特征点数εN_fp。

5 算法的分析与改进

考虑到动态划分的随机性,可能会出现模型的边界正好处于两个大栅格之间的情况,此时按上述计算方法将无法检测到该特征边界。为避免此种错误,对所有大栅格取距离栅格中心点最近点云的法向量作为栅格法向量,计算相邻栅格法向量夹角,如果大于一定阈值θ'_difference,则表示该栅格的相接处存在特征边界,随后在该区域内对特征点进行精确定位。

6 实验实例及分析

6.1 有效性验证

为了验证本文方法的有效性,选取bunny模型,采用C++语言在VS2010平台上实现精简,并且调用OpenGL库函数显示点云,具体精简结果如图3所示。

图 3. 精简结果。(a)精简35.98%;(b)精简65.23%;(c)精简78.12%;(d)精简85.41%

Fig. 3. Simplification results. (a) Reduced by 35.98%; (b) reduced by 65.23%; (c) reduced by 78.12%; (d) reduced by 85.41%

下载图片 查看所有图片

原始bunny模型有34863个数据点,检测特征点数为5642,为了避免在高精简率情况下非特征区域出现过精简现象(精简后出现空白区域),在精简率为65.23%和78.12%时,设参数ε=0.8,当精简率为85.41%时,设参数ε=0.6,对特征信息点随机采样,从精简效果图可看出,本文方法在高精简率情况下能够较好地保留模型的特征信息。

为了进一步验证本文方法精简效率,将本文方法与Geomagic内经典的点云精简方法进行比较,其中包括随机采样法、栅格法和曲率精简法,分别对shell模型进行不同程度的精简,如图4~7所示。

图 4. 随机采样法精简结果。(a)精简50%;(b)精简75%;(c)精简87.5%;(d)精简93.75%

Fig. 4. Simplification results of random sampling method. (a) Reduced by 50%; (b) reduced by 75%; (c) reduced by 87.5%; (d) reduced by 93.75%

下载图片 查看所有图片

图 5. 栅格法精简结果。(a)精简51.2%;(b)精简75.1%;(c)精简87.43%;(d)精简93.66%

Fig. 5. Simplification results of grid method. (a) Reduced by 51.2%; (b) reduced by 75.1%; (c) reduced by 87.43%; (d) reduced by 93.66%

下载图片 查看所有图片

图 6. 曲率精简法结果。(a)精简50%;(b)精简75%;(c)精简87.5%;(d)精简93.75%

Fig. 6. Simplification results of curvature method. (a) Reduced by 50%; (b) reduced by 75%; (c) reduced by 87.5%; (d) reduced by 93.75%

下载图片 查看所有图片

图 7. 本文方法精简结果。(a)精简51.5%;(b)精简75.08%;(c)精简87.53%;(d)精简93.73%

Fig. 7. Simplification results of proposed method. (a) Reduced by 51.5%; (b) reduced by 75.08%; (c) reduced by 87.53%; (d) reduced by 93.73%

下载图片 查看所有图片

由精简结果可知,随着精简率的提高,精简后模型的细微特征变模糊。特别是随机采样法和栅格法,当精简到87%左右时,模型的细节特征丢失严重,曲率精简法在高精简率情况下,平坦区域出现过精简现象,如图6(c)和图6(d)所示。相比其他三种精简方法,本文方法在精简率87.53%和93.73%时也能较好地保持模型细节特征,且避免孔洞区域的出现。

为更客观地评估精简质量,引入文献[ 23]的评估方法,求精简后shell模型与原始模型的最大误差Δ_max(S,S^*)=maxp∈Sd(p,S^*)和平均误差Δ_ave(S,S^*)=1S∑p∈Sd(p,S^*),其结果如图8所示,其中d(p,S^*)表示原始模型曲面S上点p到去噪点云曲面S^*上投影点p^*的欧氏距离。

图 8. 精简误差比较。(a)最大误差;(b)平均误差

Fig. 8. Simplified error comparison. (a) Maximum error; (b) average error

下载图片 查看所有图片

最大误差和平均误差随着精简率的提高而增大,本文方法精简后产生的最大误差相对平缓。在采样率为1∶16(精简率93.75%)时,本文方法产生的最大误差为1.502 mm,远小于其他三种方法。

6.2 算法性能分析

为了分析算法的抗噪声能力,在bunny模型的基础上分别添加了不同强度的高斯噪声并精简,计算精简75%后模型和原始模型的偏差值,以该偏差值作为衡量算法抗噪声干扰能力的指标,如表1所示。

当模型添加30 dB高斯噪声时,随机采样法、栅格法和曲率精简法精简后模型与原模型的偏差急剧增加。由于随机采样和栅格法本身就不具备抗噪声能力,当噪声强度较大时,精简后所留的点有可能是噪声点,因而偏差增加;曲率精简法由于噪声强度的增加,曲率计算的准确性降低,继而使得模型的偏差增大。相比其他三种方法,本文方法产生的模型偏差相对较平缓,具有较好的抗干扰能力,在35 dB高斯噪声的影响下,本文方法产生的平均偏差为1.82×10^-4 mm,分别为随机采样法和栅格法的40%及曲率精简法的50%。

7 结论

提出了一种动态栅格划分的点云精简方法,通过以栅格内部点云数据平整度为细分条件实行动态细分,将平坦区域压入间距较大的栅格内,特征丰富区域划分为细小栅格,查询细小栅格位置确定特征区域。针对特征点的提取,提出了一种综合曲面变化度和法向量夹角的点云特征提取方法,为提高特征点提取的准确性,引入了点云距离的高斯调节函数以降低远距离点对判断的权重,最后分析了精简的策略和算法的不足及改进。通过bunny模型的精简实验验证本文方法的有效性,通过与Geomagic内集成的三种经典的点云精简方法对比可见,本文方法能在精简率93.73%时较好地保持模型细节特征并且避免孔洞区域的出现,精简后模型的最大偏差为1.502 mm,远小于其他三种方法。对比不同强度的噪声模型的精简效果可见,本文方法精简后的模型偏差相对较平缓。当高斯噪声为35 dB时,本文方法产生的平均偏差为1.82×10^-4 mm,分别为随机采样法和栅格法的40%、曲率精简法的50%,表明本文方法的抗噪能力和稳健性优于其他三种方法。