CN102306397A - 点云数据网格化的方法 - Google Patents

Abstract

本发明公开了一种点云数据网格化的方法。该方法针对由激光扫描得到的离散点云数据，提供一个自动的、鲁棒性强的全局参数化方法，并利用参数化结果直接获取与主方向一致的能够反映模型内在几何特征的网格化结果。网格化结果有两种形式：全部由四边形组成，或者全部由三角形组成。本发明仅利用模型的点云数据而无需事先三角网格化进行点云数据网格化，处理过程完全自动，可以处理含噪声点云，无需人工干预，并且通过参数可以快速的控制四边形网格化的疏密程度，得到各种分辨率的四边形网格、或者三角形网格。

Description

点云数据网格化的方法

技术领域

本发明涉及计算机图形学和计算机视觉技术领域，特别涉及一种点云数据网格化的方法。

背景技术

由于快速而精确的激光扫描仪的发展，点云数据已经在计算机辅助设计和计算机图形学领域得到广泛的应用。通常原始点云数据不包含任何拓扑结构信息，因此大量的研究工作集中在如何从点云数据重建网格表面。但是大部分已有的工作只关注如何产生高质量的三角面片网格模型，对于三角形面片的形状和方向缺乏控制。由于四边形网格张量积特性，相对于三角形网格，四边形网格在很多领域都具有优势，例如B样条拟合、纹理映射等等。尤其是方向和主方向一致的四边形在建模时更具有优势，因为它们能够反映几何模型的对称性。

全局参数化是一种有效解决四边形网格化的方法。Bommes(BOMMES，D.，ZIMMER，H.，AND KOBBELT，L.2009.Mixed-integerquadrangulation.ACM Transactions on Graphics 28，3，1-10.)将三角形的四边形网格化问题转化为一个混合整数的二次优化问题，该方法能够生成完全的四边形网格而且很好的保持与主方向的一致性。

但是，申请人意识到现有技术存在如下技术缺陷：上述全局参数化方法仅能用于三角形网格，对于初始点云数据由于缺少点与点之间的拓扑连接关系，直接利用全局参数化方法对点云数据处理具有一定的难度。

发明内容

(一)要解决的技术问题

为解决上述缺陷，本发明提供了一种点云数据网格化的方法，以利用模型的点云数据而无需事先三角网格化来进行点云数据的四边形网格化。

(二)技术方案

根据本发明的一个方面，提供了一种点云数据网格化的方法。该点云数据网格化的方法包括：步骤A，获取点云数据；步骤B，由点云数据计算出每一点的主方向、法向量及点云数据的奇异点信息；步骤C，由点云数据每一点的坐标，该点的法向量、主方向及奇异点信息，获得每一点的标量函数参数值θ和

该标量函数满足在每一点的梯度与该点的主方向保持一致，其中奇异点处的标量函数参数值被约束为整数值；步骤D，由标量函数参数值θ和

建立与点云数据形状一致的三角形网格。

优选地，本发明点云数据网格化的方法中，步骤B包括：步骤B1；对初始点云数据获取每一点的法向量；步骤B2，由法向量计算该点的曲率张量，获取该点初始主方向；步骤B3，通过使相邻点的初始主方向的一致最大化，对每一点的初始主方向进行平滑，获取每一点的主方向；步骤B4，根据点云数据中每一点的主方向，确定点云数据中的奇异点。

优选地，本发明点云数据网格化的方法中，步骤B1包括：步骤B1a，对于点云数据的每一个点p，利用点云数据的kd树查找点p的n个近邻点；步骤B1b，假设这些近邻点来自于同一个平面T，用这些近邻点到拟合平面的残差的绝对值，再乘以权系数的积的和构造最小二乘问题，其中的权的确定是以点云数据中每一个点与近邻点的欧式距离的倒数作为权值；步骤B1c，利用最小二乘方法拟合出平面T，以这个平面的法向量作为点云中该点的法向量

优选地，n为整数，优选取15或30。

优选地，本发明点云数据网格化的方法中，步骤B2包括：步骤B2a，对每一点p建立局部坐标系，设点p的法向量为则这个p点就是局部坐标系的原点，设局部坐标系的三个方向分别为

步骤B2b，对于p的一个近邻点p_i，其法向量为

则点p处的曲率张量必须满足下述方程式的约束：

$\left[\begin{array}{ccc}{\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{u}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{u}& \stackrel{\→}{N}\·\stackrel{\→}{u}\\ {\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{v}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{v}& \stackrel{\→}{N}\·\stackrel{\→}{v}\\ {\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{w}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{w}& \stackrel{\→}{N}\·\stackrel{\→}{w}\end{array}\right]\·\left[\begin{array}{c}{\mathrm{\Δ}}_{\stackrel{\→}{p}}\·\stackrel{\→}{u}\\ {\mathrm{\Δ}}_{\stackrel{\→}{p}}\·\stackrel{\→}{v}\\ 1\end{array}\right]=\left[\begin{array}{c}\stackrel{\→}{N_i}\·\stackrel{\→}{u}\\ \stackrel{\→}{N_i}\·\stackrel{\→}{v}\\ \stackrel{\→}{N_i}\·\stackrel{\→}{w}\end{array}\right];$

步骤B2c，对点p的m个近邻点的法方向进行拟合，代入上述方程中，求解出左边矩阵，即为该点处的曲率张量，矩阵的最大两个特征值所对应的特征向量即为对应的点p的初始主方向

和

优选地，m为15。

优选地，本发明点云数据网格化的方法中，步骤B3包括：步骤B3a，定义衡量点云数据相邻点的主方向差异的函数

$E_{\mathrm{smoothing}}=\underset{e_{\mathrm{ij}}}{\mathrm{\Σ}}{(({\mathrm{\β}}_i-{\mathrm{\β}}_j)+{\mathrm{\α}}_i-{\mathrm{\α}}_j+\mathrm{n\π}/2)}^2,$

其中，α_i为所要求的主方向与一个参考方向的夹角，该参考方向取为点i的切平面上的任意一个方向，β_i为原始主方向K′与该参考方向的夹角，e_ij为k近邻图中的一条边，n为整数决定两个相邻点的主方向如何变化；步骤B3b，为了求解上述函数，首先通过近似变化来消除整数变量n，引入变量sin4α_i和cos4α_i来代替原始变量α_i，方程变为：

$E_{\mathrm{smoothing}}=\underset{e_{\mathrm{ij}}}{\mathrm{\Σ}}{\{\left(\begin{array}{c}\cos 4{\mathrm{\α}}_i\\ \sin {4\mathrm{\α}}_i\end{array}\right)-\left(\begin{array}{cc}\cos \mathrm{\β}& -\sin \\ \sin \mathrm{\β}& \cos \mathrm{\β}\end{array}\right)\left(\begin{array}{c}\cos {4\mathrm{\α}}_j\\ \sin {4\mathrm{\α}}_j\end{array}\right)\}}^2;$

其中β＝β_i-β_j，对sin4α_i和cos4α_i求导得到线性方程，求解得到sin4α_i和cos4α_i的解，通过sin4α_i和cos4α_i最终求出原始变量α_i；步骤B3c，利用夹角α_i，获取点云数据中每一点的主方向。

优选地，本发明点云数据网格化的方法中，步骤B4包括：步骤B4a，将点p的近邻点和相对应的主方向投影到点p的切平面上；步骤B4b，在切平面上对这些点按照逆时针顺序排序；步骤B4c，选择其中一个近邻点在切平面上的主方向作为参考方向，从该参考方向出发，根据上述排序后的顺序求出每两个相邻点的主方向的变化角度；步骤B4d，若所有变化角度之和位于区间[-π/2，π/2]之外，则该点被标记为奇异点。

优选地，本发明点云数据网格化的方法中，步骤C包括：步骤C1，对点云数据进行切割，转化为亏格为0，边界为一的拓扑结构；步骤C2，将k近邻图中的每一条边投影到分割线及其法向量所定义的平面上，此处分割线为步骤C1中切割所产生的分割线，若投影后的边与分割线相交，则在k近邻图中去除此边，从而得到一个新的k近邻图；步骤C3，对于定义在点云上每一个点的两个标量函数θ和

求出θ和

在每一个点的具体值，使得该点的最大主方向和最小主方向分别与两个标量函数的梯度尽量保持一致，两个标量函数θ和

为下述方程最小值的最优解：

其中，

和

分别为两个函数在每点的梯度，

和分别为该点的最大主方向和最小主方向，ω为由用户指定的控制参数化疏密程度的参数。

优选地，本发明点云数据网格化的方法中，步骤C1包括：步骤C1a，利用点云数据和莫尔斯函数计算出点云的同源基，沿着同源基对点云进行分割；步骤C2b，在k近邻图上求取每个奇异点到边界的最短路径，并按照连接路径对点云进行再次分割，获取由点云中若干点连接而成的路径。

优选地，本发明点云数据网格化的方法中，步骤C3包括：步骤C3a，定义衡量标量函数在每一点的梯度与该点的最大主方向和最小主方向之间差异的能量函数：

公式1

其中，θ_i，

分别表示点i对应的两个标量函数的值，K，K^⊥分别表示点i处的最大主方向和最小主方向，ij为k近邻图中的边，e_ij表示连接点i，j的线段的矢量，w是用户指定的控制参数化疏密程度的参数；步骤C3b，对于位于分割线上的点P，其近邻点被分割线分成两类，该点也相应的被分成两个点P⁺和P^-，以及对应的主方向K⁺和K^-。为了保证参数值θ和在分割线处的连续性，点P处的参数值必须满足以下条件：

公式2，

其中，(j，k)为整数变量，

表示方向K^-到K⁺所需要的旋转矩阵；并且，奇异点处的参数值也被限定为整数值；步骤C3c，通过公式1和公式2建立一个在分割线处连续的全局参数化方程；步骤C3d，求解该全局参数化方程，获得每一点的标量函数参数值θ和

优选地，本发明点云数据网格化的方法中，步骤D包括：步骤D1，对于每一个点及其周围k个近邻点根据参数化值

进行狄洛尼三角化；

步骤D2，取与该点相邻的三角形来建立三维空间中的对应点的连接关系；

步骤D3，如果有两个以上的三角形共有同一条边，则逐步移除多余的三角形直至每一条边最多被两个三角形共有；步骤D4，如果移除冗余三角形的过程中产生空洞，则通过对该洞进行狄洛尼三角化来填充该洞。优选地，k＝20。

优选地，本发明点云数据网格化的方法中，该方法还包括：步骤E，在建立的三角网格基础上建立四边形网格。优选地，步骤E包括：步骤E1，求取每个三角形内的等值线段；步骤E2，求取等值线段之间的交点；步骤E3，对于每一个等值线段之间的交点，按照等值线段之间的连接关系找到该交点的邻接点并连接，建立形状均匀且方向符合主方向的四边形网格。

(三)有益效果

本发明直接在点云模型上进行全局参数化和四边形网格化，仅利用模型的点云数据而无需事先三角网格化可以得到同样质量的结果，并且处理过程完全自动，无需人工干预，并且通过参数可以快速的控制四边形网格化的疏密程度，得到各种分辨率的四边形网格。

附图说明

图1为点云数据网格化方法的流程图；

图2为本发明实施例点云数据网格化方法中主方向场计算的流程图；

图3为本发明实例点云数据网格化方法中奇异点计算的流程图；

图4为本发明实施例点云数据网格化方法中全局参数化流程图；

图5为本发明实施例点云数据网格化方法中四边形网格化流程图；

图6a为本发明实施例点云数据网格化方法中重叠三角形处理示意图；

图6b为图6a中的重叠三角形被去除后的结果；

图6c为图6b中由于去除重叠三角形所产生的洞被补齐后的结果；

图7a为茶杯模型的原始点云数据；

图7b为茶杯模型的参数化结果；

图7c为根据参数化结果得到的茶杯模型的四边形网格化结果；

图8a为采用Bommes方法得到的四边形网格化结果；

图8b为图7a中四边形网格化结果的反映四边形网格质量的直方图；

图8c为采用本发明方法得到的四边形网格化结果；

图8d为图7c中四边形网格化结果的反映四边形网格质量的直方图；

图9a为零件1模型的原始点云数据；

图9b为零件1模型的参数化结果；

图9c为根据参数化结果得到的零件1模型的四边形网格化结果；

图9d为零件2模型的原始点云数据；

图9e为零件2模型的参数化结果；

图9f为根据参数化结果得到的零件2模型的四边形网格化结果。

图10a为采用Floater方法得到的参数化结果的二维显示；

图10b为本发明方法得到的参数化结果的二维显示；

图10c为Floater方法的三角网格重建结果；

图10d为本发明方法的三角网格重建结果。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

在本发明的一示例性实施例中，公开了一种点云数据网格化方法。该点云数据网格化包括：步骤A，获取点云数据；步骤B，由点云数据计算出每一点的法向量、主方向及点云数据中的奇异点信息；步骤C，由点云数据每一点的坐标，该点的主方向，法向量以及奇异点信息，获得每一点的标量函数参数值θ和

该标量函数满足在每一点的梯度与该点的主方向尽量保持一致，即下述方程最小值的最优解：

其中，

和分别为两个函数在每点的梯度，

和

分别为该点的最大主方向和最小主方向，ω为由用户指定的控制参数化疏密程度的参数，其中奇异点处的标量函数参数值被约束为整数值。步骤D，由标量函数参数值θ和

建立与所述点云数据形状一致的三角形网格，网格中每个三角形所有内角的大小近似。优选地，在上述步骤D之后还可以包括：步骤E，在建立的三角网格基础上建立四边形网格，网格中每个四边形接近于正方形。

对于其中的步骤B而言，具体包括：步骤B1；对初始点云数据获取每一点的法向量；步骤B2，由法向量计算该点的曲率张量，进而获取该点初始主方向；步骤B3，通过使相邻点的初始主方向的一致最大化，对每一点的初始主方向进行平滑，获取每一点的主方向；步骤B4，根据点云数据中每一点的主方向，确定点云数据中的奇异点。

其中，步骤B1具体包括：步骤B1a，对于点云数据的每一个点p，利用点云数据的kd树查找点p的n个近邻点；步骤B1b，假设这些近邻点来自于同一个平面T，用这些近邻点到拟合平面的残差的绝对值，再乘以权系数的积的和构造最小二乘问题，其中的权的确定是以点云数据中每一个点与近邻点的欧式距离的倒数作为权值；步骤B1c，利用最小二乘方法拟合出平面T，以这个平面的法向量作为点云中该点的法向量

。优选地，n为整数，可以取为15或30。

其中，步骤B2具体包括：步骤B2a，对每一点p建立局部坐标系，设点p的法向量为

则这个p点就是局部坐标系的原点，设局部坐标系的三个方向分别为步骤B2b，对于p的一个近邻点pi，其法向量为

则点p处的曲率张量必须满足下述方程式的约束：

$\left[\begin{array}{ccc}{\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{u}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{u}& \stackrel{\→}{N}\·\stackrel{\→}{u}\\ {\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{v}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{v}& \stackrel{\→}{N}\·\stackrel{\→}{v}\\ {\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{w}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{w}& \stackrel{\→}{N}\·\stackrel{\→}{w}\end{array}\right]\·\left[\begin{array}{c}{\mathrm{\Δ}}_{\stackrel{\→}{p}}\·\stackrel{\→}{u}\\ {\mathrm{\Δ}}_{\stackrel{\→}{p}}\·\stackrel{\→}{v}\\ 1\end{array}\right]=\left[\begin{array}{c}\stackrel{\→}{N_i}\·\stackrel{\→}{u}\\ \stackrel{\→}{N_i}\·\stackrel{\→}{v}\\ \stackrel{\→}{N_i}\·\stackrel{\→}{w}\end{array}\right];$

步骤B2c，对点p的m个近邻点的法方向进行拟合，代入上述方程中，求解出左边矩阵，即为该点处的曲率张量，矩阵的最大两个特征值所对应的特征向量即为对应的点p的初始主方向

和

优选地，m为整数，可以取为15。

其中，步骤B3具体包括：步骤B3a，定义衡量点云数据相邻点的主方向差异的函数：

$E_{\mathrm{smoothing}}=\underset{e_{\mathrm{ij}}}{\mathrm{\Σ}}{(({\mathrm{\β}}_i-{\mathrm{\β}}_j)+{\mathrm{\α}}_i-{\mathrm{\α}}_j+\mathrm{n\π}/2)}^2,$

其中，α_i为所要求的主方向与一个参考方向的夹角，该参考方向可以取为点i的切平面上的任意一个方向，β_i为原始主方向K′与该参考方向的夹角，e_ij为k近邻图中的一条边，n为整数决定两个相邻点的主方向如何变化；步骤B3b，为了求解上述函数，首先通过近似变化来消除整数变量n，引入变量sin4α_i和cos4α_i来代替原始变量α_i，方程变为

$E_{\mathrm{smoothing}}=\underset{e_{\mathrm{ij}}}{\mathrm{\Σ}}{\{\left(\begin{array}{c}\cos 4{\mathrm{\α}}_i\\ \sin {4\mathrm{\α}}_i\end{array}\right)-\left(\begin{array}{cc}\cos \mathrm{\β}& -\sin \\ \sin \mathrm{\β}& \cos \mathrm{\β}\end{array}\right)\left(\begin{array}{c}\cos {4\mathrm{\α}}_j\\ \sin {4\mathrm{\α}}_j\end{array}\right)\}}^2$

其中，β＝β_i-β_j，上述方程为二次优化问题，对sin 4α_i和cos 4α_i求导即可得到线性方程，求解即可得到sin4α_i和cos4α_i的解，通过sin4α_i和cos4α_i可以最终求出原始变量α_i；步骤B3c，利用最优方向即夹角α_i，获取点云数据中每一点的主方向。

其中，步骤B4具体包括：步骤B4a，将点p的近邻点和相对应的主方向投影到点p的切平面上；步骤B4b，在切平面上对这些点按照逆时针顺序排序；步骤B4c，选择其中一个近邻点在切平面上的主方向作为参考方向，从该参考方向出发，根据上述排序后的顺序求出每两个相邻点的主方向的变化角度；步骤B4d，若所有变化角度之和位于区间[-π/2，π/2]之外，则该点被标记为奇异点。

对于其中的步骤C而言，具体包括：步骤C1，对点云数据进行切割，转化为亏格为0，边界为一的拓扑结构；步骤C2，将k近邻图中的每一条边投影到分割线及其法向量所定义的平面上，此处分割线为步骤C1中切割所产生的分割线，若投影后的边与分割线相交，则在k近邻图中去除此边，从而得到一个新的k近邻图；步骤C3，对于定义在点云上每一个点的两个标量函数θ和

求出θ和

在每一个点的具体值，使得该点的最大主方向和最小主方向分别与两个标量函数的梯度尽量保持一致。

其中，步骤C1具体包括：步骤C1a，利用点云数据和莫尔斯函数计算出点云的同源基，沿着同源基对点云进行分割；步骤C2b，在k近邻图上求取每个奇异点到边界的最短路径，并按照连接路径对点云进行再次分割，获取由点云中若干点连接而成的路径；

其中，步骤C3具体包括：步骤C3a，定义衡量标量函数在每一点的梯度与该点的最大主方向和最小主方向之间差异的能量函数：

其中，θ_i，

分别表示点i对应的两个标量函数的值，K，K^⊥分别表示点i处的最大主方向和最小主方向，ij为k近邻图中的边，e_ij表示连接点i，j的线段的矢量，w是用户指定的控制参数化疏密程度的参数；

步骤C3b，对于位于分割线上的点P，其近邻点被分割线分成两类，该点也相应的被分成两个点P⁺和P^-，以及对应的主方向K⁺和K^-。为了保证参数值θ和

在分割线处的连续性，点P处的参数值必须满足以下条件：

其中，(j，k)为整数变量，表示方向K^-到K⁺所需要的旋转矩阵；此外，奇异点处的参数值也被限定为整数值；步骤C3c，通过公式(2)和公式(3)可以建立一个在分割线处连续的全局参数化方程；步骤C3d，求解该全局参数化方程，获得每一点的标量函数参数值θ和

对于其中的步骤D而言，具体包括：步骤D1，对于每一个点及其周围k个近邻点(k为整数，可以取k＝20)根据参数化值

进行狄洛尼三角化；步骤D2，取与该点相邻的三角形来建立三维空间中的对应点的连接关系；步骤D3，如果有两个以上的三角形共有同一条边，则逐步移除多余的三角形直至每一条边最多被两个三角形共有；步骤D4，如果移除冗余三角形的过程中产生空洞，则通过对该洞进行狄洛尼三角化来填充该洞；步骤D5，由三角形网格获得四边形网格。

其中，步骤D5具体包括：步骤D5a，求取每个三角形内的等值线段；步骤D5b，求取等值线段之间的交点；步骤D5c，对于每一个等值线段之间的交点，按照等值线段之间的连接关系找到该交点的邻接点并连接，建立形状均匀且方向符合主方向的四边形网格。

以下将在上述实施例的基础上，给出本发明的最优实施例。需要说明的，该最优的实施例仅用于理解本发明，并不用于限制本发明的保护范围。并且，最优实施例中的特征，在无特别注明的情况下，均同时适用于方法实施例和装置实施例，在相同或不同实施例中出现的技术特征在不相互冲突的情况下可以组合使用。

图1为点云数据网格化方法的流程图。如图1所示，点云数据网格化包括以下三个基本步骤：主方向场的计算，即对点云中的每一个点计算出主方向；全局参数化；四边形网格化。针对每个步骤的具体算法以下将详细说明。

图2为本发明实施例点云数据网格化方法中主方向场计算的流程图。如图2所示，主方向场的计算首先需要求取点云每一点的法向量。因为三维点云数据一般只有点的坐标信息。为求取点云的曲率张量信息，以及对点云进行局部三角化，获得点云数据中的每一个点的法方向是必须的。首先，建立kd树。在计算几何中，kd树是已经被证明的查找近邻的最快捷的数据结构之一。kd树基于点的空间位置信息，通过二分法迭代划分三维空间，实现最优存储。在kd树上，进行k近邻查找的时间复杂度为O(log2n)，这里n为点云数据的点的个数。

为求取各个点的法向量，对于点云数据的每一个点p，利用点云数据的kd树查找n个近邻点(n为整数，可以取为15或30)，假设这些近邻点来自于同一个平面T，于是可以用这些近邻点到拟合平面的残差的绝对值，再乘以权系数的积的和构造最小二乘问题，其中的权的确定是以点云数据中每一个点与近邻点的欧式距离的倒数作为权值。利用最小二乘方法拟合出平面T，以这个平面的法向量作为点云中该点的法向量

由于本实施例方法的全局参数化是受点云数据的主方向约束的，因此获取平滑且准确的主方向场是必要的。主方向场的建立包括两个步骤，如图2后两步所示，即计算各个点的曲率张量与主方向场的全局平滑。

为计算各点的曲率张量，首先对每一点建立局部坐标系，设点p的法向量为

则这个p点就是局部坐标系的原点，设局部坐标系的三个方向分别为

对于p的一个近邻点p_i，上述步骤所得法向量为则点p处的曲率张量必须满足下述方程式的约束：

$\left[\begin{array}{ccc}{\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{u}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{u}& \stackrel{\→}{N}\·\stackrel{\→}{u}\\ {\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{v}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{v}& \stackrel{\→}{N}\·\stackrel{\→}{v}\\ {\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{w}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{w}& \stackrel{\→}{N}\·\stackrel{\→}{w}\end{array}\right]\·\left[\begin{array}{c}{\mathrm{\Δ}}_{\stackrel{\→}{p}}\·\stackrel{\→}{u}\\ {\mathrm{\Δ}}_{\stackrel{\→}{p}}\·\stackrel{\→}{v}\\ 1\end{array}\right]=\left[\begin{array}{c}\stackrel{\→}{N_i}\·\stackrel{\→}{u}\\ \stackrel{\→}{N_i}\·\stackrel{\→}{v}\\ \stackrel{\→}{N_i}\·\stackrel{\→}{w}\end{array}\right]$

在本方法中我们对点p的若干个(15个)近邻点的法方向进行拟合，代入上述方程中，解出左边矩阵里的张量项，即可求出对应的点p的初始主方向K′。主方向在全局参数化过程中起约束作用。

主方向场的全局平滑。平滑过程使相邻点的主方向尽可能一致。本方法定义了一个衡量点云数据相邻点之间主方向差异的函数，通过求解该函数的最小值，可以得到平滑后的主方向。

该函数的定义如下：

$E_{\mathrm{smoothing}}=\underset{e_{\mathrm{ij}}}{\mathrm{\Σ}}{(({\mathrm{\β}}_i-{\mathrm{\β}}_j)+{\mathrm{\α}}_i-{\mathrm{\α}}_j+\mathrm{n\π}/2)}^2$

其中，α_i为所要求的主方向与一个参考方向的夹角，该参考方向可以取为点i的切平面上的任意一个方向，β_i为原始主方向K′与该参考方向的夹角，e_ij为k近邻图中的一条边。该方程可以转化为一个二次优化问题，使用最速下降法解出优化问题，最优解即平滑后的主方向与参考方向的夹角，利用此夹角即可求出平滑后的主方向K。

获取主方向场之后需要确定奇异点的位置。图3是奇异点计算的流程图。首先将点p的近邻点和相对应的主方向投影到点p的切平面上，然后在切平面上对这些点按照逆时针顺序排序。选择其中一个近邻点在切平面上的主方向作为参考方向，从该参考方向出发，根据上述排序后的顺序求出每两个相邻点的主方向的变化角度，若所有变化角度之和位于区间[-π/2，π/2]之外则该点被标记为奇异点。在全局参数化的过程中，奇异点处的参数化值被约束为整数值。

为了对任意拓扑形状的点云进行全局参数化，需要对点云进行切割，转化为亏格为0，边界为一的拓扑结构。首先利用莫尔斯函数计算出点云的同源基，沿着同源基对点云进行分割，然后在k近邻图上求取每个奇异点到边界的最短路径，并按照连接路径对点云进行再次分割。该步骤输入的为原始点云，输出数据为由点云中若干点连接而成的路径，在全局参数化之前需要根据该路径对k近邻图进行更新。

图4是本发明实施例点云数据网格化方法中全局参数化流程图。全局参数化的目的是对于定义在点云上每一个点的两个标量函数θ和

求出θ和

在每一个点的具体值，使得该点的最大主方向和最小主方向分别与两个标量函数的梯度尽量保持一致。输入数据为点云每一点的坐标，该点的主方向，法向量以及k近邻图，输出数据为θ和

在每一个点的具体值。为了计算点云的能量函数，本方法首先定义一种能量函数来衡量这种差异，定义如下：

在实际求取中，上述方程离散化为下面的形式：

其中，θ_i，

分别表示点i对应的两个标量函数的值，K，K^⊥分别表示点i处的最大主方向和最小主方向，ij为k近邻图中的边，e_ij表示连接点i，j的线段的矢量，w是用户指定的控制参数化疏密程度的参数。

由于对点云进行了分割，因此在定义点云上的标量函数的梯度的时候必须去除k近邻图中与分割线相交的边，从而保持全局参数化的连续性。对于每一个靠近分割线的点，检测k近邻图中与该点相连的边是否与分割线相交，即将该边投影到分割线及其法向量所定义的平面上，若投影后的边与分割线相交，则在k近邻图中去除此边，从而得到一个新的k近邻图。

利用上述方法得到的k近邻图的连接关系，根据公式(2)建立能量函数。全局参数化则转化为求取能量函数的最优解。为了保证全局参数化结果在分割处的连续性，在参数化的过程中需要添加额外的约束条件来保证参数值在跨越分割线两边时保持一致。对于位于分割线上的点P，其近邻点被分割线分成两类，该点也相应的被分成两个点P⁺和P^-，以及对应的主方向K⁺和K^-。为了保证参数值θ和在分割线处的连续性，点P处的参数值必须满足以下条件：

其中，(j，k)为整数变量，

表示方向K^-到K⁺所需要的旋转矩阵。

另外，奇异点处的参数值也被限定为整数值。通过公式(2)和公式(3)可以建立一个在分割线处连续的全局参数化方程，该方程的优化是一个混合整数的优化问题，利用文献(BOMMES，D.，ZIMMER，H.，ANDKOBBELT，L.2009.Mixed-integer quadrangulation.ACM Transactions onGraphics 28，3，1-10.)中提供的求解器可以快速而准确的求解。该求解器首先利用高斯消去法消去公式(3)中的多余变量，然后求解线性方程，对于限定为整数值的变量取与该变量求解结果的最接近整数值作为最终结果，然后将该变量视为常量代入原方程继续求解其他整数变量。

图5为本发明实施例点云数据网格化方法中四边形网格化流程图。点云的参数化结果可以用来建立三角形网格和四边形网格。为了建立三角形网格，对于每一个点及其周围k个近邻点(k＝20)根据参数化值

进行狄洛尼三角化，并取与该点相邻的三角形来建立三维空间中的对应点的连接关系。在分割线附近有重叠的三角形存在，则通过以下步骤来移除多余的三角形并保持三角网格的流形特征：

若有两个以上的三角形共有同一条边，则逐步移除多余的三角形直至每一条边最多被两个三角形共有；

如果移除冗余三角形的过程中产生空洞，则通过对该洞进行狄洛尼三角化来填充该洞。图6为本发明实施例点云数据网格化方法中重叠三角形处理示意图。该步骤输入的是一个原始的三角网格结果如图6a所示，图6b为对图6a中的网格移除重叠三角形后的结果，由于重叠三角形的移除会产生空洞，图6c为对所产生的空洞进行补齐后的三角网格数据。

为了得到最终的四边形网格化结果，本方法首先求取每个三角形内的等值线段，等值线网络构成了基本的四边形网格化基础，等值线的交点即四边形网格的顶点，顶点之间的连接关系则由等值线之间的连接关系决定，对于每一个等值线之间的交点，按照等值线段之间的连接关系找到该交点的邻接点并连接，即可建立形状均匀且方向符合主方向的四边形网格。网格中的每一个四边形形状接近于正方形。

用C++语言实现了本发明所描述的方法，并且在几个不同的数据集上做了实验。所有的实验都是在一台

Core^TM2Quad CPU Q6600，2.40GHz 4G内存的PC机上完成的，显示部分使用了标准的OpenGL图形函数库。

表1给出部分实验所用模型的复杂度(点的个数)和二种处理方法(MIQP和TightCocone+MIQ)所花的时间对比。其中MIQP表示直接在点云数据应用本方法的时间，TightCocone+MIQuelon表示采用传统的首先对点云数据进行网格重建然后对网格使用文献Bommes(BOMMES，D.，ZIMMER，H.，AND KOBBELT，L.2009.Mixed-integer quadrangulation.ACM Transactions on Graphics 28，3，1-10.)中的方法进行四边形网格化所需要的时间，本实验中所使用的网格重建方法为TightCocone(T.K.Dey，S.Goswami，Tight cocone：a water-tight surface reconstructor，in：Proceedings ofthe eighth ACM symposium on Solid modeling and applications，pp.127-134.)

表1

附图7给出了在激光单次扫描数据茶杯模型上得到的结果。图7(a)为原始点云数据，图7a为茶杯模型上的原始点云数据，图7b显示了参数化结果，图7c为利用参数化结果进行四边形网格化的结果。图片显示四边形分布均匀，走向符合物体的几何特征。本结果均在点云模型上处理进行的，为方便显示在网格曲面上显示结果。

附图8列出了在零件的模型上利用我们的方法在零件的点云模型上得到的结果与利用Bommes的方法在零件的网格模型上得到的结果之间的对比。图8a和图8b为Bommes的方法得到的四边形网格化结果和反映四边形网格质量的直方图，图8c和图8d为本方法得到的四边形网格化结果和反映四边形网格质量的直方图，从图中可以看出对于相同的模型，本方法仅利用模型的点云数据而无需事先三角网格化可以得到同样质量的结果。

附图9列出了在更为复杂的机械零件模型的结果。由于噪声和点云分布不均匀的情况存在使得处理更为困难，本方法仍然能够给出良好的结果。其中图9a为零件1模型的原始点云数据，图9b为零件1模型的参数化结果，图9c为根据参数化结果得到的零件1模型的四边形网格化结果，图9d为零件2模型的原始点云数据，图9e为零件2模型的参数化结果，图9f为根据参数化结果得到的零件2模型的四边形网格化结果。

附图10列出了本方法的三角形网格化结果与Floater(M.S.Floater，M.Reimers，Meshless parameterization and surface reconstruction，Comput.Aided Geom.Des.18(2001)77-92.)的方法对比。图10a和图10b分别为Floater的方法和本方法的参数化结果。Floater的方法参数域边界固定为圆形，因此扭曲度较大，本方法无需固定边界，因而能够得到变形更小的参数化结果。注意，变形更小也表现为定义域部分区域的重叠，见图10b，但丝毫不影响参数化的最后结果。图10c和图10d分别为Floater的方法和本发明方法得到的三角网格重建结果对比。可以看出，本方法能够得到更好的三角网格结果，三角网格分布均匀，狭长的三角形较少。而Floater的方法出现了大量的狭长三角形。

本发明直接在点云模型上进行全局参数化和四边形网格化。而且本方法完全自动，无需人工干预，并且通过参数可以快速的控制四边形网格化的疏密程度，得到各种分辨率的四边形网格。本发明可以用于计算机图形学各应用领域，具有较高的实际应用价值。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (14)

1.一种点云数据网格化的方法，其特征在于，该点云数据网格化的方法包括：

步骤A，获取点云数据；

步骤B，由所述点云数据确定每一点的法向量、主方向及点云数据的奇异点信息；

步骤C，由点云数据每一点的坐标，该点的主方向，法向量以及奇异点信息进行点云的参数化；

步骤D，由点云参数化结果，建立与所述点云数据形状一致的三角形网格。

2.根据权利要求1所述的点云数据网格化的方法，其特征在于，所述步骤B包括：

步骤B1；对所述初始点云数据获取每一点的法向量；

步骤B2，由所述法向量计算该点的曲率张量，获取该点初始主方向；

步骤B3，通过使相邻点的初始主方向的一致最大化，对每一点的初始主方向进行平滑，获取每一点的主方向；

步骤B4，根据点云数据中每一点的主方向，确定点云数据中的奇异点。

3.根据权利要求2所述的点云数据网格化的方法，其特征在于，所述步骤B1包括：

步骤B1a，对于点云数据的每一个点p，利用点云数据的kd树查找所述点p的n个近邻点，其中n为整数；

步骤B1b，假设这些近邻点来自于同一个平面T，用这些近邻点到拟合平面的残差的绝对值，再乘以权系数的积的和构造最小二乘问题，其中的权的确定是以点云数据中每一个点与近邻点的欧式距离的倒数作为权值；

步骤B1c，利用最小二乘方法拟合出平面T，以这个平面的法向量作为点云中该点的法向量

4.根据权利要求3所述的点云数据网格化的方法，其特征在于，所述n＝15或30。

5.根据权利要求2所述的点云数据网格化方法，其特征在于，所述步骤B2包括：

步骤B2a，对每一点p建立局部坐标系，设点p的法向量为

则这个p点就是局部坐标系的原点，设局部坐标系的三个方向分别为

步骤B2b，对于p的一个近邻点p_i，其法向量为

则点p处的曲率张量必须满足下述方程式的约束：

$\left[\begin{array}{ccc}{\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{u}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{u}& \stackrel{\→}{N}\·\stackrel{\→}{u}\\ {\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{v}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{v}& \stackrel{\→}{N}\·\stackrel{\→}{v}\\ {\▿}_{\stackrel{\→}{u}}\stackrel{\→}{N}\·\stackrel{\→}{w}& {\▿}_{\stackrel{\→}{v}}\stackrel{\→}{N}\·\stackrel{\→}{w}& \stackrel{\→}{N}\·\stackrel{\→}{w}\end{array}\right]\·\left[\begin{array}{c}{\mathrm{\Δ}}_{\stackrel{\→}{p}}\·\stackrel{\→}{u}\\ {\mathrm{\Δ}}_{\stackrel{\→}{p}}\·\stackrel{\→}{v}\\ 1\end{array}\right]=\left[\begin{array}{c}\stackrel{\→}{N_i}\·\stackrel{\→}{u}\\ \stackrel{\→}{N_i}\·\stackrel{\→}{v}\\ \stackrel{\→}{N_i}\·\stackrel{\→}{w}\end{array}\right];$

步骤B2c，对点p的m个近邻点的法方向进行拟合，代入上述方程中，求解出左边矩阵，即为该点处的曲率张量，矩阵的最大两个特征值所对应的特征向量即为对应的点p的初始主方向

和

6.根据权利要求5所述的点云数据网格化的方法，其特征在于，所述m为15。

7.根据权利要求2所述的点云数据网格化的方法，其特征在于，所述步骤B3包括：

步骤B3a，定义衡量点云数据相邻点的主方向差异的函数

$E_{\mathrm{smoothing}}=\underset{e_{\mathrm{ij}}}{\mathrm{\Σ}}{(({\mathrm{\β}}_i-{\mathrm{\β}}_j)+{\mathrm{\α}}_i-{\mathrm{\α}}_j+\mathrm{n\π}/2)}^2,$

其中，α_i为所要求的主方向与一个参考方向的夹角，该参考方向取为点i的切平面上的任意一个方向，β_i为原始主方向K′与该参考方向的夹角，e_ij为k近邻图中的一条边，n为整数决定两个相邻点的主方向如何变化；

步骤B3b，为了求解上述函数，首先通过近似变化来消除整数变量n，引入变量sin 4α_i和co54α_i来代替原始变量α_i，方程变为：

$E_{\mathrm{smoothing}}=\underset{e_{\mathrm{ij}}}{\mathrm{\Σ}}{\{\left(\begin{array}{c}\cos 4{\mathrm{\α}}_i\\ \sin {4\mathrm{\α}}_i\end{array}\right)-\left(\begin{array}{cc}\cos \mathrm{\β}& -\sin \\ \sin \mathrm{\β}& \cos \mathrm{\β}\end{array}\right)\left(\begin{array}{c}\cos {4\mathrm{\α}}_j\\ \sin {4\mathrm{\α}}_j\end{array}\right)\}}^2;$

其中β＝β_i-β_j，对sin4α_i和co54α_i求导得到线性方程，求解得到sin 4α_i和co54α_i的解，通过sin4α_i和co54α_i最终求出原始变量α_i；

步骤B3c，利用所述夹角α_i，获取点云数据中每一点的主方向。

8.根据权利要求2所述的点云数据网格化的方法，其特征在于，所述步骤B4包括：

步骤B4a，将点p的近邻点和相对应的主方向投影到点p的切平面上；

步骤B4b，在切平面上对这些点按照逆时针顺序排序；

步骤B4c，选择其中一个近邻点在切平面上的主方向作为参考方向，从该参考方向出发，根据上述排序后的顺序求出每两个相邻点的主方向的变化角度；

步骤B4d，若所有变化角度之和位于区间[-π/2，π/2]之外，则该点被标记为奇异点。

9.根据权利要求1所述的点云数据网格化的方法，其特征在于，所述步骤C包括：

步骤C1，对所述点云数据进行切割，转化为亏格为0，边界为一的拓扑结构；

步骤C2，将k近邻图中的每一条边投影到分割线及其法向量所定义的平面上，此处所述分割线为步骤C1中切割所产生的分割线，若投影后的边与分割线相交，则在k近邻图中去除此边，从而得到一个新的k近邻图；

步骤C3，对于定义在点云上每一个点的两个标量函数θ和

求出θ和

在每一个点的具体值，使得该点的最大主方向和最小主方向分别与两个标量函数的梯度尽量保持一致，两个标量函数θ和

为下述方程最小值的最优解：

其中，和

分别为两个函数在每点的梯度，和分别为该点的最大主方向和最小主方向，ω为由用户指定的控制参数化疏密程度的参数。

10.根据权利要求9所述的点云数据网格化的方法，其特征在于，所述步骤C1包括：

步骤C1a，利用点云数据和莫尔斯函数计算出点云的同源基，沿着同源基对点云进行分割；

步骤C2b，在k近邻图上求取每个奇异点到边界的最短路径，并按照连接路径对点云进行再次分割，获取由点云中若干点连接而成的路径。

11.根据权利要求9所述的点云数据网格化的方法，其特征在于，所述步骤C3包括：

步骤C3a，定义衡量标量函数在每一点的梯度与该点的最大主方向和最小主方向之间差异的能量函数：

公式1，

其中，θ_i，

分别表示点i对应的两个标量函数的值，K，K^⊥分别表示点i处的最大主方向和最小主方向，ij为k近邻图中的边，e_ij表示连接点i，j的线段的矢量，w是用户指定的控制参数化疏密程度的参数；

步骤C3b，对于位于分割线上的点P，其近邻点被分割线分成两类，该点也相应的被分成两个点P⁺和P^-，以及对应的主方向K⁺和K^-。为了保证参数值θ和

在分割线处的连续性，点P处的参数值必须满足以下条件：

公式2，

其中，(j，k)为整数变量，

表示方向K^-到K⁺所需要的旋转矩阵；并且，奇异点处的参数值也被限定为整数值；

步骤C3c，通过公式1和公式2建立一个在分割线处连续的全局参数化方程；

步骤C3d，求解该全局参数化方程，获得每一点的标量函数参数值θ和

12.根据权利要求1所述的点云数据网格化的方法，其特征在于，所述步骤D包括：

步骤D1，对于每一个点及其周围k个近邻点根据参数化值

进行狄洛尼三角化；

步骤D2，取与该点相邻的三角形来建立三维空间中的对应点的连接关系；

步骤D3，如果有两个以上的三角形共有同一条边，则逐步移除多余的三角形直至每一条边最多被两个三角形共有；

步骤D4，如果移除冗余三角形的过程中产生空洞，则通过对该洞进行狄洛尼三角化来填充该洞。

13.根据权利要求1所述的点云数据网格化的方法，其特征在于，该方法还包括：

步骤E，在建立的三角网格基础上建立四边形网格。

14.根据权利要求13所述的点云数据网格化的方法，其特征在于，所述步骤E包括：

步骤E1，求取每个三角形内的等值线段；

步骤E2，求取等值线段之间的交点；

步骤E3，对于每一个等值线段之间的交点，按照等值线段之间的连接关系找到该交点的邻接点并连接，建立形状均匀且方向符合主方向的四边形网格。

Images