基于迭代最小二乘的点云法向量估计方法

马学磊,薛河儒

(内蒙古农业大学,内蒙古 呼和浩特 010018)

1 引言

点云法向量的计算是点云数据处理中的重要环节。许多曲面重建算法需要精确的法向量信息才能生成高质量曲面。法向量估计可以应用在很多场合,如点云分割、点云精简、模型建模、特征检测和提取等。点云法向量可通过分析邻域点集的几何结构近似获取。目前,常用的点云法向量估计方法有两种。

第一种方法利用待估计点的邻域点集计算法向量。Hoppe等利用局部邻域拟合平面的方法估计法向量[1]。这种方法计算速度快,而且对光滑曲面性能较好。然而对于分片光滑的曲面,在细节特征处法向量估计结果不准确。Guennebaud等和Cazals等在回归中使用二次曲面或球体代替平面。然而球体和二次曲面仍然是光滑曲面,这些方法在尖锐特征处估计效果不理想[2,3]。Pauly等在平面拟合过程中对邻域点赋予高斯权,以达到削弱某些点在平面拟合中的作用[4]。Niloy等提出一种方法自适应改变邻域大小,然而该方法对位于特征线上或接近特征线上的点会产生各向异性的邻域集[5]。上述方法均在以待估计点为中心的邻域集上做回归以实现法向量估计。当待估计点位于边缘或尖锐特征区域,以上方法的法向量估计结果不准确。近年来,许多学者提出了特征保持的法向量估计方法。Huang等人提出了一种将点云重采样和法向量估计相结合的方法[6]。它能够对具有噪声和异常值的点云模型准确地估计法向量。然而,这种方法的输出为新的合并点云,原始点云的法向量并没有被计算。通过基于核密度估计的方法最大化目标函数,Li等人减少了位于不同表面上的邻域点的影响[7]。它仅对采样均匀的点云数据会产生较好的法线估计结果。

第二种方法利用Voronoi图计算点云法向量。Amenta 等提出在每个点的Voronoi晶格中,利用每个点与离其最远的Voronoi顶点的连线方向估计点云法向量[8]。QuYang等对每个点找到其局部Voronoi邻域,通过对这些邻域点构建二次曲面估计点云法向量[9]。以上两种方法的缺点是当点云噪声率较高时,法向量估计结果不准确。Alliez等人提出了一种更加稳定的法向量估计方法,该方法结合了PCA和Voronoi图的优点[10]。但是这种方法对具有尖锐特征的点云数据没有效果。

本文提出一种新的点云法向量估计方法,可对具有边缘或细节特征的点云法向量进行有效估计,且在外点噪声存在的情况下也能获得较好的效果。首先,使用主成分分析(Principal Component Analysis, PCA)估计点云的初始法向量;其次,针对PCA在边缘和细节特征处估计结果不准确的问题,使用层次聚类和迭代最小二乘法对PCA估计的法向量进行调整。在迭代过程中引入三个权函数控制邻域点对平面拟合的作用,各权函数带宽利用特征系数自适应计算获得。

2 主成分分析计算点云法向量

利用局部平面拟合的方法可计算点云法向量。对于点云中的每个点pi,获取pi的k近邻点集P,本文使用KD树存储点云数据以提高搜索的效率。通过近邻点集计算一个最小二乘意义上的局部平面,该问题可表示为

(1)

(2)

Cp为3×3实对称矩阵。对协方差矩阵Cp进行奇异值分解(SVD)得到特征值λi(i=0,1,2),其中λ2≥λ1≥λ0。特征值λ2,λ1,λ0对应的特征向量分别为v2,v1,v0,最小特征值λ0的特征向量v0即为平面的法向量。

3 迭代最小二乘的平面拟合

(3)

其中ε是一个对称、正定函数,在零处有唯一的最小值。平面拟合问题如下

(4)

(5)

(6)

图1是使用迭代最小二乘拟合平面的结果图。从图1可以看出,与标准最小二乘相比,迭代最小二乘受外点影响很小。算法1是具体的迭代最小二乘平面拟合算法。

图1 迭代最小二乘的平面拟合结果图

算法1:迭代最小二乘平面拟合算法

输入:P待检测平面的点,It最大迭代次数,γ迭代终止阈值。

初始化:nold←Φ;∥ 上一次计算的法线

对C进行奇异值分解得到特征值λ0≤λ1≤λ2,对应特征向量v0,v1,v2;

n=v0;∥v0为平面的法向量

fork=1toIt

nold=n

计算权函数w

对C进行奇异值分解,最小特征值对应的特征向量为平面法向量n;

convg=max(abs(nold-n)/abs(nold));

if convg<γ then

迭代终止;

end if

end for

4 使用迭代最小二乘计算点云法向量

4.1 点云特征系数

点云高斯映射是将非平坦聚类细分为多个子类,使特征点和平坦区域的点形成不同的聚类。设S=S(u,v)是三维空间一正则曲面,该曲面的法向量方向一致,将点q∈S处曲面的单位法矢量的始端平移到单位球S2={(x,y,z)⊂R3|x2+y2+z2=1}的中心,使曲面上的点与球面上的点相对应,这种对应被称为曲面S上的高斯隐射G:S→S2。曲面S的高斯隐射的像记为G(S),称G(S)为S的高斯图,单位球S2称为高斯球。

高斯映射可以把几何空间中曲面的法向量映射到单位球上。曲面法向量在平滑区域内变化平缓;在细节或尖锐特征区域内,法向量会发生显著变化。因此,利用法向量的高斯映射可以确定曲面的细节和尖锐特征区域。曲面平滑区域内的点及其相邻点的高斯映射在单位球上形成一个聚类簇;而细节或尖锐特征区域的点经高斯隐射后则可能产生多个聚类簇,每个聚类簇分别对应曲面上的不同的区域。平面和柱形曲面的高斯映射分别对应高斯球上一个点和圆。

曲面点云的高斯映射如图2所示。图2(a)为两个相交平面,其中绿色点是位于两相交平面边缘上的点,红色点和蓝色点是绿色点的邻域点,分别位于不同平面上。图2(b)为图2(a)中红色点和蓝色点的法向量经高斯映射后的结果。绿色点位于两相交平面的边缘处,所以它的邻域点经高斯隐射后在单位球上会形成两个聚类簇(对应图2(b)单位球上的蓝色和红色点)。

图2 曲面点云的高斯映射结果图

曲面点云的法向量经高斯隐射后,单位球上两点a和b的测地距离为法向量na和nb之间的夹角,即

d(a,b)=arccos(na,nb)

(7)

测地距离d描述了两个法向量的相似程度。利用层次聚类算法对单位球上的点进行聚类。层次聚类的基本思想是:首先把每个点当作一类;其次逐次合并两个类形成一个大类,至到任意两个类之间的距离大于给定的阈值。使用平均距离D度量两类A和B之间的距离

(8)

使用特征系数fi表示点pi是否为特征点。若fi=0,则表示pi为非特征点;若fi=1,则表示pi为特征点;若0

点pi的特征系数fi的计算步骤为:

1)首先找到点pi的邻域点集合,对各邻域点利用PCA计算的法向量进行高斯映射。

2)对高斯球上的点进行层次聚类,并统计聚类簇总数、主簇和非主簇。对每个聚类簇,如果该簇的大小超过τ,则该簇为主簇。τ为邻域集大小与聚类簇总数的比值。

3)计算特征系数。如果层次聚类后聚类簇个数为1,则特征系数fi=0;若主簇个数大于1,则特征系数fi=1;若主簇个数为1,则特征系数按式(9)计算:

fi=min([(2*p2)/p1,1])

(9)

其中p1和p2分别为主簇和非主簇内各点到离其最近的六个点的距离均值的平方和。

4.2 法向量的调整

(10)

5 试验与分析

本文对迭代最小二乘估计点云法向量进行了验证。图3为不同方法对两种类型点云数据的法向量估计侧面结果图。图3(a)是未添加噪声的点云数据,对点云数据图3(a)添加噪声如图3(b)所示。图3(c)、 (d)和(e)分别为主成分分析(PCA)、随机采样一致性算法(RANSAC)和本文方法对未添加噪声的点云数据(图3(a))的法向量估计结果。图3(f)、(g)和(h)是对添加噪声的点云数据(图3(b))使用主成分分析(PCA)、RANSAC和本文方法的法向量估计结果。

图3 不同方法对两种类型点云数据的法向量估计侧面结果图

从图3(e)和(h)可以看出,对未添加噪声和添加噪声的两种点云数据,本文方法比主成分分析(PCA) 和RANSAC对点云法向量的估计更加准确。在图3(e)和(h)中,边缘处点的法向量朝向基本一致,其它点处的法向量也更加紧凑,而且本文方法在边缘处的估计效果基本不受噪声影响。从图3(c)、(d)、(f)、(g)中可以看出,在点云边缘处,主成分分析法(PCA) 和RANSAC对法向量的估计基本失效,PCA的估计结果更易受到噪声的干扰。

为了定量评价算法性能,使用根均方误差(RMS),定义如下

(11)

(12)

表1和表2 分别是不同方法对未添加噪声和添加噪声点云法向量估计的NBP和RMS。可以看出,对于两种不同类型的点云数据,利用本文方法比主成分分析法(PCA)和随机采样一致性算法(RANSAC)对法向量的估计更加接近真实法向量,而且估计错误点数更少。

表1 各种方法对法向量估计的NBP和RMS (不含异常值)

表2 各种方法对法向量估计的NBP和 RMS (含异常值)

6 结束语

本文提出一种基于迭代最小二乘法的法向量估计方法。通过迭代最小二乘拟合局部平面的方法对点云初始法向量进行调整,进而获得最终的法向量估计结果。本文方法克服了传统方法在细节特征区域法向量估计结果不准确的缺点,同时不受外点噪声的影响。使用本文方法得到的比较准确的法向量估计结果,对后续三维点云模型的点云分割、特征检测和提取、曲面重建等具有重要意义。