专利名称：点云数据的主曲率和主方向估计方法
技术领域：
本发明涉及微分几何、计算数学、计算机图形学和计算机视觉技术领域的一种利
用三维激光扫描仪进行实物测量得到点云数据，并根据点云数据来进行主曲率和主方向计算的方法。在虚拟现实、电脑游戏、自然场景模拟、城市景观设计、数据压縮、特征提取、实物
3D重建等领域具有重要的应用价值。
背景技术：
随着激光扫描仪精度的提高，扫描得到的信息越来越丰富，扫描得到的模型数据越来越庞大。人们利用这些庞大的数据进行特征提取、数据压縮或者进行三维重建。但是，这些工作的实现，往往需要对一些微分几何量进行估计，其中最重要的估计包括主曲率和主方向的估计。
目前主曲率和主方向的估计方法大致可以分为三大类 第一类是进行局部曲面拟合，求出二次或者三次曲面，再根据微分几何理论求出主曲率和主方向，如2004年Goldfeather提出的三次曲面拟合方法； 第二类是先对点云数据进行网格化，再取一个或者二个环内的近邻点进行主曲率或者主方向的估计，如1995年Taubin提出的基于网格数据的方法； 第三大类是直接对点云数据进行微分几何特征量的计算。前两类方法在过去的一些应用中效果很好，但是随着扫描仪精度的提高，数据规模的日益增大，这两类方法的时间空间开销较大，他们的应用就受到了制约。于是人们试图直接从点云数据进行微分几何特征的计算，也就是第三类方法。对于目前的成果来说，直接从点云进行计算时，有些方法只是利用点的位置信息，没有利用各个点的法向量信息，这样会使得方法的鲁棒性较差；也有一些方法用了法向量信息，但是他们往往结合第一类方法，把法向量作为一个约束条件加以利用，增加了主曲率和主方向模型求解的计算时间和存储空间的开销。
发明内容
本发明欲解决的技术问题是离散点云数据的主曲率和主方向的准确计算，本发明的目的在于，针对现实世界中由激光扫描得到的离散点云数据，提供一个对主曲率和主方向进行准确、快速计算的、鲁棒的点云数据的主曲率和主方向估计方法。 为实现上述目的，本发明的技术解决方案是提供一种点云数据的主曲率和主方向的估计方法，称为法截线拟合法，该主曲率和主方向估计步骤包括 步骤1 :利用激光扫描仪扫描直接采集点云数据并对点云数据预处理，按照点云
数据中每个点的坐标进行空间划分，实现三维空间的二分查找树的数据存储结构称为kd树(k-dimensional tree);所述每个点的坐标，都是采用激光扫描仪产生的原始坐标。
步骤2 :对于点云数据的每一个点，利用点云数据的kd树查找15个或30个近邻点，根据最小二乘方法把这些近邻点拟合出一个平面，以这个平面的法向量作为点云法向量的初始估计值，然后通过加权平均算法修正点云数据的各个点的法向量估计，是以点云
4数据中每一个点与近邻点的欧式距离的倒数作为权值，所述欧式距离的倒数，若距离值为O，则不计算这个近邻点。所述根据最小二乘方法把这些近邻点拟合出一个平面，其中拟合平面时需要构造带权最小二乘问题，这个最小二乘问题是这些近邻点到拟合平面的残差的绝对值，再乘以权系数的积的和的最小值问题。所述权系数，是以点云数据中每一个点与近邻点的欧式距离的倒数作为权值。 步骤3 :对于点云数据的每一个点，利用其法向量、切平面构造局部三维直角坐标系；所述局部三维直角坐标系，其构造方法是，对于点云数据中一个点P，若点P的法向量为N = (nx,p， ny,p， nz,p)，则这个点p就是局部坐标系的原点，X， Y， Z三个坐标轴分别为X = (—sin p， cos p， 0) ，y = (cos y cos p， cos —n p， _ sin ，z = N= (nx,p，ny,p，nz,p)，其
中p-arctan("^/"w)表示法向量N的第一方向角，v = arccos (nz,p)表示法向量的第三方向角的余角。 步骤4 :对于点云数据的每一个点，利用点云数据的kd树查找15个或30个近邻点； 步骤5 :对于查找到的近邻点，通过三维坐标变换，把这些近邻点的原始坐标和这些近邻点的法向量都转化为局部坐标系的坐标，是把这些近邻点的坐标减去局部坐标系原点的坐标就是这些近邻点在局部坐标系中的坐标，这些近邻点的法向量坐标与局部坐标系的三个坐标轴分别作数量积，算出这些近邻点的法向量在局部坐标系中的坐标。
步骤6 :利用点云数据的每一个点及其法向量、一个近邻点、近邻点的法向量构造近似三角形，根据正弦定理给出点云的法截线的法曲率的近似表达式；所述构造近似三角形，该三角形由角角边定理确定设点云中的待计算法曲率的点为P，它的一个近邻点为点P和的法向量分别为N和Mi，向量N和Mi的夹角为13 ，向量N和向量^J,.的夹角为
a，线段|pq」为边，则角P，角a，边|pq」确定了这个三角形。所述用正弦定理给出点云的法截线的法曲率的近似表达式为《=-sin/ + (l /^Isino;)，其中C表示第i近邻点所对应的法截线的法曲率，P为近邻点法向量与局部坐标系Z轴的夹角，a为近邻点的定位向量与Z轴夹角的补角，局部坐标系的原点为p，近邻点为qi， |pqi|表示点A与点p的欧式距离。所述^=-sin/ + (l/^,lsinoO，这个表达式在实际计算中采用近似计算形式
*二 = —"v + M:2 + "W)，其中，/^ = (x, + +V^ +凡2 ,nz = nz,i，近令卩点
qi的局部坐标为(Xi，yi，Zi)，近邻点qi的法向量的局部坐标为(r^,i，ny,i，r^i)。 步骤7 :在局部坐标系中，利用法曲率，根据欧拉公式(Euler Equation)构造
非线性最优化问题，通过三角形公式进行恒等变换，把这个非线性最优化问题转化为线
性拟合，求出韦恩伽汀矩阵(Weingarten矩阵)的各个元素；所述根据欧拉公式(Euler
Equation)构造非线性最优化问题，设需要求点云数据中一个点的主曲率和主方向，则
4J^。
g[M。s2^") + Min^ + e)-A:f是构造的非线性最优化问题，其中、，k2是待求的主曲
率，未知数e为局部坐标系中X轴与最大主曲率对应的主方向的夹角，9i为近邻点qi在局部坐标系中xOy面的投影点的定位向量与局部坐标系中X轴的夹角。所述把这个非线性最优化问题转化为线性拟合，就是把所有近邻点的法曲率作为因变量观测值，所有近邻点与局部坐标系的x轴的方向角的余弦值的平方、方向角的2倍的正弦值、方向角的正弦值的平方作为自变量的观测值，作三元线性拟合。三个拟合系数值依次作为二阶对称矩阵Weingarten矩阵的上三角的三个元素。 步骤8 :利用矩阵的奇异值分解(SVD分解)求出Weingarten矩阵的特征值和特征向量； 步骤9 :利用Weingarten矩阵的特征值和特征向量求出主曲率和主方向。所述求出Weingarten矩阵的特征值和特征向量，若两个特征值相等，则该点为脐点，特征向量取为(l，O)和(O，l)。所述利用Weingarten矩阵的特征值求出主曲率，是用两个特征值作为主曲率。所述利用Weingarten矩阵的特征向量求主方向，是分别用两个特征向量的两个分量作为组合系数，作局部坐标系的X轴和Y轴的线性组合，作为两个主方向。
本发明的有益效果是提出对第三类技术改进的新方法，是在点云数据上直接计算主曲率和主方向。利用点云的坐标信息和点云的法向量信息。本发明与前人方法的区别主要体现在，法向量被直接用确定法截线，通过法截线的法曲率来进行主曲率估计，而不是作为非线性规划问题的约束条件，此外本发明还把非线性最优化问题归结为求解一个线性系统。因此，本发明计算的时间和空间性能也较以往方法优越，实验表明，主方向的计算结果也比其他方法准确。我们利用三维激光扫描仪获取实物的表面信息，并根据扫描数据利用常规方法计算出各个点的法向量，再给出法截线的法曲率估计式，最后经过线性拟合和韦恩伽汀(Weingarten)矩阵计算求出主曲率和主方向。本发明所获得的计算结果能用于计算机图形学各应用领域，包括特征增强、表面简化、网格优化、特征提取、数据压縮、骨架提取和3D实物重建、虚拟现实、电脑游戏、自然场景模拟、城市景观设计、数据压縮等领域具有重要的应用价值。利用本发明，可以快速、方便、准确地计算出点云数据的主曲率和主方向。
图1算法流程2近邻点及其法向量所确定的三角形图3局部空间直角坐标系和法截线示意4圆环面图及其添加噪声之后的效果5圆环面主曲率估计结果的误差比较图6无噪声手模型的主方向估计结果图7有噪声手模型的主方向估计结果图8无噪声大象模型的主方向估计结果图9有噪声大象模型的主方向估计结果图io树枝模型点云数据图11树枝模型点云数据的主方向图12树枝模型点云数据树干表面的主方向图13树枝模型点云数据树干分叉处的主方向图14树枝模型点云数据树干凹陷处的主方向
具体实施例方式
下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。应指出的是，
6所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。
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近邻查找的时间复杂度为0 (log2n)，这里n为点云数据的点的个数。
其次，计算各个点的法向量。对于点云数据的每一个点，利用点云数据的kd树查找15个或30个近邻点，假设这些近邻点来自于同一个平面，于是可以用这些近邻点到拟合平面的残差的绝对值，再乘以权系数的积的和构造最小二乘问题，其中的权的确定是以点云数据中每一个点与近邻点的欧式距离的倒数作为权值。利用最小二乘方法拟合出这个平面，以这个平面的法向量作为点云的法向量的初始估计值，然后通过加权平均算法修正点云的法向量估计，权值等于以点云数据中每一个点与近邻点的欧式距离的倒数。
3逐点计算主曲率和主方向 快速而且准确地计算出扫描得到的点云数据中各个点的主曲率和主方向是本发
明中的主要目标。计算的方法是逐点迭代进行的。计算一个点(记为P)的主曲率和主方
向包括以下七个步骤 3. 1建立局部坐标系 建立局部坐标系，如图3显示了这个局部空间直角坐标系局部三维直角坐标系，设点P的法向量为N = (nx,p， riy,p， nz,p)，则这个p点就是局部坐标系的原点，X， Y， Z三个坐标轴分别为夂=(—sinP，cosp，0) ，F = (cosy cos p， cosy sin p， —sin y)，Z = N = (nx,p，ny,p，nz,p)，其中V = arccos(nz,p) ，p = arctan("w/"w)。
3. 2搜索近邻点 由于预处理阶段建立了 kd树，以点p为查询点，基于欧式距离查找k近邻，比如k=15近邻或k = 30近邻，根据数据扫描精度确定近邻点的个数，扫描精度较低则取15个近邻点，否则取30个近邻点，设定查找近邻的误差阈值为0. O，这里记qi(i = 1，2， ， m)为查询得到的m个近邻点。
3. 3近邻点坐标局部化
7
求出所有近邻点在局部坐标系的坐标。本发明的方法可以描述为下面几个步骤。假设待求主曲率的点为P，近邻点为qi，以定位向量^J,.的坐标作为点&的局部坐标，这里坐标转化的实质是对近邻点做了 一个点的平移变换。 把近邻点的法向量转化为局部坐标系的坐标。本发明的方法是预处理阶段求出的各个近邻点的法向量(记为Mi)在局部坐标系中的方向余弦作为法向量的在局部坐标系中的坐标。 3. 4计算各个近邻点对应的法截线的法曲率 假定待求主曲率和主方向的点、近邻点以及他们的法向量确定了一条法截线，现在估计这条法截线在点P的法曲率。 待求主曲率和主方向的点p、近邻点Qi以及他们的法向量N和Mi确定一个近似三角形，由图2示出的近邻点0、p、qi及其法向量N和Mi所确定的三角形。确定的方法是依据角角边定理，P为近邻点法向量与局部坐标系Z轴的夹角，a为近邻点的定位向量pqi与Z轴夹角的补角，线段pqi为边。在这个近似三角形，依据正弦定理取法曲率为
= - sin 〃 + (| /;仏| sin a)， 其中&表示第i近邻点所对应的法截线的法曲率，I pQi I表示点Qi与点p的欧式距离。为了便于计算，把坐标代入上面式子A！ =-sin" + (l/^,.lsina)，运用三角公式可以把这个式子转化为下面的式子，作为法曲率的估计式
A:=-气+ + :2 . 7< +凡2 ) 其中，、=(A +凡 + V< +凡2 ,=，近邻点的局部坐标为(Xi， yi，Zi)，近邻点&的法向量的局部坐标为(nx,i， riy,i， nz,i)。
3. 5线性回归求出Weingarten矩阵 根据欧拉公式(Euler Equation)构造非线性最优化问题，设需要求点云数据中一个点P的主曲率和主方向g，^ ，则^';n。
g[M。S2W") + *2Sin2(《+ 。-^是构造非线性最优化
问题，其中、，k2是待求的与主方向^，^对应的主曲率，e为未知数，是局部坐标系中X轴与最大主曲率对应的主方向ei的夹角，点Qi是近邻点&的在切平面S上的投影点，9 i为向量pQi与局部坐标系中X轴的夹角。图3示出法截线和主方向和qi、 Qi、 p这些点、角度9 .、 9 。
1 、 o 本发明把上述非线性规划问题々n。 咖2(《+ 0+ ft2 sii2(《+ 。-《f转化为求解下
面的最小二乘问题， min -/ ||
cos2《2cos《sin《sin 《、r 其中，Mmx3 =cos2《2cos《sin《sin 《《
cos2 <9W;Zcos sin
k丄cos2 9 +k2sin2 9 ， B = (k厂k》cos 9 sin 9 ， C = k^sin2 9 +k2cos2 9 这一步的求解可以转化为求解超定线性方程组Mii =R。然后，用求出来的A，B，C的值来构造Weingarten矩阵
8
<formula>formula see original document page 9</formula>
3. 6计算Weingarten矩阵的特征值和特征向量 求出上一步的Weingarten矩阵的特征值和特征向量，对特征值的大小进行排序， 把大的特征值记为、，它也是该点的最大曲率，对应的单位长度的特征向量称为较大特征 值对应的特征向量，另一个较小的特征向量记为IV其值为该点的最小曲率，对应的单位长 度的特征向量称为较小特征值对应的特征向量。 假如该点的两个特征值相等，即表示该点为脐点，而脐点的主方向是平行于切平 面的任意方向，方向不是惟一的。本发明对这种特殊点的处理方法是取标准正交基(i，o) 和(O，l)作为特征向量。 3. 7利用特征值和特征向量求出主曲率和主方向 上一步计算出的特征向量是二维的，是相对于局部坐标系中xOy面上的向量。需 要把它们转化到全局坐标系。转化的方法是把最大曲率对应的特征向量的各个分量作为组 合系数，求出局部坐标系的X轴和Y轴的线性组合作为最大曲率对应的主方向，即为主方向 ei ;同样地，把最小曲率对应的特征向量的各个分量作为组合系数，求出局部坐标系的X轴 和Y轴的线性组合作为最小曲率对应的主方向，即为e" 主曲率的值与坐标系的选取无关，不需要转化，即为上一步计算出来的特征值、 和k2。 4结果输出 对于计算的结果，本发明提供两种形式输出。第一种方式为把点的坐标、法向量、 主方向、主曲率以文件形式存储，再由专业软件进行显示。第二种方式为把主方向、主曲率 信息转化为矢量信息或者颜色信息，结合点云数据显示出来。
5实验结果与结论 我们将此方法应用于解析数据和实际扫描得到的不同复杂度的数据，与之前的两
种主要的主曲率主方向估计方法进行了比较。通过实验证明，本发明对主曲率和主方向的
估计更加准确可靠，鲁棒性也更强。 5. 1解析曲面模拟数据的对比实验 解析数据的实验用来说明本发明提出的新方法比之前的各种方法更加准确可靠。
解析数据来自一个圆环解析曲面，它的解析式是已知的，本发明用的圆环面方程为
r(u，v) = ((R+rcosu) cos v， (R+rcosu)sin v， rsi皿)
(0《u《2ji，0《v《2ji) 其中参数R = 2，r = 1。从这个解析曲面，按照二维均匀分布获取5000个随机采 样点，根据微分几何公式可以计算出各个点的主曲率和主方向(称为采样点的实际主曲率 和实际主方向)。然后，把这5000个采样点视作扫描得到的点云数据。对每一个点p，沿着 法向量方向添加噪声(p +^ )，这里^为点p的法向量，a服从[-mh， mh]的均匀分布，m是 任意两点之间的距离的中位数，h为噪声水平，实验中用过的噪声水平为h二O. 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平均曲率估计绝对误差对比表
噪声 水平高斯曲率误差平均曲率误差Goldfeath方法本发明的方法Goldfeath方法本发明的方法
0. 20.0971380.0883880. 1416800.192398
0,40. 1224850. 1085300.1794230. 206980
0.60. 1591270.1395150. 2469970.238619
0.80.1970230.1673080. 3033470.272985
1.00.2644660. 2156750. 4522470.343714 从表1可以看出，本发明与Go 1 df eath方法相比，除了在噪声水平较低时平均曲率 误差略大以外，其余情况误差更小，说明计算更加准确。对于噪声较大的情况，本发明优势 更加明显。 5. 2点云数据的对比实验 把本发明与Goldfeath方法应用于实际扫描数据，用了两个模型进行测试。
第一个模型是一个"手"的数据。该模型的高度为202个长度单位，有6258个随机 采样点，数据本身带有较准确的法向量。图6为无噪声手模型的主方向估计结果，显示初始 数据图以及本发明与Goldfeath方法的主方向估计结果。给数据点沿着法向量方向加入噪 声，噪声水平为2.0个单位，则法向量与点的位置信息不是完全正确，求出的主方向显示在 图7为有噪声手模型的主方向估计结果。图6(a)为原始数据，图6(b)为本发明的实验计算 结果，图6 (c)为Goldfeath方法的计算结果；图7(a)为原始数据；图7(b)为本发明的计算
10结果，图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为树枝模型原始的点云数据。图ll为 树枝模型点云数据各个点的主方向图。图12显示了树枝模型点云数据树干表面各个点的 主方向，其特征是主曲率较大者对应的主方向体现了树干的横截面方向，而主曲率较小者 对应的主方向刻画了树干的生长方向。图13显示了树枝模型点云数据树枝分叉点处的主 方向特征。图14显示了树枝模型点云数据树干表面凹陷处的主方向特征，这些主曲率较大 者对应的主方向指向凹陷的中心。 本发明提出的主曲率和主方向计算方法的特征在于利用点云数据和各个点的法 向量计算点的主曲率和主方向。 上述实验结果和利用点云数据求取主方向和主曲率的方法，可以用于计算机图形 学各应用领域，具有高可信度、操作简单、应用前景广的特点。 以上所述，仅为本发明中的具体实施方式
，但本发明的保护范围并不局限于此，任 何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在 本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。
权利要求
一种点云数据的主曲率和主方向估计方法，其特征在于，该主曲率和主方向估计步骤包括步骤1利用激光扫描仪扫描直接采集点云数据并对点云数据预处理，按照点云数据中每个点的坐标进行空间划分，实现三维空间的二分查找树的数据存储结构称为kd树(k-dimensional tree)；步骤2对于点云数据的每一个点，利用点云数据的kd树查找15个或30个近邻点，根据最小二乘方法把这些近邻点拟合出一个平面，以这个平面的法向量作为点云法向量的初始估计值，然后通过加权平均算法修正点云数据的各个点的法向量估计；步骤3对于点云数据的每一个点，利用其法向量、切平面构造局部三维直角坐标系；步骤4对于点云数据的每一个点，利用点云数据的kd树查找15个或30个近邻点；步骤5对于查找到的近邻点，通过三维坐标变换，把这些近邻点的原始坐标和这些近邻点的法向量都转化为局部坐标系的坐标；步骤6利用点云数据的每一个点及其法向量、一个近邻点、近邻点的法向量构造近似三角形，根据正弦定理给出点云的法截线的法曲率的近似表达式；步骤7在局部坐标系中，利用法曲率，根据欧拉公式(Euler Equation)构造非线性最优化问题。通过三角形公式进行恒等变换，把这个非线性最优化问题转化为线性拟合，求出韦恩伽汀矩阵(Weingarten矩阵)的各个元素；步骤8利用矩阵的奇异值分解(SVD分解)求出Weingarten矩阵的特征值和特征向量；步骤9利用Weingarten矩阵的特征值和特征向量求出主曲率和主方向。
2. 按权利要求1所述的方法，其特征在于，所述每个点的坐标，都是采用激光扫描仪产生的原始坐标。
3. 按权利要求1所述的方法，其特征在于，所述加权平均算法修正点云数据的各个点的法向量估计，是以点云数据中每一个点与近邻点的欧式距离的倒数作为权值，所述欧式距离的倒数，若距离值为0，则不计算这个近邻点。
4. 按权利要求1所述的方法，其特征在于，所述根据最小二乘方法把这些近邻点拟合出一个平面，其中拟合平面时需要构造带权最小二乘问题，这个最小二乘问题是这些近邻点到拟合平面的残差的绝对值，再乘以权系数的积的和的最小值问题。
5. 按权利要求4所述的方法，其特征在于，所述权系数，是以点云数据中每一个点与近邻点的欧式距离的倒数作为权值。
6. 按权利要求l所述的方法，其特征在于，所述局部三维直角坐标系，其构造方法是，对于点云数据中一个点P，若点P的法向量为N二 (nx,p， riy,p， n,.p)，则这个点P就是局部坐标系的原点，X， Y， Z三个坐标轴分别为X = (-sinp， cosp， 0)，y = (cosy cos p， cos"—，-sin y) ，Z = N = (nx,p， ny,p， nz, p)，其中p = arcta—J"")表示法向量N的第一方向角，v = arccos(nz,p)表示法向量的第三方向角的余角。
7. 按权利要求1所述的方法，其特征在于，所述把这些近邻点的原始坐标和这些近邻点的法向量都转化为局部坐标系的坐标，是把这些近邻点的坐标减去局部坐标系原点的坐标就是这些近邻点在局部坐标系中的坐标，这些近邻点的法向量坐标与局部坐标系的三个坐标轴分别作数量积，算出这些近邻点的法向量在局部坐标系中的坐标。
8. 按权利要求1所述的方法，其特征在于，所述构造近似三角形，该三角形由角角边定理确定设点云中的待计算法曲率的点为P，它的一个近邻点为qi点P和qi的法向量分别为N和Mi，向量N和Mi的夹角为13 ，向量N和向量^,.的夹角为a ，线段|pq」为边，则角|3 ，角a ，边|pq」确定了这个三角形。
9. 按权利要求1所述的方法，其特征在于，所述用正弦定理给出点云的法截线的法曲率的近似表达式为<formula>formula see original document page 3</formula>其中表示第i近邻点所对应的法截线的法曲率，P为近邻点法向量与局部坐标系Z轴的夹角，a为近邻点的定位向量与Z轴夹角的补角，局部坐标系的原点为P，近邻点为qi， |pq」表示点qi与点P的欧式距离。
10. 按权利要求9所述的方法，其特征在于，所述近似表达式<formula>formula see original document page 3</formula>在实际计算中采用近似计算形式<formula>formula see original document page 3</formula>其中，<formula>formula see original document page 3</formula>近令卩点qi的局部坐标为(Xi，yi，Zi)，近邻点qi的法向量的局部坐标为(r^,i，ny,i，r^i)。
11. 按权利要求1所述的方法，其特征在于，所述根据欧拉公式(EulerEquation)构造非线性最优化问题，设需要求点云数据中 一 个点的主曲率和主方向，则<formula>formula see original document page 3</formula>是构造的非线性最优化问题，其中、，k2是待求的主曲 率，未知数e为局部坐标系中x轴与最大主曲率对应的主方向的夹角，9i为近邻点qi在局部坐标系中x0y面的投影点的定位向量与局部坐标系中X轴的夹角。
12. 按权利要求1所述的方法，其特征在于，所述把这个非线性最优化问题转化为线性拟合，就是把所有近邻点的法曲率作为因变量观测值，所有近邻点与局部坐标系的x轴的方向角的余弦值的平方、方向角的2倍的正弦值、方向角的正弦值的平方作为自变量的观测值，作三元线性拟合。三个拟合系数值依次作为二阶对称矩阵Weingarten矩阵的上三角的三个元素。
13. 按权利要求1所述的方法，其特征在于，所述求出Weingarten矩阵的特征值和特征向量，若两个特征值相等，则该点为脐点，特征向量取为(l，O)和(O，l)。
14. 按权利要求l所述的方法，其特征在于，所述利用Weingarten矩阵的特征值求出主曲率，是用两个特征值作为主曲率。
15. 按权利要求l所述的方法，其特征在于，所述利用Weingarten矩阵的特征向量求主方向，是分别用两个特征向量的两个分量作为组合系数，作局部坐标系的X轴和Y轴的线性组合，作为两个主方向。
全文摘要
本发明涉及一种点云数据的主曲率和主方向估计方法步骤包括预处理，法截线曲率估计，拟合出Weingarten矩阵，求出Weingarten矩阵的特征值和特征向量，求出主曲率和主方向。本发明仅利用激光扫描仪的扫描数据和估计的法向量，得到忠实于原始实物的主曲率和主方向。该方法通过最小二乘线性拟合和矩阵的特征值特征向量求出主曲率和主方向，算法简单，计算结果准确，时间复杂度是高效的。该方法称为法截线拟合法，其计算结果在虚拟现实、电脑游戏、自然场景模拟、城市景观设计、数据压缩、特征提取、实物3D重建等领域具有重要的应用价值。
文档编号G01B11/24GK101751695SQ20081023932
公开日2010年6月23日 申请日期2008年12月10日 优先权日2008年12月10日
发明者张晓鹏, 李红军, 程章林 申请人:中国科学院自动化研究所