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

近年来点云的法向估计受到了越来越多的关注，点云法向估计可以简单分为两类：基于局部邻域拟合法和基于Voronoi/Dalaunay方法。基于局部邻域拟合法首先由Hoppe等人^[5]首次提出，给定点p，通过对点p的k近邻进行最小二乘拟合得到点p的切平面，应用主成分分析法(PCA) 求解由邻域点构成的协方差矩阵，协方差矩阵最小特征值所对应的特征向量即为平面的法向量，同时也为点p的法向，此方法也称为PCA法向估计法。当点云模型处处光滑，PCA方法能比较准确地估计点云法向，当物体中包含尖锐特征时，PCA方法估计的法向在特征点处不准确。为了能够准确估计尖锐曲面的法向，学者们做了大量研究^[6-13]，Guennebaud等人^[6]利用代数球拟合邻域点来估计法向；Pauly等人^[7]通过对当前点的邻域赋予高斯权重使距当前点越近的邻域点对法向量作用越大，距离越远的点作用越小。以上这些改进的PCA方法一定程度上提高了点云法向的准确性，但是在特征点处，拟合曲面的法向量各向同性，从而估算的法向量被平滑。为了使拟合曲面的邻域为各向异性邻域点，Mederos等人^[8]通过对$k$近邻点云叠加权重因子，其权重因子包括邻域点到当前点的距离因子和邻域点到拟合平面的残差因子，用加权的邻域点拟合曲面，使不在同一曲面的邻域点对拟合曲面的作用减小；Li等人^[9]根据随机采样策略提出了鲁棒统计点云法向估计法，该方法通过在一个大邻域(邻域$k$=500) 中利用随机采样策略选择一个小邻域($k$=16) 拟合曲面，当大邻域中的点到拟合曲面的残差最小时所确定的邻域为所求邻域，该邻域拟合曲面对应的法向量为当前点的法向量，文献[9]估算的法向量能有效抑制外点和较准确的估计特征点的法向；Zhang等人^[10]采用低秩子空间聚类法将特征点的各向同性邻域映射到子空间，K均值聚类法在子空间中将各向同性邻域分割为多个各向异性子邻域，以最优子邻域拟合曲面来估算法向量。文献[9-10]均能比较准确地估算点云法向量，且对噪声有较好的鲁棒性，然而，文献[9-10]方法非常耗时，很难将这些方法实时应用。文献[11]首次提出基于Voronoi/Delaunay方法，当点云中不含噪声或者噪声较少时Voronoi图对点云法向估计效果较好。文献[12]将Voronoi单元格极点扩展到大Delaunay球来估计点云法向。文献[13]将基于曲面拟合方法和Voronoi方法相结合能获得更稳定的结果。然而以上几种基于Voronoi/Delaunay方法均不能很好地估计尖锐特征曲面的法向。

针对以上问题，提出基于高斯映射的特征曲面点云法向量估算方法，对点云的$k$邻域采用PCA方法粗略估算点云的法向量和候选特征点，将特征点及其各向同性邻域映射到高斯球，用K均值聚类法对高斯球上的数据分割成多个子集，以最优子集对应的各向异性邻域拟合曲面来精确估算特征点的法向量。

给定点集$\boldsymbol{P}=\{ {p_1}, {p_2}, \ldots, {p_n}\} $，其中$n$为点云总数，点${p_i}$的最近k邻域表示为${N_b}({p_i})$，其中${N_b}$表示邻域。对任意一点${p_i}$用其$k$邻域拟合的最小二乘平面的表示为

$ Pl\left( {n,d} \right) = \arg \min \sum\limits_{{p_i} \in {N_b}} {{{\left( {n{p_i} - d} \right)}^2}} $

(1)

式中，$\boldsymbol{n}$为平面$Pl$的法向量，$\boldsymbol{n}$须要满足${\left\| \boldsymbol{n} \right\|_2}=1$；$d$为邻域点到拟合平面的距离。式(1) 可以转化为对式(2) 中半正定协方差矩阵$\boldsymbol{C}$进行特征值分解，即

$ \mathit{\boldsymbol{C = }}{\left[ {\begin{array}{*{20}{c}} {p - {p_1}}\\ {p - {p_2}}\\ \vdots \\ {p - {p_k}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {p - {p_1}}\\ {p - {p_2}}\\ \vdots \\ {p - {p_k}} \end{array}} \right] $

(2)

式中，$\boldsymbol{C}$最小特征值的特征向量可被当作${p_i}$的法向量，这就是所谓的PCA方法。协方差矩阵$\boldsymbol{C}$定义了局部曲面的几何信息，为一个对称的半正定矩阵。$\boldsymbol{C}$可以分解为3个特征向量${\boldsymbol{v}_2}$，${\boldsymbol{v}_1}$和${\boldsymbol{v}_0}$，3个特征向量对应的特征值分别是${\lambda _2}$、${\lambda _1}$和${\lambda _0}$，其中${\lambda _0} \le {\lambda _1} \le {\lambda _2}$。最小特征值对应的特征向量为平面的法向量即$\boldsymbol{n}={\boldsymbol{v}_0}$。

最小特征值${\lambda _0}$表示点沿着曲面法矢的变化情况，因此可以利用各点的最小特征值的变化来估计高曲率区域。点$p_i$的曲面变化表示为

$ {\omega _i} = \frac{{{\lambda _0}}}{{{\lambda _0} + {\lambda _1} + {\lambda _2}}} $

(3)

当${\omega _i}=0$时，表示$p_i$为平面点，${\omega _i}$值越大，曲面越尖锐。根据${\omega _i}$值判断点$p_i$是否处于高曲率区域，如果${\omega _i}$大于给定阈值$Thresh$, 表示点${p_i}$处于高曲率区域，本文称${p_i}$为候选特征点。以经典模型Smooth-feature点云模型为例来显示PCA方法估计的法矢和特征点，其结果如图 1所示。图 1(a)为Smooth-feature模型的原型，图 1(b)(c)显示该模型各点云对应的法矢，其中图 1(c)为图 1(b)矩形框中折角面的放大图，红色线段表示的法向量，且方向调整为一致方向，大小归一化。图 1(d)红色点表示该模型的特征点。

对于光滑曲面上的点，PCA方法能比较准确地估算点云法向量；对特征曲面，特征点(两个或者多个不连续曲面相交的边界点及其周围的点) 的$k$近邻位于不同的曲面上，这种邻域称之为各向同性邻域，该邻域拟合的曲面与点云所在的局部曲面偏差较大，从而估计的法向与标准法向偏差较大。如图 1(b)(c)矩形框所示的法向量。图 2为一个立方体模型点云示意图，图 2(a)中几个不同平面交界的边界点及边界点周围的点的$k$邻域位于不同的平面上，如黄色点所表示的特征点，黄色点的$k$近邻位于3个不同的平面上，然而，黄色点理想的邻域点应该只在同一个平面上(为绿色、蓝色和红色任意一种颜色表示的点)，这种邻域称之为各向异性邻域。为了能更准确地估算点云法向量，将特征点及其各向同性邻域映射到高斯球，在高斯球上对数据点进行分割，将各向同性邻域分割为多个各向异性子邻域，在最优子邻域内拟合曲面，再次以PCA方法精确估算特征点云法向量。

设$T=T (u, v) \subset {{\rm{\boldsymbol{R}}}^3}$为一正则曲面，该曲面的法矢方向一致，将点$x \in T$处曲面的单位法向量的始端平移到单位球

$ {T^2} = \left\{ {\left( {x,y,z} \right) \subset {{\bf{R}}^3}\left| {{x^2} + {y^2} + {z^2} = 1} \right.} \right\} $

的中心，使曲面上的点与球面上的点相对应，这种对应被称为曲面$T$的高斯映射$G:T \to {T^2}$。曲面$T$的高斯映射的像记为$G\left (T \right)$，被称作$T$的高斯图，单位球${T^2}$则被称为高斯球。

高斯球具有以下性质：

如果一连通曲面的所有点都是平点，则该曲面的高斯图为高斯球上的一点；反之，如果一连通曲面的高斯图为高斯球上的一点，则该曲面的所有点都是平点。

在实际应用中，估算的点云法向量与标准法向量总是存在一定的偏差，所以连通曲面的平点映射为高斯图上的一个密集分布的数据簇，如图 3(a)所示。如果法向量估算比较准确，对于一个折角面，两曲面交界的边界点的邻域位于两个不连续的曲面上，其边界特征点邻域映射为高斯球上两个密集的独立数据簇；如果法向量估计不准确，特征点的邻域映射为高斯球上稀疏分布的数据点，如图 3(b)示意图所示，紫色点表示折角曲面中一个面上点云映射结果，红色点表示另一个面点云映射结果。然而，高斯球上来自同一连续曲面的数据点的相似性大于来自不同曲面的相似性，对高斯球上的数据点采用聚类法对其分类，将来自同一曲面的数据点分为一类，从而将特征点的各向同性邻域分割为多个各向异性子邻域。

由于PCA方法估算的法向量方向不一致，采用文献[5]方法将法向量调整为一致方向，且对法向量进行归一化。设点集$\boldsymbol{P}=\{ {p_1}, {p_2}, \ldots, {p_n}\} $的特征点${p_i}$及其${k_1}$邻域点云表示为${\boldsymbol{P}_{{\rm{sub}}}}=\{ {p_j}, j=0, 1, \ldots, {k_1}\} $。点集${\boldsymbol{P}_{{\rm{sub}}}}$的高斯图表示为

$ \begin{array}{*{20}{c}} {G\left( {{P_{{\rm{sub}}}}} \right) = \left\{ {Q\left( p \right),p \in {\mathit{\boldsymbol{P}}_{{\rm{sub}}}}} \right\} = }\\ {\left\{ {{q_i},i = 0,1, \cdots ,{k_1}} \right\}} \end{array} $

式中，$\boldsymbol{Q}(p)$为点$p$的法矢，${\boldsymbol{q}_i}$为点集${\boldsymbol{P}_{{\rm{sub}}}}$中第$j$个点的单位法向量，且方向一致。由于事先无法知道每个特征点的邻域位于多少个不同曲面上，因而无法直接确定将高斯球上数据分为多少类。均值漂移聚类法和层次聚类法是两种应用比较广泛的无参聚类方法，无需事先指定聚类数目就能自动将数据分为若干类，然而，这两种方法均不能有效的应用于此。例如，对于一个折角面，特征点的${k_1}$邻域位于两个不同曲面上，因而需将其各向同性邻域分割为两个子类，但是无参聚类法将高斯球上的数据类分割为3个及更多的数据簇，有些子数据簇的个数只有3个或者更少，从而分割的结果不准确。采用K均值聚类法将对高斯球上的数据分类，K均值聚类法能将高斯球上的数据分成$K$个子集，使子集中的点在空间上相似。$K$均值聚类法对高斯球上数据分割步骤如下：

1) 聚类数目初始化$K$=2；

2) 聚类中心初始化，计算高斯球上各数据点到聚类中心的距离，将欧氏距离最近的点归为一类，数据被分为$K$类；

3) 选择最大子类对应的邻域点为最优子邻域，最小二乘法对该邻域拟合切平面；

4) 计算最优邻域点到切平面的冗余度(冗余度表示点到切平面的距离)，并计算冗余度之和${R_{{\rm{sum}}}}$；

5) 判断${R_{{\rm{sum}}}}$是否小于$TH$，如果${R_{{\rm{sum}}}} < TH$, 切平面对应的法向量为特征点的法向量，转步骤6)；否则，$K={\rm{ }}K + 1$，重复步骤2)-5)，直到$K$=5时。当$K$>5转步骤6)；

6) 结束。

由于事先不能确定特征点的各向同性邻域位于几个曲面上，但特征点的邻域至少位于2个不同曲面上，因而K均值聚类初始化时$K$=2，将数据分为两类。K均值聚类法还需要初始化聚类中心，聚类结果依赖于初始化聚类中心。最简单的聚类中心初始化方法是随机选取$K$个数据作为初始聚类中心，但聚类结果不稳定。本文采用K-means ++算法中选择种子点的方法，该方法选择初始化聚类中心的基本思想是：初始的聚类中心之间的相互距离要尽可能的远。一般来说，当选择的初始聚类中心的相互距离大，则聚类中心点来自于不同的曲面，从而聚类结果更准确。此外，步骤5) 中，设置$K$=5时终止循环，其目的是避免出现邻域分割时不收敛情况，导致得到的邻域比较小，拟合的切平面不准确。

采用C++语言编程在Microsoft visual studio 2008上实现法向估计算法，并且调用OpenGL库函数显示点云，实验所用的计算机配置为Intel Core 2.30 GHz CPU，1.19 GB内存。为了验证算法的有效性，测试估计法矢与标准法矢偏差(估计误差)，将估计的法矢作为点云处理的输入来验证特征保留效果。并且将本文方法与文献[5]的PCA方法、文献[7-8]的加权邻域拟合法进行比较较，这些方法是经典且应用广泛的散乱点云法矢估计方法。

本文算法有4个参数需要设定，分别是$Thresh, k, {k_1}$和$TH$。$Thresh$为特征点选择阈值，需要用户设定。$k$为PCA方法粗略估算点云法向量的邻域大小，一般取$k$=864，本文取$k=16, {k_1}$为特征点映射到高斯球上的邻域，${k_1} > k$，默认取${k_1}=5k=80$。$TH$是一个关键阈值，若$TH$取较大阈值，对于两曲面交界的边界特征点邻域分割较准确，从而估算的法向量更准确；对两个以上相交的转角面周围的特征点的法向量估算不准确。相反，若$TH$取较小值，对转角面上的特征点的法向量估算更准确。然而，全局的$TH$不能满足要求，需要局部阈值。本文选取阈值$TH$的方法是：以最小二乘法对平坦点的邻域(邻域点均为非特征点) 拟合平面，计算各邻域到拟合平面的冗余度之和${R^i}_{sum}$(${R^i}_{sum}$表示第$i$个平坦点的${R}_{sum}$)，以各平坦点的${R^i}_{sum}$的平均值作为$TH$的阈值。

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

$ \begin{array}{l} \mathit{\boldsymbol{RSM = }}\sqrt {\frac{1}{{\left| S \right|}}\sum\limits_{x \in \mathit{\boldsymbol{S}}} {f{{\left( {\left\langle {{\mathit{\boldsymbol{n}}_{{\rm{ref}}}},{\mathit{\boldsymbol{n}}_{{\rm{est}}}}} \right\rangle } \right)}^2}} } \\ f\left( {\left\langle {{\mathit{\boldsymbol{n}}_{{\rm{ref}}}},{\mathit{\boldsymbol{n}}_{{\rm{est}}}}} \right\rangle } \right) = \left\{ \begin{array}{l} \left\langle {{\mathit{\boldsymbol{n}}_{{\rm{ref}}}},{\mathit{\boldsymbol{n}}_{{\rm{est}}}}} \right\rangle \;\;\;\;\left\langle {{\mathit{\boldsymbol{n}}_{{\rm{ref}}}},{\mathit{\boldsymbol{n}}_{{\rm{est}}}}} \right\rangle < \tau \\ \frac{{\rm{ \mathsf{ π} }}}{2}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. \end{array} $

(4)

式中，${\boldsymbol{n}_{{\rm{ref}}}}、{\boldsymbol{n}_{{\rm{est}}}}$分别为单位化的参考法向和估计法向，|$S$|为曲面$S$的点云数量，$\tau $为阈值，当参考法向与估计法向的点积大于阈值$\tau $时，认为该点的法向与标准法向偏差为90°。文献[8-10]均取法向夹角阈值为10°，即$\tau $=0.984 8。

因为复杂尖锐特征曲面难以得到标准法向，所以用MATLAB合成规则的立方体点云模型(不含噪声)，根据点云坐标计算其标准法向。合成模型模型的点云数量分别为6 147，其特征点为两平面交叉的边界特征点和三平面交叉的角点。为了测试估计法向对噪声的敏感性，对理想模型添加不同程度的噪声，计算理想模型及各噪声模型的RMS。当估计的法向与标准法向偏差大于10°时称该点为坏点，4种方法比较结果如图 4所示，第1行为理想数据，第2、3行是理想模型添加了高斯白噪声的结果，其信噪比SNR分别是40 dB、30 dB，图中红色点表示坏点。表 1为两个模型的RMS和和坏点数量(BPN)。

当模型不含噪声时本文方法几乎没有坏点，文献[8]在特征点处有少数的坏点。当SNR=40 dB时，文献[8]有少数坏点，PCA方法和文献[7]的坏点主要是特征点及其周围的点。当数据被噪声严重干扰时(SNR=30 dB)，如图 4中第3行所示，立方体已经严重失真，角点和边界点模糊。PCA和文献[7]方法估计法向有一半是坏点(BPN=306 4，详见表 1)，本文方法得到BPN，RMS最小，其坏点主要分布在角点周围。总之，PCA方法和文献[7]方法在特征点出估计的法向量与标准法向量的偏差较大，文献[7]虽然对邻域加权改进了PCA方法，但是文献[7]对邻域加权只考虑了距离权重，因而改进后取得的效果不佳，文献[8]对邻域加权时不仅考虑了距离对当前点的作用，还考虑了邻域点到拟合曲面的残差值，因而在特征点处估算的法向量更准确，但是文献[8]方法需人工确认权重带宽，且权重带宽为全局参数，不适合每个不同的特征点。本文方法将特征点的各向同性邻域进行分割，以各向异性邻域拟合曲面，从而得到的法向量更准确。

为了更清楚估计法矢的准确性，对估计的法矢可视化，并且将估计的法矢作为点云曲面重建算法输入参数来测试特征保留效果。

对两个经典模型(八面体和Fandisk模型，其中八面体模型添加了噪声，SNR=40 dB) 的法向可视化，模型上的线条表示一致朝外的法向量，八面体模型的特征点包括多个平面交接的角点和两个平面交接的过渡边界点，角点处的曲面比较尖锐，过渡边相对平缓。Fandisk模型由平面、球面、圆柱面等基面组成，其过渡边有尖锐边和平缓边等。为了更清楚地显示模型，将模型三角化加光照显示，如图 5(a)、图 6(a)所示。图 5、图 6中，从左到右列分别是点云三角化曲面模型，PCA、文献[7]、文献[8]和本文方法估计的法向。准确的法向应该垂直于点云所在的局部曲面，两个不同曲面上点云法向夹角为两曲面的夹角，局部区域同一曲面上的法向接近平行，在两个曲面交界的边界处有清晰的边界。图 5、图 6中(b)(c) 的边界几乎无法辨认，图 5(d)、图 6(d)有部分特征点的法向误差大，如角点的法向不垂直于其所在的曲面。图 7(a)为图 6(a)矩形框的区域的放大图，图 7(b)-(e)为图 6(a)曲面对应方法估算的法向量。

将估计的法向应用到点云处理中，根据点云数据进行曲面重建。在逆向工程应用领域中，通过3维扫描仪扫描实物获得其3维点云，通过点云数据处理及曲面重建生成实物的CAD模型，可以实现产品的缺陷修复和再制造。在点云数据处理过程中许多点云处理算法都依赖准确的点云法向量，若点云法向量估算不准确，直接影响曲面重建效果，造成再制造产品与原始产品的较大偏差。为了验证本文法向量估算的有效性，将估算的法向量应用到点云曲面重建中。隐式移动最小二乘曲面拟合法是一种应用比较广泛的曲面重建算法^[14]，若点云法向量估算准确该方法能够比较准确的重建尖锐特征曲面，即在特征区域能够进行特征保留，即准确的法向是该算法特征保留的前提。若法向在特征点处被平滑，则曲面重建后尖锐特征也被平滑，从而使重建的模型与原始模型在特征区域处生成较大误差。本文扫描模型A实物，获得该实物的点云，再对该模型采用文献[14]方法对其进行曲面重建，以不同方法估算的法向量作为该算法的输入参数之一，为了消除其他因素的影响，曲面重建算法中除了法向不一样，其他参数值均相同。计算重建模型与原始模型的误差，在不同法向量的参数下曲面重建误差如图 8所示，从图 8误差比较结果可知本文方法的法向量作为曲面重建输入参数重建误差在特征区域最小，进一步验证了本文方法估算法向量的有效性。

本文对特征曲面的法向估计进行了研究。用PCA方法对点云的$k$邻域粗略估计曲面法向和特征点，PCA方法估计的点云法矢在平坦区域比较准确，但是在特征区域，特征点的$k$邻域为各向同性，估算的法向量不准确。将特征点的各向同性邻域映射到高斯球，在高斯球上采用K均值聚类法对数据进行分类得到各向异性子邻域，以各向异性邻域拟合曲面精确估算特征点的法向量。本文方法能够比较准确地估计特征曲面的法向，最小误差接近0，且通过与其他方法估算的法向量的定量(误差) 和定性(可视化和应用) 比较验证了本文方法估算法向量的有效性。然而，对尖锐特征模型点云法向量估算，本文方法对规则模型具有较好的适用性，对复杂模型的鲁棒性有待提高，在未来工作中将研究复杂模型的尖锐特征曲面的法向量估算。