点云数据的主曲率和主方向估计方法

摘要

本发明涉及一种点云数据的主曲率和主方向估计方法步骤包括：预处理，法截线曲率估计，拟合出Weingarten矩阵，求出Weingarten矩阵的特征值和特征向量，求出主曲率和主方向。本发明仅利用激光扫描仪的扫描数据和估计的法向量，得到忠实于原始实物的主曲率和主方向。该方法通过最小二乘线性拟合和矩阵的特征值特征向量求出主曲率和主方向，算法简单，计算结果准确，时间复杂度是高效的。该方法称为法截线拟合法，其计算结果在虚拟现实、电脑游戏、自然场景模拟、城市景观设计、数据压缩、特征提取、实物3D重建等领域具有重要的应用价值。

说明书

技术领域

本发明涉及微分几何、计算数学、计算机图形学和计算机视觉技术领域的一种利用三维激光扫描仪进行实物测量得到点云数据，并根据点云数据来进行主曲率和主方向计算的方法。在虚拟现实、电脑游戏、自然场景模拟、城市景观设计、数据压缩、特征提取、实物3D重建等领域具有重要的应用价值。

背景技术

随着激光扫描仪精度的提高，扫描得到的信息越来越丰富，扫描得到的模型数据越来越庞大。人们利用这些庞大的数据进行特征提取、数据压缩或者进行三维重建。但是，这些工作的实现，往往需要对一些微分几何量进行估计，其中最重要的估计包括主曲率和主方向的估计。

目前主曲率和主方向的估计方法大致可以分为三大类：

第一类是进行局部曲面拟合，求出二次或者三次曲面，再根据微分几何理论求出主曲率和主方向，如2004年Goldfeather提出的三次曲面拟合方法；

第二类是先对点云数据进行网格化，再取一个或者二个环内的近邻点进行主曲率或者主方向的估计，如1995年Taubin提出的基于网格数据的方法；

第三大类是直接对点云数据进行微分几何特征量的计算。前两类方法在过去的一些应用中效果很好，但是随着扫描仪精度的提高，数据规模的日益增大，这两类方法的时间空间开销较大，他们的应用就受到了制约。于是人们试图直接从点云数据进行微分几何特征的计算，也就是第三类方法。对于目前的成果来说，直接从点云进行计算时，有些方法只是利用点的位置信息，没有利用各个点的法向量信息，这样会使得方法的鲁棒性较差；也有一些方法用了法向量信息，但是他们往往结合第一类方法，把法向量作为一个约束条件加以利用，增加了主曲率和主方向模型求解的计算时间和存储空间的开销。

发明内容

本发明欲解决的技术问题是离散点云数据的主曲率和主方向的准确计算，本发明的目的在于，针对现实世界中由激光扫描得到的离散点云数据，提供一个对主曲率和主方向进行准确、快速计算的、鲁棒的点云数据的主曲率和主方向估计方法。

为实现上述目的，本发明的技术解决方案是提供一种点云数据的主曲率和主方向的估计方法，称为法截线拟合法，该主曲率和主方向估计步骤包括：

步骤1：利用激光扫描仪扫描直接采集点云数据并对点云数据预处理，按照点云数据中每个点的坐标进行空间划分，实现三维空间的二分查找树的数据存储结构称为kd树(k-dimensional tree)；所述每个点的坐标，都是采用激光扫描仪产生的原始坐标。

步骤2：对于点云数据的每一个点，利用点云数据的kd树查找15个或30个近邻点，根据最小二乘方法把这些近邻点拟合出一个平面，以这个平面的法向量作为点云法向量的初始估计值，然后通过加权平均算法修正点云数据的各个点的法向量估计，是以点云数据中每一个点与近邻点的欧式距离的倒数作为权值，所述欧式距离的倒数，若距离值为0，则不计算这个近邻点。所述根据最小二乘方法把这些近邻点拟合出一个平面，其中拟合平面时需要构造带权最小二乘问题，这个最小二乘问题是这些近邻点到拟合平面的残差的绝对值，再乘以权系数的积的和的最小值问题。所述权系数，是以点云数据中每一个点与近邻点的欧式距离的倒数作为权值。

步骤3：对于点云数据的每一个点，利用其法向量、切平面构造局部三维直角坐标系；所述局部三维直角坐标系，其构造方法是，对于点云数据中一个点p，若点p的法向量为N＝(n_x，p，n_y，p，n_z，p)，则这个点p就是局部坐标系的原点，X，Y，Z三个坐标轴分别为Z＝N＝(n_x，p，n_y，p，n_z，p)，其中表示法向量N的第一方向角，ψ＝arccos(n_z，p)表示法向量的第三方向角的余角。

步骤4：对于点云数据的每一个点，利用点云数据的kd树查找15个或30个近邻点；

步骤5：对于查找到的近邻点，通过三维坐标变换，把这些近邻点的原始坐标和这些近邻点的法向量都转化为局部坐标系的坐标，是把这些近邻点的坐标减去局部坐标系原点的坐标就是这些近邻点在局部坐标系中的坐标，这些近邻点的法向量坐标与局部坐标系的三个坐标轴分别作数量积，算出这些近邻点的法向量在局部坐标系中的坐标。

步骤6：利用点云数据的每一个点及其法向量、一个近邻点、近邻点的法向量构造近似三角形，根据正弦定理给出点云的法截线的法曲率的近似表达式；所述构造近似三角形，该三角形由角角边定理确定：设点云中的待计算法曲率的点为p，它的一个近邻点为q_i点p和q_i的法向量分别为N和M_i，向量N和M_i的夹角为β，向量N和向量的夹角为α，线段|pq_i|为边，则角β，角α，边|pq_i|确定了这个三角形。所述用正弦定理给出点云的法截线的法曲率的近似表达式为其中k_nⁱ表示第i近邻点所对应的法截线的法曲率，β为近邻点法向量与局部坐标系Z轴的夹角，α为近邻点的定位向量与Z轴夹角的补角，局部坐标系的原点为p，近邻点为q_i，|pq_i|表示点q_i与点p的欧式距离。所述这个表达式在实际计算中采用近似计算形式其中，n_z＝n_z，i，近邻点q_i的局部坐标为(x_i，y_i，z_i)，近邻点q_i的法向量的局部坐标为(n_x，i，n_y，i，n_z，i)。

步骤7：在局部坐标系中，利用法曲率，根据欧拉公式(Euler Equation)构造非线性最优化问题，通过三角形公式进行恒等变换，把这个非线性最优化问题转化为线性拟合，求出韦恩伽汀矩阵(Weingarten矩阵)的各个元素；所述根据欧拉公式(Euler Equation)构造非线性最优化问题，设需要求点云数据中一个点的主曲率和主方向，则是构造的非线性最优化问题，其中k₁，k₂是待求的主曲率，未知数θ为局部坐标系中X轴与最大主曲率对应的主方向的夹角，θ_i为近邻点q_i在局部坐标系中xOy面的投影点的定位向量与局部坐标系中X轴的夹角。所述把这个非线性最优化问题转化为线性拟合，就是把所有近邻点的法曲率作为因变量观测值，所有近邻点与局部坐标系的x轴的方向角的余弦值的平方、方向角的2倍的正弦值、方向角的正弦值的平方作为自变量的观测值，作三元线性拟合。三个拟合系数值依次作为二阶对称矩阵Weingarten矩阵的上三角的三个元素。

步骤8：利用矩阵的奇异值分解(SVD分解)求出Weingarten矩阵的特征值和特征向量；

步骤9：利用Weingarten矩阵的特征值和特征向量求出主曲率和主方向。所述求出Weingarten矩阵的特征值和特征向量，若两个特征值相等，则该点为脐点，特征向量取为(1，0)和(0，1)。所述利用Weingarten矩阵的特征值求出主曲率，是用两个特征值作为主曲率。所述利用Weingarten矩阵的特征向量求主方向，是分别用两个特征向量的两个分量作为组合系数，作局部坐标系的X轴和Y轴的线性组合，作为两个主方向。

本发明的有益效果是提出对第三类技术改进的新方法，是在点云数据上直接计算主曲率和主方向。利用点云的坐标信息和点云的法向量信息。本发明与前人方法的区别主要体现在，法向量被直接用确定法截线，通过法截线的法曲率来进行主曲率估计，而不是作为非线性规划问题的约束条件，此外本发明还把非线性最优化问题归结为求解一个线性系统。因此，本发明计算的时间和空间性能也较以往方法优越，实验表明，主方向的计算结果也比其他方法准确。我们利用三维激光扫描仪获取实物的表面信息，并根据扫描数据利用常规方法计算出各个点的法向量，再给出法截线的法曲率估计式，最后经过线性拟合和韦恩伽汀(Weingarten)矩阵计算求出主曲率和主方向。本发明所获得的计算结果能用于计算机图形学各应用领域，包括特征增强、表面简化、网格优化、特征提取、数据压缩、骨架提取和3D实物重建、虚拟现实、电脑游戏、自然场景模拟、城市景观设计、数据压缩等领域具有重要的应用价值。利用本发明，可以快速、方便、准确地计算出点云数据的主曲率和主方向。

附图说明

图1算法流程图

图2近邻点及其法向量所确定的三角形

图3局部空间直角坐标系和法截线示意图

图4圆环面图及其添加噪声之后的效果图

图5圆环面主曲率估计结果的误差比较

图6无噪声手模型的主方向估计结果

图7有噪声手模型的主方向估计结果

图8无噪声大象模型的主方向估计结果

图9有噪声大象模型的主方向估计结果

图10树枝模型点云数据

图11树枝模型点云数据的主方向

图12树枝模型点云数据树干表面的主方向

图13树枝模型点云数据树干分叉处的主方向

图14树枝模型点云数据树干凹陷处的主方向

具体实施方式

下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

1、方法概述(overview of approach)

如图1示出本发明整个算法的流程，本发明算法的主要步骤包括：

1)、数据预处理，创建kd树，计算点云数据各个点的法向量；

2)、逐点计算主曲率和主方向，包括7个子步骤：(a)建立局部坐标系、(b)搜索近邻点、(c)近邻点坐标局部化、(d)计算每个近邻点的对应的法截线的法曲率、(e)线性拟合出韦恩伽汀(Weingarten)矩阵的元素、(f)求出Weingarten矩阵的特征值和特征向量、(g)求出主曲率和主方向。

3)、输出计算结果，存储点的主曲率和主方向，显示点的主方向。

2、数据预处理

三维点云数据一般只有点的坐标信息。要估计主曲率和主方向，应该先估计各个点的法向量。而求法向量和主曲率、主方向都需要用到近邻点，从计算几何理论可以得知，查找近邻点的最快捷方法之一就是建立点云数据的3维空间的二分查找树，简称为kd树(k-dimensional tree)，这是查找效率很高的数据存储结构。于是，预处理阶段需要完成kd树的建立和法向量的计算。

首先，建立kd树。在计算几何中，kd树是已经被证明的查找近邻的最快捷的数据结构之一。kd树基于点的空间位置信息，通过二分法迭代划分三维空间，实现最优存储。在kd树上，进行k近邻查找的时间复杂度为O(log₂n)，这里n为点云数据的点的个数。

其次，计算各个点的法向量。对于点云数据的每一个点，利用点云数据的kd树查找15个或30个近邻点，假设这些近邻点来自于同一个平面，于是可以用这些近邻点到拟合平面的残差的绝对值，再乘以权系数的积的和构造最小二乘问题，其中的权的确定是以点云数据中每一个点与近邻点的欧式距离的倒数作为权值。利用最小二乘方法拟合出这个平面，以这个平面的法向量作为点云的法向量的初始估计值，然后通过加权平均算法修正点云的法向量估计，权值等于以点云数据中每一个点与近邻点的欧式距离的倒数。

3逐点计算主曲率和主方向

快速而且准确地计算出扫描得到的点云数据中各个点的主曲率和主方向是本发明中的主要目标。计算的方法是逐点迭代进行的。计算一个点(记为p)的主曲率和主方向包括以下七个步骤：

3.1建立局部坐标系

建立局部坐标系，如图3显示了这个局部空间直角坐标系：局部三维直角坐标系，设点p的法向量为N＝(n_x，p，n_y，p，n_z，p)，则这个p点就是局部坐标系的原点，X，Y，Z三个坐标轴分别为Z＝N＝(n_x，p，n_y，p，n_z，p)，其中ψ＝arccos(n_z，p)，

3.2搜索近邻点

由于预处理阶段建立了kd树，以点p为查询点，基于欧式距离查找k近邻，比如k＝15近邻或k＝30近邻，根据数据扫描精度确定近邻点的个数，扫描精度较低则取15个近邻点，否则取30个近邻点，设定查找近邻的误差阈值为0.0，这里记q_i(i＝1，2，...，m)为查询得到的m个近邻点。

3.3近邻点坐标局部化

求出所有近邻点在局部坐标系的坐标。本发明的方法可以描述为下面几个步骤。假设待求主曲率的点为p，近邻点为q_i，以定位向量的坐标作为点q_i的局部坐标，这里坐标转化的实质是对近邻点做了一个点的平移变换。

把近邻点的法向量转化为局部坐标系的坐标。本发明的方法是预处理阶段求出的各个近邻点的法向量(记为M_i)在局部坐标系中的方向余弦作为法向量的在局部坐标系中的坐标。

3.4计算各个近邻点对应的法截线的法曲率

假定待求主曲率和主方向的点、近邻点以及他们的法向量确定了一条法截线，现在估计这条法截线在点p的法曲率。

待求主曲率和主方向的点p、近邻点q_i以及他们的法向量N和M_i确定一个近似三角形，由图2示出的近邻点O、p、q_i及其法向量N和M_i所确定的三角形。确定的方法是依据角角边定理，β为近邻点法向量与局部坐标系Z轴的夹角，α为近邻点的定位向量pq_i与Z轴夹角的补角，线段pq_i为边。在这个近似三角形，依据正弦定理取法曲率为

$k_n^i=-\sin \beta ÷\left(\right|{\mathrm{pq}}_i\left|\sin \alpha \right),$

其中k_nⁱ表示第i近邻点所对应的法截线的法曲率，|pq_i|表示点q_i与点p的欧式距离。为了便于计算，把坐标代入上面式子运用三角公式可以把这个式子转化为下面的式子，作为法曲率的估计式

$k_n^i=-n_{\mathrm{xy}}÷(\sqrt{n_{\mathrm{xy}}^2+n_z^2}·\sqrt{x_i^2+y_i^2})$

其中，近邻点q_i的局部坐标为(x_i，y_i，z_i)，近邻点q_i的法向量的局部坐标为(n_x，i，n_y，i，n_z，i)。

3.5线性回归求出Weingarten矩阵

根据欧拉公式(Euler Equation)构造非线性最优化问题，设需要求点云数据中一个点p的主曲率和主方向则是构造非线性最优化问题，其中k₁，k₂是待求的与主方向对应的主曲率，θ为未知数，是局部坐标系中X轴与最大主曲率对应的主方向e₁的夹角，点Q_i是近邻点q_i的在切平面S上的投影点，θ_i为向量pQ_i与局部坐标系中X轴的夹角。图3示出法截线和主方向和q_i、Q_i、p这些点、角度θ_i、θ。

本发明把上述非线性规划问题转化为求解下面的最小二乘问题，

$\underset{\mu }{\min }{\left|\right|\mathrm{M\mu }-R\left|\right|}^2$

其中，$M_{m\times 3}=\left(\begin{array}{ccc}{\cos }^2{\theta }_1& 2\cos {\theta }_1\sin {\theta }_1& {\sin }^2{\theta }_1\\ .& .& .\\ .& .& .\\ .& .& .\\ {\cos }^2{\theta }_i& 2{\cos \theta }_i\sin {\theta }_i& {\sin }^2{\theta }_i\\ .& .& .\\ .& .& .\\ .& .& .\\ {\cos }^2{\theta }_m& 2\cos {\theta }_m\sin {\theta }_m& {\sin }^2{\theta }_m\end{array}\right),$$R_{m\times 1}=\left(\begin{array}{c}{k}_n^1\\ .\\ .\\ .\\ k_n^i\\ .\\ .\\ .\\ k_n^m\end{array}\right),$μ＝(A，B，C)^T，A＝k₁cos²θ+k₂sin²θ，B＝(k₂-k₁)cosθsinθ，C＝k₁sin²θ+k₂cos²θ

这一步的求解可以转化为求解超定线性方程组Mμ＝R。然后，用求出来的A，B，C的值来构造Weingarten矩阵

$W=\left(\begin{array}{cc}A& B\\ B& C\end{array}\right)$

3.6计算Weingarten矩阵的特征值和特征向量

求出上一步的Weingarten矩阵的特征值和特征向量，对特征值的大小进行排序，把大的特征值记为k₁，它也是该点的最大曲率，对应的单位长度的特征向量称为较大特征值对应的特征向量，另一个较小的特征向量记为k₂，其值为该点的最小曲率，对应的单位长度的特征向量称为较小特征值对应的特征向量。

假如该点的两个特征值相等，即表示该点为脐点，而脐点的主方向是平行于切平面的任意方向，方向不是惟一的。本发明对这种特殊点的处理方法是取标准正交基(1，0)和(0，1)作为特征向量。

3.7利用特征值和特征向量求出主曲率和主方向

上一步计算出的特征向量是二维的，是相对于局部坐标系中xOy面上的向量。需要把它们转化到全局坐标系。转化的方法是把最大曲率对应的特征向量的各个分量作为组合系数，求出局部坐标系的X轴和Y轴的线性组合作为最大曲率对应的主方向，即为主方向e₁；同样地，把最小曲率对应的特征向量的各个分量作为组合系数，求出局部坐标系的X轴和Y轴的线性组合作为最小曲率对应的主方向，即为e₂。

主曲率的值与坐标系的选取无关，不需要转化，即为上一步计算出来的特征值k₁和k₂。

4结果输出

对于计算的结果，本发明提供两种形式输出。第一种方式为把点的坐标、法向量、主方向、主曲率以文件形式存储，再由专业软件进行显示。第二种方式为把主方向、主曲率信息转化为矢量信息或者颜色信息，结合点云数据显示出来。

5实验结果与结论

我们将此方法应用于解析数据和实际扫描得到的不同复杂度的数据，与之前的两种主要的主曲率主方向估计方法进行了比较。通过实验证明，本发明对主曲率和主方向的估计更加准确可靠，鲁棒性也更强。

5.1解析曲面模拟数据的对比实验

解析数据的实验用来说明本发明提出的新方法比之前的各种方法更加准确可靠。解析数据来自一个圆环解析曲面，它的解析式是已知的，本发明用的圆环面方程为：

r(u，v)＝((R+rcosu)cos v，(R+rcosu)sin v，rsinu)

(0≤u≤2π，0≤v≤2π)

其中参数R＝2，r＝1。从这个解析曲面，按照二维均匀分布获取5000个随机采样点，根据微分几何公式可以计算出各个点的主曲率和主方向(称为采样点的实际主曲率和实际主方向)。然后，把这5000个采样点视作扫描得到的点云数据。对每一个点p，沿着法向量方向添加噪声这里为点p的法向量，a服从[-mh，mh]的均匀分布，m是任意两点之间的距离的中位数，h为噪声水平，实验中用过的噪声水平为h＝0.1，0.2，...，1.0。图4显示了圆环面图及其添加噪声之后的效果图，其中图4(a)为初始圆环，图4(b)为添加噪声的圆环面。

与本发明进行对比实验的方法是Goldfeather立方曲面拟合方法和Taubin方法。2004年公开发表的Goldfeather立方曲面拟合方法是目前公布的最准确方法，这种方法与本发明的法截线拟合方法在机理上是不同的，它实际上是对待求主曲率的点的局部进行三次曲面拟合，然后根据解析方法求出主曲率和主方向，这种方法把法向量作为约束条件加以考虑，因此比以往的所有估计方法都更加准确，鲁棒性也更好。1995年公开发布的Taubin方法则利用近邻点构造多边形计算曲率张量，求出主曲率和主方向。这种方法是曲率估计的典型算法，其效果并不比2003年之前公布的改进算法或者新算法差多少，这些算法不用法向量进行曲率估计，计算结果的准确性较差。因此，本发明以Taubin方法作为这些算法的一个代表。

对比实验在个人台式电脑上执行，电脑的配置为Intel(R)Core(TM)2CPU，4400@2G，1.99GHz，2.0G内存。编程语言为C++和Opengl glut3.7。

实验显示的是两个主曲率构造的几何量，即平均曲率和高斯曲率。图5显示了曲率高斯曲率和平均曲率估计结果的误差比较。

从图5(a)高斯曲率估计误差和图5(b)平均曲率估计误差可以看出，本发明的估计结果比Taubin方法误差要小许多，噪声水平越高，优越性越显著。从图形上看本发明的方法的估计结果与Goldfeath方法的估计结果比较接近，因此把实验结果的数据列表如下。

表1平均曲率估计绝对误差对比表

从表1可以看出，本发明与Goldfeath方法相比，除了在噪声水平较低时平均曲率误差略大以外，其余情况误差更小，说明计算更加准确。对于噪声较大的情况，本发明优势更加明显。

5.2点云数据的对比实验

把本发明与Goldfeath方法应用于实际扫描数据，用了两个模型进行测试。

第一个模型是一个“手”的数据。该模型的高度为202个长度单位，有6258个随机采样点，数据本身带有较准确的法向量。图6为无噪声手模型的主方向估计结果，显示初始数据图以及本发明与Goldfeath方法的主方向估计结果。给数据点沿着法向量方向加入噪声，噪声水平为2.0个单位，则法向量与点的位置信息不是完全正确，求出的主方向显示在图7为有噪声手模型的主方向估计结果。图6(a)为原始数据，图6(b)为本发明的实验计算结果，图6(c)为Goldfeath方法的计算结果；图7(a)为原始数据；图7(b)为本发明的计算结果，图7(c)为Goldfeath方法的计算结果数据。从这些图可以看出，在噪声水平较低时，本发明的方法与Goldfeath方法都有较好的结果，本发明的方法比Goldfeath方法的优势没有特别显著，但是当噪声水平加大时，Goldfeath方法估计的主方向明显变得零乱。比如图7中，拇指的背部，本发明的方法同图6中的相差不大，依旧体现原有次序，而Goldfeath方法却较凌乱。这个实验说明本发明比Goldfeath方法更加准确，鲁棒性也更好。

另一个模型数据为“大象”。这个点云模型高度为75个长度单位，包含6859个随机采样点。实验方法同“手”模型。图8是无噪声大象模型及其在不同方法下得到的主方向估计结果，图9是有噪声大象模型及其在不同方法下得到的主方向估计结果，其为噪声水平为2.0的数据。观察图9中大象的背部，可以看出，本发明的方法同无噪声的情况相差不大，依旧体现原有方向，而Goldfeath方法比较凌乱。这个实验进一步说明本发明比Goldfeath方法更加准确，更鲁棒。图8(a)为原始数据；图8(b)为本发明的数据，图8(c)为Goldfeath方法的数据，图9(a)为原始数据；图9(b)为本发明的数据，图9(c)为Goldfeath方法的数据。

5.3本发明在复杂点云数据(树枝模型)中的应用

树的结构比一般的点云模型复杂许多，主要体现在以下三点：由于遮挡，点的分布不均匀，有时成团，有时成线，有时成簇；表面变化较大，树干的粗细变化较大；枝的生长变化丰富，导致连接关系复杂。因此，树枝的主方向估计很有难度。

本发明用到的树枝点云模型是树高12米，除了主枝外还有四个分枝。该点云模型数据有6950个点。实验效果如图10至图14。图10为树枝模型原始的点云数据。图11为树枝模型点云数据各个点的主方向图。图12显示了树枝模型点云数据树干表面各个点的主方向，其特征是主曲率较大者对应的主方向体现了树干的横截面方向，而主曲率较小者对应的主方向刻画了树干的生长方向。图13显示了树枝模型点云数据树枝分叉点处的主方向特征。图14显示了树枝模型点云数据树干表面凹陷处的主方向特征，这些主曲率较大者对应的主方向指向凹陷的中心。

本发明提出的主曲率和主方向计算方法的特征在于利用点云数据和各个点的法向量计算点的主曲率和主方向。

上述实验结果和利用点云数据求取主方向和主曲率的方法，可以用于计算机图形学各应用领域，具有高可信度、操作简单、应用前景广的特点。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。