CN101702237B - 一种基于二次误差测度的层位面拟合方法 - Google Patents

Abstract

一种基于二次误差测度的层位面拟合方法，包括以下步骤：数据输入及预处理；拟合区域确定；计算采样点的二次误差测度矩阵；将采样点分配给每个待调整的网格顶点；利用二次误差测度矩阵调整网格顶点。本发明提出一种基于二次误差测度的层位面拟合方法，可以根据输入的采样点和断层面，对已重构的层位面进行变形，使得变形后的层位面更符合输入点的形状。相对于现有的网格变形方法，本发明使用简便，速度快，可以较好地解决网格模型局部变形的问题。

Description

一种基于二次误差测度的层位面拟合方法

技术领域

本发明属于计算机图形学领域和地质建模领域，具体涉及一种基于二次误差测度的层位面拟合方法。

背景技术

三维网格可用来表示自由曲面，广泛应用于计算机图形学、虚拟现实等领域。随着扫描技术和重建技术的发展，越来越多的模型使用网格来表示。

地质曲面通常分为层位面和断层面，以三角网格来描述。层位面通常被断层面切分成若干个连通或不连通的层位面片。由于断层面的存在，层位面位于断层面附近的左右两部分几何上往往并不连续。如果重建算法是以整张层位面为单位重建，重建出来的三角网格曲面只能反映层位面的全局轮廓。对于更加精细的局部信息，特别是断层面对层位面的切分效果，还需要对层位面片在内边界区域进行拟合处理，也就是对部分与采样点形状不符的网格区域进行变形。

一个好的网格变形方法应当具备操作简单、变形结果符合用户期望的特点。早期主要的网格变形技术是自由变形技术，其算法简单，变形过程几何直观性强。但由于是逐顶点的操作，对于局部细节丰富的模型需要大量的人工调整。这类方法包括形状插值法、自由体变形法、基本骨骼变法、细节保持法等。上世纪90年代出现了多分辨率的网格变形技术，克服了自由变形技术在局部细节编辑上的弱点，支持整体控制和局部细节编辑。但需要构造几何模型的累进表示，或因需半规整重采样而受到限制。

针对上述应用，原有的网格变形操作要么过于复杂、使用不便，要么局部变形效果不好，难以满足该需求，因此需要设计一种操作简单、变形效果好的网格局部变形方法。

发明内容

本发明要解决的问题是：带断层约束的层位面重建时，在断层附近的部分易受断层影响，现有的相应算法要么过于复杂、使用不便，要么局部变形效果不好，难以满足需求。

本发明的技术方案为：一种基于二次误差测度的层位面拟合方法，包括以下步骤：

步骤1：数据输入及预处理

输入的数据包括层位面和断层面以及层位面的初始采样点，层位面和断层面以三角网格表示；

预处理包括层位面与断层面的交线的求取，定义层位面的内边界为层位面与断层面相交处的交线，层位面的外边界为其原来的边界，内边界是由于断层面切割或其它层位面穿刺而造成的，此类断层面和层位面统称为约束曲面，断层面切割层位面有四种类型，分别称为不完全切割、半完全切割-1、半完全切割-2、完全切割，针对以上四种类型应用网格求交算法，得到交线数据；

步骤2：拟合区域确定

拟合时不只对内边界上的网格顶点进行调整，内边界上某网格顶点调整后，对其一定范围内的邻近网格顶点也进行相应调整，调整范围的确定方法为：对于交线上的每个网格顶点V_i，找到其在断层面约束下的最近采样点NearSV，设V_i与NearSV之间距离d＝||V_i-NearSV||，在断层面约束下搜索以V_i为球心d为球半径的范围内的所有网格顶点，这些网格顶点就是需调整的邻近网格顶点，遍历完交线上所有的网格顶点后得到的调整范围内的邻近网格顶点之和，即为此次拟合时内边界附近所有调整顶点；

步骤3：计算采样点的二次误差测度矩阵

在三维欧氏空间中，一个平面可以表示为ax+by+cz+d＝0，平面法向量为n＝(a，b，c)^T，平面方程可写为n^Tv+d＝0，空间一点v＝(x，y，z)^T到该平面距离平方为

D²(v)＝(ax+by+cz+d)²

     ＝(n^Tv+d)²

     ＝(n^Tv+d)(n^Tv+d)

     ＝(v^Tn+d)(n^Tv+d)

     ＝v^T(nn^T)v+2(dn)^Tv+d²

可以记为三元组

Q＝(A，b，c)

Q为二次误差测度矩阵，其中

$A={\mathrm{nn}}^T=\left[\begin{array}{ccc}{a}^2& \mathrm{ab}& \mathrm{ac}\\ \mathrm{ab}& b^2& \mathrm{bc}\\ \mathrm{ac}& \mathrm{bc}& c^2\end{array}\right],$ $b=\mathrm{dn}=\left[\begin{array}{c}\mathrm{ad}\\ \mathrm{bd}\\ \mathrm{cd}\end{array}\right],$ c＝d²

具体针对每个采样点s_i，首先找到其最近k个采样点，k由用户指定，在寻找过程中，保证采样点和其最近点都位于断层面及其延长面的同一边，找到k个最近采样点后，用主成分分析法计算s_i的单位法矢(a，b，c)，d＝(a，b，c)*s_i，则采样点的二次误差测度矩阵可以如下计算：

a²＝a*a；ab＝a*b；ac＝a*c；ad＝a*d；

         b²＝b*b；bc＝b*c；bd＝b*d；

                  c²＝c*c；cd＝c*d；

                           d²＝d*d；

步骤4：将采样点分配给每个待调整的网格顶点

将采样点向以v_i为顶点的所有三角网格做投影，如果某个采样点的投影点落在这些三角网格内，则将该采样点分配给v_i，用于v_i的调整；

步骤5：利用二次误差测度矩阵调整网格顶点

采用二次误差矩阵计算某网格顶点到多个采样点所在切平面的距离平方的和，并利用该式求得网格顶点拟合的最佳位置，调整网格顶点v_i使v_i到所有分配给它的n个采样点切平面的距离的平方和最小：

即 $D^2\left(v_i\right)=\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}(v_i^T{A}_j{v}_i+2b_j^T{v}_i+c_j)=v_i^T\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{A}_j{v}_i+2\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{b}_j^T{v}_i+\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{c}_i=v_i^{T}A{v}_i+2b{v}_i+c$ 最小，

D²(v)为二次表达式，找到其最小值是一个线性问题，对D²(v)求导可得：

$\▿D^2\left(v\right)=2\mathrm{Av}+2b$

解上式，得：

v′_i＝-A^-1b

由采样点二次误差测度矩阵解出v_i的新位置，即新顶点v_i’的位置，对原网格顶点v_i进行调整；对所有要调整的网格顶点使用上述过程进行调整后，则完成了一次拟合。

步骤1中，当层位面未完全穿透断层面，首先在层位面上选定未完全穿透的区域进行延伸，延伸后使选定的区域与断层面形成完全切割类型，最后再应用网格求交算法得到相应交线。

步骤1中，如果交线被两个层位面片共享，则拟合前需根据切割类型在层位面中添加一份交线数据；在层位面添加交线数据就是将交线复制一份，根据切割类型，复制时交线的两个端点的复制方式不一样：

1)、不完全切割的交线的两个端点都不复制，即拟合和层位面在此处还保持连续；

2)、半完全切割的情况下，位于层位面内部的端点不复制，而位于层位面边界的端点复制一份；

3)、完全切割的情况下，交线的两个端点都复制一份。

进一步的，本发明进行多次拟合操作提高拟合效果。

本发明提出了点云的二次误差测度计算方法；将采样点分配给每个要变形的网格顶点，这是本发明的关键之一，并且，本发明网格顶点的变形不受法矢限制，更加光顺。本发明在调整网格顶点位置时，未限定只能延网格顶点的法矢方向调整，而是通过最小化网格面得到采样点的二次误差，进而调整网格顶点位置，这样的好处是可以获得更加光顺的变形结果。

本发明提出一种基于二次误差测度的层位面拟合方法，可以根据输入的采样点和断层面，对已重构的层位面进行变形，使得变形后的层位面更符合输入点的形状。相对于现有的网格变形方法，本发明使用简便，速度快，可以较好地解决网格模型局部变形的问题。

附图说明

图1为本发明方法的流程图。

图2为断层面切割层位面的四种类型，(a)为不完全切割、(b)为半完全切割-1、(c)为半完全切割-2、(d)为完全切割。

图3为层位面未完全穿透断层面，(a)为进行曲面延伸前，(b)为曲面延伸后。

图4为本发明确定拟合区域的示意图。

图5为本发明拟合效果图。

具体实施方式

本发明方法如图1所示，包括数据输入及预处理；拟合区域确定；计算采样点的二次误差测度矩阵；将采样点分配给每个待调整的网格顶点五个步骤，利用二次误差测度矩阵调整网格顶点。

下面结合附图进一步对本发明加以说明：

步骤1：数据输入及预处理

本发明用来处理重建后的层位面。因此输入的数据包括层位面和断层面以及层位面的初始采样点，层位面和断层面以三角网格描述。层位面通常被断层切分成若干相连或不相连的区域。

预处理主要包括交线的求取。定义层位面的内边界为其与断层面相交处的交线，层位面的外边界为其原来的边界。内边界主要是由于断层面切割或其它层位面穿刺等情况造成的，此类断层面和层位面可统称为约束曲面。断层面切割层位面会有四种类型，如图2(a)、(b)、(c)、(d)所示，分别称为不完全切割、半完全切割-1、半完全切割-2、完全切割。针对以上四种类型应用网格求交算法，得到交线数据，如果交线被两个层位面片共享，则拟合前需根据切割类型在层位面中添加一份交线数据。在层位面添加交线数据就是将交线复制一份，根据切割类型，复制时交线的两个端点的复制方式不一样：

1)、不完全切割的交线的两个端点都不复制，即拟合和层位面在此处还保持连续；

2)、半完全切割的情况下，位于层位面内部的端点不复制，而位于层位面边界的端点复制一份；

3)、完全切割的情况下，交线的两个端点都复制一份。

当出现如图3所示的切割情况时，层位面未完全穿透断层面，首先在层位面上选定区域进行延伸，选定未完全穿透部分，延伸后使选定区域与断层面形成完全切割类型，最后再应用网格求交算法得到相应交线，如图3(b)。

步骤2：拟合区域确定

拟合时不可能只对内边界上的顶点进行调整，而是调整时有个波动影响范围，即内边界上某点调整后必须对其一定范围内的邻近顶点也进行相应调整，调整范围的确定方法如图4所示，对于交线上的每个网格顶点V_i，找到在断层面约束下的最近采样点NearSV，设V_i与NearSV之间距离d＝||V_i-NearSV||，在断层面约束下搜索以V_i为球心d为球半径的范围内的所有网格顶点，这些网格顶点就是需调整的顶点，遍历完交线上所有的网格顶点后得到的调整范围之和即为此次拟合时内边界附近所有调整顶点。图4a)中找到的最近采样点NearSV距离V_i较近，所以球范围内顶点只有4个，也就是需调整的顶点有4个，而在图4(b)中点NearSV与V_i距离较远，调整范围就相应扩大为6个。

步骤3：计算采样点的二次测度误差矩阵

在三维欧氏空间中，一个平面可以表示为ax+by+cz+d＝0，平面法向量为n＝(a，b，c)^T，平面方程可写为n^Tv+d＝0，空间一点点v＝(x，y，z)^T到该平面距离平方为

D²(v)＝(ax+by+cz+d)²

     ＝(n^Tv+d)²

     ＝(n^Tv+d)(n^Tv+d)

     ＝(v^Tn+d)(n^Tv+d)

     ＝v^T(nn^T)v+2(dn)^Tv+d²

可以记为三元组

Q＝(A，b，c)

Q为二次误差测度矩阵，其中

$A={\mathrm{nn}}^T=\left[\begin{array}{ccc}{a}^2& \mathrm{ab}& \mathrm{ac}\\ \mathrm{ab}& b^2& \mathrm{bc}\\ \mathrm{ac}& \mathrm{bc}& c^2\end{array}\right],$ $b=\mathrm{dn}=\left[\begin{array}{c}\mathrm{ad}\\ \mathrm{bd}\\ \mathrm{cd}\end{array}\right],$ c＝d²

具体针对每个采样点s_i，首先找到其最近k个采样点，k可以由用户指定，通过最小化网格面得到采样点的二次误差。在寻找过程中，应保证采样点和其最近采样点都位于断层面及其延长面的同一边。找到k个最近采样点后，用主成分分析法计算s_i的单位法矢(a，b，c)。d＝(a，b，c)*s_i，则采样点的二次测度误差矩阵可以如下计算：

a²＝a*a；ab＝a*b；ac＝a*c；ad＝a*d；

         b²＝b*b；bc＝b*c；bd＝b*d；

                  c²＝c*c；cd＝c*d；

                           d²＝d*d；

步骤4：将采样点分配给每个待调整的网格顶点

并不是所有的采样点对每个网格顶点都起影响作用，必须确定调整某网格顶点时要使用哪些采样点。调整网格顶点v_i时，其对网格形状的影响仅局限在以其为顶点的所有三角网格内。因此，对调整给顶点有所以得采样点也应该位于该范围内。具体而言，本发明将采样点向以v_i为顶点的所有三角网格做投影。如果某个采样点的投影点落在这些三角网格内，则将该采样点分配给v_i时，也即意味着调整v_i时，要用到这些采样点。

步骤5：利用二次误差测度矩阵调整网格顶点

调整网格顶点v_i时，希望v_i到所有分配给它的n个采样点切平面的距离的平方和最小。

即 $D^2\left(v_i\right)=\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}(v_i^T{A}_j{v}_i+2b_j^T{v}_i+c_j)=v_i^T\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{A}_j{v}_i+2\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{b}_j^T{v}_i+\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{c}_i=v_i^{T}A{v}_i+2b{v}_i+c$ 最小。

显然，D²(v)是二次表达式，找到其最小值是一个线性问题。对D²(v)求导可得：

$\▿D^2\left(v\right)=2\mathrm{Av}+2b$

解上式，得：

v′_i＝-A^-1b

解出v_i的新位置，即新顶点v_i’的位置，就可以将其更新，完成网格顶点的调整。

对所有要调整的网格顶点使用上述过程进行调整后，则完成了一次拟合。要取得好的拟合效果可能需要进行多次拟合操作。

图5(a)为拟合前的层位面，图5(b)为本发明方法拟合后的层位面，可见本发明对于更加精细的局部信息，特别是断层面对层位面的切分有十分好的拟合效果。

Claims (4)

1.一种基于二次误差测度的层位面拟合方法，其特征是包括以下步骤：

步骤1：数据输入及预处理

输入的数据包括层位面和断层面以及层位面的初始采样点，层位面和断层面以三角网格表示；

预处理包括层位面与断层面的交线的求取，定义层位面的内边界为层位面与断层面相交处的交线，层位面的外边界为其原来的边界，内边界是由于断层面切割或其它层位面穿刺而造成的，此类断层面和层位面统称为约束曲面，断层面切割层位面有四种类型，分别称为不完全切割、半完全切割-1、半完全切割-2、完全切割，针对以上四种类型应用网格求交算法，得到交线数据；

步骤2：拟合区域确定

拟合时不只对内边界上的网格顶点进行调整，内边界上某网格顶点调整后，对其一定范围内的邻近网格顶点也进行相应调整，调整范围的确定方法为：对于交线上的每个网格顶点V_i，找到其在断层面约束下的最近采样点NearSV，设V_i与NearSV之间距离dist＝||V_i-NearSV||，在断层面约束下搜索以V_i为球心dist为球半径的范围内的所有网格顶点，这些网格顶点就是需调整的邻近网格顶点，遍历完交线上所有的网格顶点后得到的调整范围内的邻近网格顶点之和，即为此次拟合时内边界附近所有调整顶点；

步骤3：计算采样点的二次误差测度矩阵

在三维欧氏空间中，一个平面表示为n^Tv+d＝0，其中n为平面单位法向量，表示为n＝(a，b，c)^T，d为一标量，空间一点v＝(x，y，z)^T到该平面距离平方为

D²(v)＝(ax+by+cz+d)²

＝(n^Tv+d)²

＝(n^Tv+d)(n^Tv+d)

＝(v^Tn+d)(n^Tv+d)

＝v^T(nn^T)v+2(dn)^Tv+d²

记为三元组

Q＝(A，B，C)

Q为二次误差测度矩阵，其中

$A={\mathrm{nn}}^T=\left[\begin{array}{ccc}{a}^2& \mathrm{ab}& \mathrm{ac}\\ \mathrm{ab}& b^2& \mathrm{bc}\\ \mathrm{ac}& \mathrm{bc}& c^2\end{array}\right],$ $B=\mathrm{dn}=\left[\begin{array}{c}\mathrm{ad}\\ \mathrm{bd}\\ \mathrm{cd}\end{array}\right],$ C＝d²

具体针对每个采样点s_i，首先找到其最近k个采样点，k由用户指定，在寻找过程中，保证采样点和其最近点都位于断层面及其延长面的同一边，找到k个最近采样点后，用主成分分析法计算s_i的单位法矢n_i，记为(a_i，b_i，c_i)^T，则d_i＝(a_i，b_i，c_i)^T·s_i，则采样点的二次误差测度矩阵如下计算：

a_i ²＝a_i*a_i；a_ib_i＝a_i*b_i；a_ic_i＝a_i*c_i；a_id_i＝a_i*d_i；

b_i ²＝b_i*b_i；b_ic_i＝b_i*c_i；b_id_i＝b_i*d_i；

c_i ²＝c_i*c_i；c_id_i＝c_i*d_i；

$d_i^2=d_i*d_i;$

步骤4：将采样点分配给每个待调整的网格顶点

将采样点向以v_i为顶点的所有三角网格做投影，如果某个采样点的投影点落在这些三角网格内，则将该采样点分配给v_i，用于v_i的调整；

步骤5：利用二次误差测度矩阵调整网格顶点

采用二次误差矩阵计算某网格顶点到多个采样点所在切平面的距离平方的和，并利用该式求得网格顶点拟合的最佳位置，调整网格顶点v_i使v_i到所有分配给它的num_i个采样点切平面的距离的平方和最小：

即

$D^2\left(v_i\right)=\underset{j=0}{\overset{{\mathrm{num}}_i-1}{\mathrm{\Σ}}}({v_i}^T{A}_j{v}_i+2{B_j}^T{v}_i+C_j)={v_i}^T\underset{j=0}{\overset{{\mathrm{num}}_i-1}{\mathrm{\Σ}}}{A}_j{v}_i+2\underset{j=0}{\overset{{\mathrm{num}}_i-1}{\mathrm{\Σ}}}{B_j}^T{v}_i+\underset{j=0}{\overset{{\mathrm{num}}_i-1}{\mathrm{\Σ}}}{C}_j={v_i}^T{A}_i{v}_i+2B_i{v}_i+C_i$ 最小，

D²(v)为二次表达式，找到其最小值是一个线性问题，对D²(v)求导可得：

$\▿D^2\left(v\right)=2A_{i}v+2B_i$

解上式，得：

$v_i^{\′}=-A_i^{-1}{B}_i$

由采样点二次误差测度矩阵解出v_i的新位置，即新顶点v_i’的位置，对原网格顶点v_i进行调整；对所有要调整的网格顶点使用上述过程进行调整后，则完成了一次拟合。

2.根据权利要求1所述的一种基于二次误差测度的层位面拟合方法，其特征是步骤1中，当层位面未完全穿透断层面，首先在层位面上选定未完全穿透的区域进行延伸，延伸后使选定的区域与断层面形成完全切割类型，最后再应用网格求交算法得到相应交线。

3.根据权利要求1所述的一种基于二次误差测度的层位面拟合方法，其特征是步骤1中，如果交线被两个层位面片共享，则拟合前需根据切割类型在层位面中添加一份交线数据；在层位面添加交线数据就是将交线复制一份，根据切割类型，复制时交线的两个端点的复制方式不一样：

1)、不完全切割的交线的两个端点都不复制，即拟合和层位面在此处还保持连续；

2)、半完全切割的情况下，位于层位面内部的端点不复制，而位于层位面边界的端点复制一份；

3)、完全切割的情况下，交线的两个端点都复制一份。

4.根据权利要求1-3任一项所述的一种基于二次误差测度的层位面拟合方法，其特征是进行多次拟合操作提高拟合效果。

Images