CN101833790B - 一种基于波动方程的各向异性的四边形网格生成方法 - Google Patents

Abstract

本发明公开了一种基于波动方程的各向异性的四边形网格生成方法。包括如下步骤：1)根据用户需求和模型特点，生成相应的特征约束、方向场和各向异性的密度场；2)根据特征约束、方向场和各向异性的密度场，在网格表面上构造二维驻波；3)从二维驻波中形成模型的四边形剖分；4)在四边形剖分的基础上构造各向异性的全局参数化，然后得到最终的四边形化网格。本发明可以生成各向异性的四边形网格，同时能够对四边形网格的形状，密度，方向，特征对齐，奇异点分布等需求进行灵活的控制和优化。

Description

一种基于波动方程的各向异性的四边形网格生成方法

技术领域

本发明涉及四边形网格生成方法，尤其涉及一种基于波动方程的各向异性的四边形网格生成方法。

背景技术

四边形网格在许多领域有着广泛的应用，比如纹理贴图，有限元模拟，B样条表面重建等。四边形网格生成的主要难点在于网格质量严重依赖于生成时对模型特征和用户多样需求的全局考虑。一般来说，我们从以下方面评价四边形网格的质量：

(1)面片质量：每个四边形面片需要尽量的是正方形或者长方形。

(2)方向：四边形网格的边需要对应于网格的主曲率方向或者是用户指定的方向。

(3)特征对齐：模型上的特征线需要被连续的四边形网格的边表示，以尽量的减少原模型和四边形化后的模型之间的Hausdorff距离。

(4)全局结构：四边形网格中奇异点的位置和数量需要对应模型全局的几何结构和用户需求。

目前四边形网格生成算法主要分以下几种：

1.在模型上沿着给定的方向场跟踪，生成曲线，这些线相交就得到了四边形网格。(可参考ALLIEZ，P.，COHEN-STEINER，D.，DEVILLERS，O.，L′EVY，B.，AND DESBRUN，M.2003.Anisotropic polygonal remeshing.ACM Trans.Graph.22，3，485-493)

2.全局参数化的方法，通过在网格表面构造一个参数化域，使得参数化域的坐标轴与用户指定的方向对齐。(可参考RAY，N.，LI，W.C.，L′EVY，B.，SHEFFER，A.，AND ALLIEZ，P.2006.Periodic global parameterization.ACM Trans.Graph.25，4，1460-1485.)

3.网格分片的方法，基于用户的要求(方向，特征，大小等)，将模型分成一些大的片，然后把这些片转化成为四边形。(可参考DONG，S.，BREMER，P.-T.，GARLAND，M.，PASCUCCI，V.，AND HART，J.C.2006.Spectral surfacequadrangulation.ACM Trans.Graph.25，3，1057-1066)

发明内容

本发明的目的在于针对现有技术的不足，提出了一种基于波动方程的各向异性的四边形网格生成方法。

为实现上述目的，本发明采用的技术方案如下：

基于波动方程的各向异性的四边形网格生成方法包括如下步骤：

1)根据用户需求和模型特点，生成相应的特征约束、方向场和各向异性的密度场；

2)根据特征约束、方向场和各向异性的密度场，在网格表面上构造二维驻波；

3)从二维驻波中形成模型的四边形剖分；

4)在四边形剖分的基础上构造各向异性的全局参数化，然后得到最终的四边形化网格。

所述的根据用户需求和模型特点，生成相应的特征约束、方向场和各向异性的密度场的步骤为：

1)网格上二面角大于45度的边会被认为是特征边，相连的特征边会被组合成特征线约束，同时用户直接在模型上指定特征约束；

2)使用网格上的主曲率方向作为初始值，用户在模型上指定的方向和特征约束的方向为边界条件，光顺去噪后得到方向场；

3)根据方向场的旋度约束能量和用户提出的边界约束能量、各向异性约束能量和光顺约束能量，优化得到各向异性的密度场。

所述的根据特征约束、方向场和各向异性的密度场，在网格表面上构造二维驻波的步骤为：

对于网格上任意边e_ij＝<p_i，p_j>，顶点p_i，p_j之间的二维驻波相位角关系定义如下： ${\mathrm{\&alpha;}}_{\mathrm{ij}}={\mathrm{\&theta;}}_j-{\mathrm{\&theta;}}_i\&ap;{\mathrm{\&omega;}}_x\stackrel{\&RightArrow;}{X}\left[e_{\mathrm{ij}}\right]\&CenterDot;(p_j-p_i)$ ${\mathrm{\&beta;}}_{\mathrm{ij}}={\mathrm{\&phi;}}_j-{\mathrm{\&phi;}}_i\&ap;{\mathrm{\&omega;}}_y\stackrel{\&RightArrow;}{Y}\left[e_{\mathrm{ij}}\right]\&CenterDot;(p_j-p_i)$

其中(θ_i，φ_i)和(θ_j，φ_j)是顶点p_i，p_j上的二维驻波相位角，ω_x和ω_y是在方向

上的波速；顶点p_i，p_j之间的二维驻波振幅关系定义如下：

f(p_j)＝cos(θ_j)cos(φ_j)≈cos(θ_i+α_ij)cos(φ_i+β_ij)＝c_ij·f(p_i)

其中f(p_i)为顶点p_i上的二维驻波振幅，c_ij为顶点p_i，p_j之间的二维驻波振幅张量，根据振幅关系在网格表面上构造二维驻波，二维驻波的极值点连续分布在特征线上，相邻极值点之间的连线方向对应于方向场的方向，周期正比于密度场。

所述的从二维驻波中形成模型的四边形剖分的步骤为：

1)找出二维驻波场的极大点、极小点和鞍点；

2)从每个二维驻波场的鞍点出发，沿着二维驻波场梯度上升/下降最快的方向在模型表面上前进，直到遇到二维驻波场的极大/极小点，得到剖分图，采用剖分图的对偶图形成模型的四边形剖分。

所述的在四边形剖分的基础上构造各向异性的全局参数化，然后得到最终的四边形化网格的步骤为：

1)为剖分得到的每个四边形区域确定参数化映射空间；

2)根据参数化映射空间，为任意相邻的四边形区域确定参数化坐标之间的转换函数；

3)利用相邻四边形区域参数化坐标之间的转换函数，构造各向异性的全局参数化；

4)根据全局参数化的结果生成最终的四边形化网格。

本发明与背景技术相比较，其优点在于：

1)可以生成各向异性的四边形网格，同时网格中的奇异点分布和数量将被自动优化。

2)可以对四边形网格的形状，密度，方向，特征对齐的需求进行灵活的控制。用户需求将被纳入一个波动方程的求解框架中统一优化。

3)可以很好的处理复杂的模型(高亏格或者多样的特征约束)和极端情况(很小的环或者柄)。

附图说明

下面结合附图和实施例对本发明作进一步说明：

图1是本发明方法的流程图；

图2是生成的特征约束、方向场和各向异性的密度场；

图3是在上面的特征约束、方向和密度下求解得到的二维驻波；

图4是形成的模型四边形剖分；

图5是全局参数化后的四边形网格。

具体实施方式

基于波动方程的各向异性的四边形网格生成方法包括如下步骤：

1)根据用户需求和模型特点，生成相应的特征约束、方向场和各向异性的密度场；

2)根据特征约束、方向场和各向异性的密度场，在网格表面上构造二维驻波；

3)从二维驻波中形成模型的四边形剖分；

4)在四边形剖分的基础上构造各向异性的全局参数化，然后得到最终的四边形化网格。

所述的根据用户需求和模型特点，生成相应的特征约束、方向场和各向异性的密度场的步骤为：

1)网格上二面角大于45度的边会被认为是特征边，相连的特征边会被组合成特征线约束，同时用户直接在模型上指定特征约束；

2)使用网格上的主曲率方向作为初始值，用户在模型上指定的方向和特征约束的方向为边界条件，光顺去噪后得到方向场；

3)根据方向场的旋度约束能量和用户提出的边界约束能量、各向异性约束能量和光顺约束能量，优化得到各向异性的密度场。

各向异性的密度场是定义在网格边上的一组变量(μ，η)。网格上任一面片T上方向场的旋度算子定义如下：

$\mathrm{curl}\left(\stackrel{\&RightArrow;}{X}\right)\&ap;\frac{1}{\left|T\right|}{\&Integral;}_{\&PartialD;T}\stackrel{\&RightArrow;}{X}\&CenterDot;\mathrm{ds}\&ap;\frac{1}{\left|T\right|}\underset{e=e_i,e_j,e_k}{\mathrm{\&Sigma;}}\stackrel{\&RightArrow;}{X}\left[e\right]\&CenterDot;\stackrel{\&RightArrow;}{e}$

其中是|T|面片T的面积，e_i，e_j，e_k是面片T是3条边，

是边e上的

方向，

是边的方向。密度场μ需要满足的旋度条件是：

${\stackrel{\&RightArrow;}{X}}^{\&perp;}\&CenterDot;\&dtri;\widetilde{\mathrm{\&mu;}}=\mathrm{curl}\left(\stackrel{\&RightArrow;}{X}\right)$

其中

是

逆时针旋转90度后的方向，

是

在面片T内的梯度，定义如下：

$\&dtri;\widetilde{\mathrm{\&mu;}}=\underset{e=e_i,e_j,e_k}{\mathrm{\&Sigma;}}{\widetilde{\mathrm{\&mu;}}}_e\&dtri;B_{T,e}$

其中是边e上的B_T，e＝(B_T，i(e)，B_T，j(e)，B_T，k(e))是边e的中点在面片T顶点(i，j，k)上对应的重心坐标，

是重心坐标在面片T内的梯度。

密度场η的旋度条件可以用上面同样的步骤得到。

所以各向异性的密度场(μ，η)的旋度能量定义如下：

$E_{\mathrm{curl}(\widetilde{\mathrm{\&mu;}},\widetilde{\mathrm{\&eta;}})}=\underset{T\&Element;K}{\mathrm{\&Sigma;}}{(\underset{e=e_i,e_j,e_k}{\mathrm{\&Sigma;}}{\widetilde{\mathrm{\&mu;}}}_e{\&dtri;B}_{T,e}-\frac{1}{\left|T\right|}\underset{e=e_i,e_j,e_k}{\mathrm{\&Sigma;}}\stackrel{\&RightArrow;}{X}[e]\&CenterDot;\stackrel{\&RightArrow;}{e})}^2+$

$\underset{T\&Element;K}{\mathrm{\&Sigma;}}{(\underset{e=e_i,e_j,e_k}{\mathrm{\&Sigma;}}{\widetilde{\mathrm{\&eta;}}}_e\&dtri;B_{T,e}-\frac{1}{\left|T\right|}\underset{e=e_i,e_j,e_k}{\mathrm{\&Sigma;}}\stackrel{\&RightArrow;}{Y}[e]\&CenterDot;\stackrel{\&RightArrow;}{e})}^2$

我们同时可以对密度场(μ，η)进行多样的约束。这里我们提供3种：边界约束，各向异性约束，光顺约束。

边界约束指的是用户可以指定特定边e_i上的密度μ_i为s_i，e_j边的密度η_j为t_j。我们把这样的约束定义为如下能量：

$E_c(\widetilde{\mathrm{\&mu;}},\widetilde{\mathrm{\&eta;}})=\underset{i}{\mathrm{\&Sigma;}}{({\widetilde{\mathrm{\&mu;}}}_i-\log \left(s_i\right))}^2+\underset{j}{\mathrm{\&Sigma;}}{({\widetilde{\mathrm{\&eta;}}}_j-\log \left(t_j\right))}^2$

各向异性约束指的是用户可以指定特定边e上密度μ，η的比例为μ_e/η_e＝γ_e。我们把这样的约束定义为如下能量：

$E_a(\widetilde{\mathrm{\&mu;}},\widetilde{\mathrm{\&eta;}})=\underset{e}{\mathrm{\&Sigma;}}{({\widetilde{\mathrm{\&mu;}}}_e-{\widetilde{\mathrm{\&eta;}}}_e-\log \left({\mathrm{\&gamma;}}_e\right))}^2$

光顺约束指的是用户可以控制密度场的光滑程度，我们使用如下能量来衡量：

$E_s(\widetilde{\mathrm{\&mu;}},\widetilde{\mathrm{\&eta;}})={\left|\right|\mathrm{LP}\widetilde{\mathrm{\&mu;}}\left|\right|}^2+{\left|\right|\mathrm{LP}\widetilde{\mathrm{\&eta;}}\left|\right|}^2$

其中L是Laplacian矩阵，P是顶点相对边中点的插值矩阵。

把上面的能量结合起来，我们通过优化如下总能量，得到所需的各向异性的密度场(μ，η)：

$E_{\mathrm{total}}(\widetilde{\mathrm{\&mu;}},\widetilde{\mathrm{\&eta;}})=E_{\mathrm{curl}}(\widetilde{\mathrm{\&mu;}},\widetilde{\mathrm{\&eta;}})+{\mathrm{\&omega;}}_c{E}_c(\widetilde{\mathrm{\&mu;}},\widetilde{\mathrm{\&eta;}})+{\mathrm{\&omega;}}_a{E}_a(\widetilde{\mathrm{\&mu;}},\widetilde{\mathrm{\&eta;}})+{\mathrm{\&omega;}}_s{E}_s(\widetilde{\mathrm{\&mu;}},\widetilde{\mathrm{\&eta;}})$

其中ω_c，ω_a，ω_s是加权系数。

所述的根据特征约束、方向场和各向异性的密度场，在网格表面上构造二维驻波的步骤为：

对于网格上任意边e_ij＝<p_i，p_j>，顶点p_i，p_j之间的二维驻波相位角关系定义如下： ${\mathrm{\&alpha;}}_{\mathrm{ij}}={\mathrm{\&theta;}}_j-{\mathrm{\&theta;}}_i\&ap;{\mathrm{\&omega;}}_x\stackrel{\&RightArrow;}{X}\left[e_{\mathrm{ij}}\right]\&CenterDot;(p_j-p_i)$ ${\mathrm{\&beta;}}_{\mathrm{ij}}={\mathrm{\&phi;}}_j-{\mathrm{\&phi;}}_i\&ap;{\mathrm{\&omega;}}_y\stackrel{\&RightArrow;}{Y}\left[e_{\mathrm{ij}}\right]\&CenterDot;(p_j-p_i)$

其中(θ_i，φ_i)和(θ_j，φ_j)是顶点p_i，p_j上的二维驻波相位角，ω_x和ω_y是在方向

上的波速；顶点p_i，p_j之间的二维驻波振幅关系定义如下：

f(p_j)＝cos(θ_j)cos(φ_j)≈cos(θ_i+α_ij)cos(φ_i+β_ij)＝c_ij·f(p_i)

其中f(p_i)为顶点p_i上的二维驻波振幅，c_ij为顶点p_i，p_j之间的二维驻波振幅张量，根据振幅关系在网格表面上构造二维驻波，二维驻波的极值点连续分布在特征线上，相邻极值点之间的连线方向对应于方向场的方向，周期正比于密度场。得到方向场各向异性的密度场(μ，η)和特征约束之后，我们构造网格上的二维驻波。对于任一条边e_ij＝<p_i，p_j>，相邻顶点p_i，p_j上的驻波相位角为(θ_i，φ_i)和(θ_j，φ_j)，边e_ij沿方向

的波速为ω_x＝π/μ和ω_y＝π/η。因此p_i，p_j的相位角之间存在如下关系：

${\mathrm{\&alpha;}}_{\mathrm{ij}}={\mathrm{\&theta;}}_j-{\mathrm{\&theta;}}_i\&ap;{\mathrm{\&omega;}}_x\stackrel{\&RightArrow;}{X}\left[e_{\mathrm{ij}}\right]\&CenterDot;(p_j-p_i)$ ${\mathrm{\&beta;}}_{\mathrm{ij}}={\mathrm{\&phi;}}_j-{\mathrm{\&phi;}}_i\&ap;{\mathrm{\&omega;}}_y\stackrel{\&RightArrow;}{Y}\left[e_{\mathrm{ij}}\right]\&CenterDot;(p_j-p_i)$

假设顶点p_i上的驻波振幅为f(p_i)，相邻顶点p_i，p_j上的振幅存在如下关系：

f(p_j)＝cos(θ_j)cos(φ_j)≈cos(θ_i+α_ij)cos(φ_i+β_ij)＝c_ij·f(p_i)

其中 $c_{\mathrm{ij}}=\left(\begin{array}{c}\cos \left({\mathrm{\&alpha;}}_{\mathrm{ij}}\right)\cos \left({\mathrm{\&beta;}}_{\mathrm{ij}}\right)\\ -\cos \left({\mathrm{\&alpha;}}_{\mathrm{ij}}\right)\sin \left({\mathrm{\&beta;}}_{\mathrm{ij}}\right)\\ -\sin \left({\mathrm{\&alpha;}}_{\mathrm{ij}}\right)\cos \left({\mathrm{\&beta;}}_{\mathrm{ij}}\right)\\ \sin \left({\mathrm{\&alpha;}}_{\mathrm{ij}}\right)\sin \left({\mathrm{\&beta;}}_{\mathrm{ij}}\right)\end{array}\right),$ $f\left(p_i\right)=\left(\begin{array}{c}\cos \left({\mathrm{\&theta;}}_i\right)\cos \left({\mathrm{\&phi;}}_i\right)\\ \cos \left({\mathrm{\&theta;}}_i\right)\sin \left({\mathrm{\&phi;}}_i\right)\\ \sin \left({\mathrm{\&theta;}}_i\right)\cos \left({\mathrm{\&phi;}}_i\right)\\ \sin \left({\mathrm{\&theta;}}_i\right)\sin \left({\mathrm{\&phi;}}_i\right)\end{array}\right).$

为了保证求得的是二维驻波，我们仍需要2个额外约束：

f(p_i)·f(p_i)＝||cos(θ_i)cos(φ_i)||²+||cos(θ_i)sin(φ_i)||²+||sin(θ_i)cos(φ_i)||²+||sin(θ_i)sin(φ_i)||²＝1

f^t(p_i)Hf(p_i)＝cos(θ_i)cos(φ_i)·sin(θ_i)sin(φ_i)-cos(θ_i)sin(φ_i)·sin(θ_i)cos(φ_i)＝0

其中

是一个4×4矩阵，g_i是一个4维的只有i行是1其他行全0向量。

考虑网格上的所有边，我们通过优化如下能量求得驻波：

$E_{\mathrm{wave}}(f\left(p_1\right),...,f\left(p_n\right))=\underset{e_{\mathrm{ij}}\&Element;K}{\mathrm{\&Sigma;}}\left(\begin{array}{c}{(e_1\&CenterDot;f\left(p_j\right)-c_{\mathrm{ij}}\&CenterDot;f\left(p_i\right))}^2\\ +{(e_1\&CenterDot;f\left(p_i\right)-c_{\mathrm{ji}}\&CenterDot;f\left(p_j\right))}^2\end{array}\right)$

$+{\mathrm{\&omega;}}_1(\underset{p_i\&Element;P}{\mathrm{\&Sigma;}}{(f\left(p_i\right)\&CenterDot;f\left(p_i\right)-1)}^2+\underset{p_i\&Element;P}{\mathrm{\&Sigma;}}{\left(f^t\left(p_i\right)\mathrm{Hf}\left(p_i\right)\right)}^2)$

特征控制需要使得二维驻波抽出的Morse-Small Complex恰好经过特征线，特征线上的点的相位角(θ，φ)需要满足以下条件：θ＝kπφ＝kπ

所以我们用sin(θ)sin(φ)＝0来进行特征控制。

所述的从二维驻波中形成模型的四边形剖分的步骤为：

1)找出二维驻波场的极大点、极小点和鞍点；

2)从每个二维驻波场的鞍点出发，沿着二维驻波场梯度上升/下降最快的方向在模型表面上前进，直到遇到二维驻波场的极大/极小点，得到剖分图，采用剖分图的对偶图形成模型的四边形剖分。

所述的在四边形剖分的基础上构造各向异性的全局参数化，然后得到最终的四边形化网格的步骤为：

1)为剖分得到的每个四边形区域确定参数化映射空间；

2)根据参数化映射空间，为任意相邻的四边形区域确定参数化坐标之间的转换函数；

3)利用相邻四边形区域参数化坐标之间的转换函数，构造各向异性的全局参数化；

4)根据全局参数化的结果生成最终的四边形化网格。我们使用了各向异性的密度场来指导驻波的生成，所以这里需要构造一个各向异性的全局参数化。对于剖分得到的每个四边形区域U_i，我们首先使用该区域所有顶点上的密度场的平均值(w_i，h_i)作为U_i的参数化映射空间，对于任意相邻的两个四边形区域U_α和U_β，我们定义转换函数φ_βα，用来把的U_β中的参数化坐标转换到U_α中。该函数由U_α和U_β的参数化映射空间决定。获得所有相邻区域之间的转换函数之后，使用这些转换函数替换DONG，S.，BREMER，P.-T.，GARLAND，M.，PASCUCCI，V.，AND HART，J.C.2006.Spectral surface quadrangulation.ACM Trans.Graph.25，3，1057-1066.里全局参数化方法中的转换函数，构造我们各向异性的全局参数化。为了保持特征约束，我们需要在参数化中额外为特征线添加硬约束：对于特征线上的每个顶点，我们要求它的参数化坐标始终保持在其所属区域的参数化边界上。得到全局参数化后，我们根据用户的要求在每个四边形区域的参数化空间里均匀细分，得到最终的四边形化结果。

实施例：

本发明提出一种基于波动方程的各向异性的四边形网格生成方法，流程图如图1所示，下面介绍处理一个实际模型的例子：

1)根据网格的二面角信息和用户需求得到相应的特征约束，用特征约束和用户需求为边界条件，光顺主曲率方向得到方向场，通过优化密度场能量方程，从方向场中生成各向异性的密度场，见图2；

2)根据特征约束、方向场和各向异性的密度场，通过优化驻波能量方程，得到网格表面上的二维驻波，见图3；

3)找到二维驻波中的极大点，极小点和鞍点，由鞍点出发，沿着梯度最大/最小的方向前进，得到对模型的剖分图，采用剖分图的对偶图形成模型的四边形剖分，见图4；

4)在得到四边形剖分后，为每个四边形区域确定参数化映射空间，然后为任意相邻的四边形区域确定，构造各向异性的全局参数化，根据用户的要求在每个四边形区域的参数化空间里均匀细分，得到最终的四边形化网格，见图5。

Claims (4)

1.一种基于波动方程的各向异性的四边形网格生成方法，其特征在于包括如下步骤：

1)根据用户需求和模型特点，生成相应的特征约束、方向场和各向异性的密度场；

2)根据特征约束、方向场和各向异性的密度场，在网格表面上构造二维驻波；

3)从二维驻波中形成模型的四边形剖分；

4)在四边形剖分的基础上构造各向异性的全局参数化，然后得到最终的四边形化网格。

2.根据权利要求1所述的一种基于波动方程的各向异性的四边形网格生成方法，其特征在于所述的根据用户需求和模型特点，生成相应的特征约束、方向场和各向异性的密度场的步骤为：

1)网格上二面角大于45度的边会被认为是特征边，相连的特征边会被组合成特征线约束，同时用户直接在模型上指定特征约束；

2)使用网格上的主曲率方向作为初始值，用户在模型上指定的方向和特征约束的方向为边界条件，光顺去噪后得到方向场；

3)根据方向场的旋度约束能量和用户提出的边界约束能量、各向异性约束能量和光顺约束能量，优化得到各向异性的密度场。

3.根据权利要求1所述的一种基于波动方程的各向异性的四边形网格生成方法，其特征在于所述的从二维驻波中形成模型的四边形剖分的步骤为：

1)找出二维驻波场的极大点、极小点和鞍点；

2)从每个二维驻波场的鞍点出发，沿着二维驻波场梯度上升/下降最快的方向在模型表面上前进，直到遇到二维驻波场的极大/极小点，得到剖分图，采用剖分图的对偶图形成模型的四边形剖分。

4.根据权利要求1所述的一种基于波动方程的各向异性的四边形网格生成方法，其特征在于所述的在四边形剖分的基础上构造各向异性的全局参数化，然后得到最终的四边形化网格的步骤为：

1)为剖分得到的每个四边形区域确定参数化映射空间；

2)根据参数化映射空间，为任意相邻的四边形区域确定参数化坐标之间的转换函数；

3)利用相邻四边形区域参数化坐标之间的转换函数，构造各向异性的全局参数化；

4)根据全局参数化的结果生成最终的四边形化网格。

Images