尖锐特征曲面散乱点云法向估计

袁小翠,吴禄慎,陈华伟

(1.南昌工程学院 江西省精密驱动与控制重点实验室，江西 南昌 330099;2.南昌大学 机电工程学院，江西 南昌 330031)



尖锐特征曲面散乱点云法向估计

袁小翠1*,吴禄慎2,陈华伟2

(1.南昌工程学院 江西省精密驱动与控制重点实验室，江西 南昌 330099;2.南昌大学 机电工程学院，江西 南昌 330031)

针对现有算法对尖锐特征曲面点云法矢估计不准确，点云处理时容易丢失曲面细节特征等问题，提出一 种尖锐特征曲面散乱点云法向估计法。该方法用主成分分析法粗估计点云法向；然后，根据各邻域点的空间欧氏距离和法向距离对各邻域法向加权，用加权邻域法向之和来更新当前点的法向；最后，测试估计法向与标准法向的误差，评价估计法矢的准确性，并且将估计的法向应用到点云数据处理中来比较特征保留效果。实验结果表明：本文方法能够准确地估计尖锐特征曲面的法向，最小误差接近0。另外，该方法对噪声有较好的鲁棒性，点云处理时能保留曲面的尖锐特征。相比于其他特征曲面法向估计法，所提出的方法估计的法向误差更小、速度更快、耗时更少。

散乱点云；尖锐特征；法向估计；逆向工程；特征保留；主成分分析

1　引　言

法向量是点云的重要属性之一，点云法向量的有效估计是逆向工程中点云数据处理的基础，除了精确和高质量的点云绘制方法主要依赖于法向量，此外，许多有效的点云处理方法也都需要准确的法向量作为输入，如点云去噪[1]、分割[2]、数据精简[3]、曲面重建等[4]。对尖锐特征面，若不能准确估计特征区域(两个或者多个曲面交界的过渡区域)的法向量，则点云处理时容易丢失曲面的细节特征，导致重建后的曲面难以恢复原始模型的几何特征。

近年来点云的法向估计受到了越来越多的关注，点云法向估计方法可分为局部邻域拟合法和Voronoi/Dalaunay方法两类。局部邻域拟合法由文献[5]首次提出，此方法也称为主成分分析(Principal Component Analysis,PCA)法向估计法。为了更准确地估计点云法向，文献[6]利用代数球拟合邻域点来估计法向；文献[7]通过对当前点的邻域赋予高斯权重，使得与当前点距离越近的邻域点对法向量作用越大，越远的点作用越小；文献[8]使用统计学中的集成技术改进PCA方法的鲁棒性。以上这些改进方法虽然一定程度上改善了点云的法向估计结果，但是PCA方法实际上是一种低通滤波的方法，特征点处的法向被平滑处理了。为了能够准确估计尖锐曲面的法向，学者们做了大量研究[9-16]，文献[14]首次提出基于Voronoi/Delaunay方法，当点云中不含噪声或者噪声较少时Voronoi图对点云法向估计效果较好。文献[15]将Voronoi单元格极点扩展到大Delaunay球，来估计点云法向。文献[16]将回归方法和Voronoi方法相结合，以获得更稳定的结果。然而以上几种基于Voronoi/Delaunay方法均不能很好地估计尖锐特征曲面的法向。

针对以上问题，本研究提出一种邻域法向迭代加权的点云法向估计法，可以在点云处理时保留曲面的尖锐特征，相比于其他特征曲面法向估计法，所提出的方法估计的法向误差更小、速度更快、耗时更少。

2　法向估计

假设点云的采样曲面处处光滑，任何点的局部邻域均可以用平面进行拟合，给定点集X={x1,x2,…,xt}，其中t为点云总数，点xi的k邻域表示为Nb(xi)，其中Nb为邻域。对任意一点xi，用其k邻域拟合的最小二乘平面表示为：

(1)

式中：n为平面Pl的法向量，须满足‖n‖2=1；d为邻域点到拟合平面的距离。式(1)可以转化为对中半正定协方差矩阵C进行特征值分解，C最小特征值的特征向量可被当作xi的法向量，这即是PCA方法。

(2)

C可以分解为3个特征向量v2，v1和v0，3个特征向量对应的特征值分别是λ2，λ1和λ0，其中λ2≥λ1≥λ0。最小特征值对应的特征向量为平面的法向量，即n=v0。若曲面处处光滑，则PCA方法能快速准确地估计点云法向，然而当曲面包含尖锐特征时，在特征点处，由于用来拟合平面的k邻域点位于多个曲面上，则拟合平面被平滑。

文献[7]通过对当前点的邻域赋予高斯权重，使距当前点越近的邻域点对拟合平面的作用越大，距离越远作用越小，并将式(1)改写为：

(3)

文献[9]在文献[7]的基础上叠加了残差因子，又将式(3)改进为：

Pl(n,d)=argmin∑ρ(d+(x-xi)Tn)ωd(xi),

(4)

(5)

(6)

式中u=∑ωd(xi)ωr(ri)(x-xi)/∑ωd(xi)ωr(ri)。

式(5)是通过对邻域点迭代加权来逐步减小不在同一曲面的邻域点对拟合平面的作用来求解法向的，相对于文献[7]，式(5)得到的法向在特征点处改善明显。文献[10]提出了自适应特征保留法向估计法，其主要在式(5)的基础上加入了法向偏差，当前点的法向与邻域点的偏差越大，该邻域点对当前点拟合的平面作用越小，其改进为：

(7)

式(7)中有σd，σr和σn3个带宽，但是都无须用户指定，均能随着曲面的变化自适应更新。

文献[9-10]估计法向都是通过对邻域点迭代加权来拟合平面的，减小不在同一连续曲面邻域点的作用，使得最终用来拟合平面的邻域具有各项异性。然而每一次迭代的过程中都需拟合平面、求解协方差的特征向量、邻域点到拟合平面的距离，迭代次数越大计算越耗时。文献[10]引入了法向偏差作为拟合平面的权重因子，邻域点的法向与当前点的法向夹角越小，该邻域点对当前点的作用越大，反之越小。当邻域点的法向夹角大于给定带宽时，该邻域点被当做外点丢弃。在平坦区域，一般方法都能够准确估计出点云法向，只是在特征区域估计的法矢被平滑，此时可以用平坦区域获得的准确法矢来修正特征点的法矢。对于特征区域的点，其邻域位于一个或者多个曲面上，同一曲面上的法矢相似，不同曲面上点的法矢偏差较大，而且曲面越尖锐，不同曲面上对应点的法向夹角越大；因此，这里采用邻域中准确法矢迭代加权对不准确法矢进行修正。与当前点法矢夹角越大的邻域点对当前点的法矢修正作用越小，而且距离越远的点作用也越小。经过多次迭代，逐步修正不正确的法矢，本文法矢估计法可表示为：

(8)

式中σd和σn分别为距离带宽和法向偏差带宽。

3　实验结果分析

采用C++语言在Microsoft visual studio 2008上实现法向估计算法，并且调用OpenGL库函数显示点云，实验所用的计算机配置为Intel Core 2.30 GHz CPU，1.19 GB内存。为了验证算法的有效性，测试估计法矢与标准法矢偏差(估计误差)，将估计的法矢作为点云处理的输入来验证特征保留效果。将本文方法与文献[5]的PCA方法、文献[10]的自适应迭代邻域拟合法、文献[11]的鲁棒法向估计法进行定量和定性比较。

3.1参数选择

本文算法有5个参数，分别是PCA方法估计初始化法向时的k邻域、式(8)中k1邻域、距离带宽σd、法向偏差带宽σn和迭代次数t。文献[5]中用PCA方法估计法向时拟合平面的k邻域大小取8～64比较合适，本文取k=16。在平坦区域，PCA方法能准确估计点云法向，在特征点处估计的法向被平滑。若点xi是特征点，点xi的k1邻域中有许多点处于平坦区域，则用平坦区域的点来修正特征区域点的法矢，因此，式(8)中用来加权的k1邻域应该大于PCA估计法向的k邻域，这里取k1=3k=48。σd表示点xi到其邻域的欧氏距离带宽，σd值越大包括的邻域点越多，为了确保当前点的邻域都在带宽范围内，σd取点xi到邻域点的最大距离。σn表示法向偏差带宽，带宽内的邻域对当前法向的修正作用大，带宽外的邻域对xi法向的修正作用小，法向偏差越大，作用越小。当曲面夹角比较大(或者过渡面比较平缓)时，σn应该取较小值；当曲面的夹角比较尖锐时σn应该取较大值。实验发现：对于过渡边比较平滑和尖锐的曲面，带宽分别取0.2和0.3时结果比较理想。实际物体的曲面包括不同尖锐程度的曲面，σn取全局值难以满足要求，然而设定自适应带宽较复杂，计算耗时。对过渡边比较平滑的曲面，若σn取值偏大，或者尖锐曲面的σn取值偏小，则须增加迭代次数来补偿过大或过小带宽带来的不足。当曲面比较光滑时，t取较小值；曲面越尖锐，t越大，实验发现对过渡边比较平缓和尖锐的模型，迭代次数分别为5，10时，模型为阈值法向基本调整到最佳状态。当t增大时，法向几乎不再变化，为了提高计算效率，折衷带宽值和迭代次数t将带宽设为全局值σn=0.25，最大迭代次数为t=15。

3.2误差分析

在设备研发方面，2019年我们的目标是日日顺物流有一样设备在全国是最大的，在智能化设备有一个领域是全国最大的。

为了测试估计法矢的准确性，计算估计法向与标准法向的偏差，采用法向偏差均方根(Root Mean Square，RMS)来表示估计误差，RMS表示估计法向与标准(参考)法向的平均偏移量，RMS值越大，估计法向的误差越大，文献[9-12]均采用RMS来测量估计法向的误差。

式中：nref，nest分别为单位化的参考法向和估计法向；|S|为曲面S的点云数量；τ为阈值，当参考法向与估计法向的点积大于阈值τ时，认为该点的法向与标准法向偏差为90°。文献[9-12]法向夹角阈值均取为10°，即τ=0.984 8。

因为复杂曲面难以得到标准法向，所以用Matlab合成规则点云模型(不含噪声)，根据点云坐标计算其标准法向。合成模型分别是立方体和圆柱体，两模型的点云数量分别为6 147和28 832。立方体模型包含尖锐特征点，如两平面交叉的边界特征点和三平面交叉的角点。圆柱体包含的特征点比较平缓，其特征点是上下两平面和圆柱面的交叉点。为了测试估计法向对噪声的敏感性，对理想模型添加不同程度的噪声，计算理想模型及各噪声模型的RMS。当估计的法向与标准法向偏差大于10°时称该点为坏点，4种方法比较结果如图1(彩图见期刊电子版)和图2(彩图见期刊电子版)所示，图中：SSNR为信噪比；红色点表示坏点。表1所示为两个模型的RSM和坏点数量(Bad Point Number, BPN)。

(a) PCA　(b) 文献[11]　 (c) 文献[10]　 (d) 本文方法(a) PCA　(b)Reference　 (c)Reference　 (d) Proposed method　　　[11]　　　[10]图1　立方体模型的法向估计结果Fig.1　Normal estimation results of cube point cloud model

可知，对于立方体模型，当模型不含噪声时，文献[10]和本文方法没有坏点，文献[11]在特征点处有少数的坏点。当SSNR=40 dB时，本文方法的BPN为0，RSM接近0，文献[11]和文献[10]均有少数坏点，PCA方法的坏点主要是特征点及其周围的点。当数据被噪声严重干扰时(SSNR=30 dB)，如图1中第三行所示，立方体已经严重失真，角点和边界点模糊。PCA方法估计法向有一半是坏点(BPN为3 064，详见表1)，其次是文献[11]，本文方法得到BPN和RSM最小，其坏点主要分布在角点周围。对圆柱体模型，本文方法得到的BPN和RSM比其他方法的都小，几种方法结果比较见图2和表1。

为了清楚比较各坏点的法向与标准法向的偏差，统计立方体和圆柱体模型在SSNR=30 dB时不同偏差角度下坏点数量，称之为法向偏差直方图，结果如图3(彩图见期刊电子版)所示。对于立方体模型，在法向偏差为[10°，20°]区域，PCA，文献[11]，文献[10]及本文方法分别有65.5%，75.4%，83%和97.7%的坏点，在法向偏差为[20°, 30°]区域，分别有14.5%，8.4%，9.2%和2.3%的坏点。同样，对圆柱体模型，在角度偏差为[10°, 20°]区域，本文方法的坏点数的比例最大。本文方法估计的法向其坏点主要分布在低法向偏差区域，即估计得到的法向更准确。

(a) PCA　(b)文献[11]　 (c) 文献[10]　 (d) 本文(a) PCA　(b)Reference　(c)Reference　(d) Proposed method　　　[11]　　　[10]图2　圆柱体模型法向估计结果Fig.2　Normal estimation results of cylinder point cloud model

模型SSNR/dB　　　PCA　　　文献[11]　　文献[10]　　　本文　BPNRMSBPNRMSBPNRMSBPNRMS立方体无噪10980．6642230．11660000．03104011120．6702451510270．123500．05003030641．111911670．69007390．55164550．4376圆柱体无噪2490．3321340．0654140．057840．03094014310．36151780．1329760．093670．04243084600．872340630．604932190．541624050．4730

(a) 立方体模型坏点直方图 　　　　　(b)圆柱体模型坏点直方图(a) Histogram of BPN for cube model 　　　 (b) Histogram of BPN for cylinder model图3　模型中坏点直方图分布Fig.3　Histogram of BPN for point cloud models

3.3法矢可视化和应用

3.3.1法矢可视化

对两个经典模型(八面体和Fandisk模型，其中八面体模型添加了SSNR=40 dB的噪声)进行法向可视化，如图4～6(彩图见期刊电子版)所示，红色线表示一致的法矢方向，八面体模型的特征点包括多个平面交接的角点和两个平面交接的过渡边界点，角点处的曲面比较尖锐，过渡边则相对平缓。Fandisk模型由平面、球面、圆柱面等基面组成，其过渡边有尖锐边和平缓边等。为了更清楚的显示模型，将模型三角化加光照显示，如图4(a)和5(a)所示。准确的法向应该垂直于点云所在的局部曲面，两个不同曲面上的点云法向夹角为两曲面的夹角，局部区域同一曲面上的法向接近平行，在两个曲面交界的边界处有清晰的边界。图4(d)和(e)以及图5(d)和(e)的边界清晰，图4(c)和图5(c)有部分特征点的法向误差大，如角点的法向不垂直于其所在的曲面。图6(a)为图5(a)红色矩形框的区域的放大图，图6(b)～(e)为图6(a)曲面对应方法估算的法向量。

(a)三角化曲面　(b)PCA　　(c) 文献[11]　 　　　 (d)文献[10]　 　　　 (e)本文　　(a) Wraped surface 　　 (b) PCA 　　 　 (c) Reference [11]　 　 (d) Reference[10]　 (e) Proposed method图4　八面体模型法矢可视化Fig.4　Normal vector visualization of Octahedron point model

(a)三角化曲面　(b)PCA　　　(c) 文献[11]　 　　　 (d)文献[10]　 　 (e)本文　　(a) Wraped surface 　　(b) PCA 　　 　(c) Reference [11]　 　 (d) Reference[10]　 (e) Proposed method图5　Fandisk模型法矢可视化Fig.5　Normal vector visualization of Fandisk point model

(a)放大后的三角化曲面　　(b)PCA　　 (c) 文献[11]　 　　 (d)文献[10]　　 (e)本文　　　　(a) Amplified surface　　　(b) PCA 　　 (c) Reference [11]　　 (d) Reference[10]　　 (e) Proposed method 图6　局部区域放大后的法矢Fig.6　Amplification and visualization of normal vector

将估计的法向应用到点云处理中，对噪声模型进行去噪。双边滤波去噪法能在去噪的同时保留曲面的几何特征，具有特征保留的性质，但是获得准确的法向是双边滤波特征保留的前提，若法向在特征点处被平滑，则双波滤波后尖锐特征也被平滑，文献[18]对经典双边滤波算法进行了改进，其特征保留效果更优。为了消除其他因素的影响，双边滤波算法中除了法向不一样外，其他参数值均相同。基于不同法向的轴承模型去噪结果如图7所示。PCA估计的法向在特征点被平滑，从而7(c)中的边界被平滑，文献[10]和本文方法估计的法向在特征点处比较准确，去噪后可保留清晰的边界。

图7　基于不同法向的轴承模型去噪Fig.7　Denoising results of bearing model with different normals

3.4耗时分析

显然，PCA方法估计法矢最简单，而且其只对每个点的邻域拟合一次平面、从而耗时最小。文献[11]方法需从一个大邻域(邻域大小k=500)中采用随机采样一致法寻找一个最优切平面，计算非常耗时。文献[10]方法在每一次迭代中都需拟合一个平面，计算邻域点到拟合平面的距离，采用层次聚类法估计特征系数。本文方法只在迭代初始化时拟合平面，其后的每一次迭代只对邻域点加权，且迭代次数与文献[10]相同，相比文献[10]和文献[11]，本文方法更简单，耗时更小。4种方法对以上模型的法向估计耗时如表2所示，可见PCA方法最快，本文方法次之，文献[11]最慢。

表2　4种方法耗时比较

4　结　论

本文对尖锐特征曲面的法向估计进行了研究。用PCA方法粗估计曲面法向。考虑到PCA方法估计的法矢在平坦区域比较准确，在特征区域被平滑，又采用邻域点中的正确法矢来修正特征区域不准确的法矢。因此，对邻域法向加权，其权重包括距离权重和法向偏差权重，使距当前点越远的邻域点对当前法向修正作用越小，反之越大，与当前点法向偏差越小的邻域点的作用越大，修正的法矢是邻域法矢的迭代加权和。实验结果表明： 本文方法能够准确估计尖锐特征曲面的法向，最小误差等于零，估计的法向应用到点云处理时能很好地保留曲面几何特征，相比于其他特征曲面法向估计法，本文方法简单且快速。然而本文方法的法向偏差带宽是根据实验经验设定的，将来需进一步研究自适应带宽以减少迭代次数来提高算法的执行效率。

[1]吴禄慎,史皓良,陈华伟. 基于特征信息分类的三维点数据去噪[J]. 光学 精密工程, 2016, 24(6): 1465-1473.

WU L SH, SHI H L, CHEN H W. Denoising of three-dimensional point cloud based on classfication of feature information [J].Opt.PrecisionEng., 2016, 24(6): 1465-1473.(in Chinese)

[2]WANG Y H, HAO W, NING X J,etal.. Automatic segmentation of urban point clouds based on the Gaussian map[J].PhotogrammetricRecord, 2013, 28(144): 342-361.

[3]袁小翠, 吴禄慎, 陈华伟. 特征保持点云数据精简[J]. 光学 精密工程, 2015, 23(9): 2666-2676.

YUAN X C, WU L SH, CHEN H W. Feature preserving point cloud simplification [J].Opt.PrecisionEng., 2015, 23(9): 2666-2676.(in Chinese)

[4]TANG P, HUBER D, AKINCI B,etal.. Automatic reconstruction of as-built building information models from laser-scanned point clouds: A review of related techniques [J].AutomationinConstruction, 2010, 19(7): 829-843.

[5]HOPPE H, De ROSE T, DUCHAMP T,etal.. Surface reconstruction from Unorganized Points [J].ComputerGraphics, 1992, 26(2): 71-78.

[6]GUENNEBAUD G, GROSS M,ZURICH E. Algebraic point set surfaces [J].ACMTransactionsonGraphics,2007, 26(3):23.

[7]PAULY M, GROSS M, KOBBELT L P. Efficient simplification of point-sampled surfaces [C].IEEEVisualization2002,WashingtonDC:IEEEComputerSocietyPress, 2002:163-170.

[8]YOON M, LEE Y, LEE S,etal.. Surface and normal ensembles for surface reconstruction[J].Computer-AidedDesign, 2007, 39(5): 408-420.

[9]MEDEROS B, VELHO L, FIGUIREDO L H. Robust smoothing of noisy point clouds [C].ProceedingsoftheSIAMConferenceonGeometricDesignandComputing,Seattle:SIAM, 2003: 405-416.

[10]WANG Y, FENG H Y, DELORME F,etal.. An adaptive normal estimation method for scanned point clouds with sharp features [J].Computer-AidedDesign, 2013, 45(11): 1333-1348.

[11]LI B, SCHNABEL R, KLEIN R,etal.. Robust normal estimation for point clouds with sharp features [J].Computers&Graphics, 2010, 34(2): 94-106.

[12]苏志勋, 栗志扬, 王小超. 基于法向修正及中值滤波的点云平滑[J]. 计算机辅助设计与图形学学报, 2010, 22(11): 1892-1898.

SU ZH X, LI ZH Y, WANG X CH. Denoising of point-sampled model based on norm al mollification and median filtering [J].JournalofComputer-AidedDesign&ComputerGraphics, 2010, 22(11): 1892-1898.(in Chinese)

[13]ZHANG J, CAO J, LIU X,etal.. Point cloud normal estimation via low-rank subspace clustering[J].Computers&Graphics, 2013, 37(6): 697-706.

[14]AMENTA N, BERN M. Surface reconstruction by Voronoi filtering [J].Discrete&ComputationalGeometry, 1999, 22(4): 481-504.

[15]DEY T K, GOSWAMI S. Provable surface reconstruction from noisy samples [C].IEEESymposiumonParallelandLarge-DataVisualizationandGraphics,Piscataway:IEEE, 2001: 19-27.

[16]ALLIEZ P, COHEN-STEINER D, TONG Y,etal.. Voronoi-based variational reconstruction of unoriented point sets [C].Proceedingsofthe5thEurographicsSymposiumonGeometryprocessing,Geneve:Wiley, 2007: 39-48.

[17]HUBER P J.RobustStatistics[M]. Berlin: Springer, 2011.

[18]ZHENG Y, FU H, AU O C,etal.. Bilateral normal filtering for mesh denoising[J].IEEETransactionsonVisualizationandComputerGraphics, 2011, 17(10): 1521-1530.

袁小翠(1988-)，女，江西抚州人，博士，讲师，主要研究方向为图像处理与逆向工程。E-mail：yuanxc2012@163.com

吴禄慎(1953-)，男，江西乐平人，硕士，教授，博士生导师，1978年于北京航空航天大学获得学士学位，1990年于清华大学获得硕士学位，主要从事面外“moire”法、三维光学图像测量与逆向工程的研究。E-mail：wulushen@163.com

(版权所有未经许可不得转载)

Normal estimation of scattered point cloud with sharp feature

YUAN Xiao-cui1*,WU Lu-shen2,CHEN Hua-wei2

(1.Jiangxi Province Key Laboratory of Precision Drive & Control,NanchangInstituteofTechnology,Nanchang330099,China;2.SchoolofMechanicalandElectricalEngineering,NanchangUniversity,Nanchang330031,China)*Correspondingauthor,E-mail:yuanxc2012@163.com

A novel method was proposed to estimate the normal for a scattered point cloud with sharp features to overcome the shortcomings that existing methods are unable to reliably estimate the normal of point cloud model and lead to the smoothed sharp features. With proposed method, the normal of point cloud was estimated with principal component analysis method. Then, different values were weighted on neighborhood normals according to spatial distance and normal distance of current points of the neighborhood, and the revised or current normals were updated by the sum of weighted neighborhood normal. Finally, the average deviation between standard normal and estimated normal was measured and the accuracy of estimated normal was evaluated. The estimated normal was applied to point cloud processing to verify the feature-preserving property. The proposed method was validated. The results demonstrate that proposed method accurately estimates the normal for data with noise and the least average deviation is close to 0. Moreover, the method has good robustness to the niose, and it keeps the original geometry well when the normal is used as input of the point cloud processing. Comparing with other sharp feature preserving normal estimation methods, the proposed method shows smaller average deviation, higher processing speeds and less computation time.

scattered point cloud；sharp feature; normal estimation; reverse engineering; feature preserving; principal component analysis

2016-06-15；

2016-07-19.

国家自然科学基金资助项目(No.51365037,No.51065021)

1004-924X(2016)10-2581-08

TP391.7

Adoi:10.3788/OPE.20162410.2581