CN102254351B - 一种基于三维光顺交叉标架场的六面体网格生成方法 - Google Patents

Abstract

本发明公开了一种基于三维光顺交叉标架场的六面体网格生成方法。包括如下步骤：1)根据用户需求和模型特点，利用Spherical Harmonics表示生成满足相应表面对齐约束和光顺约束的三维交叉标架场；2)根据用户需求和标架场，在体网格上构造三维驻波参数化；3)根据参数化结果抽取每个六面体单元的顶点和中心点坐标，建立它们的拓扑连接关系，得到最终的六面体网格。本发明提出了一种利用基于Spherical Harmonic表示的三维光顺交叉标架场生成高质量的六面体网格，同时能够对六面体网格的尺寸，特征对齐等需求进行灵活的控制和优化。

Description

一种基于三维光顺交叉标架场的六面体网格生成方法

技术领域

本发明涉及六面体网格生成方法，尤其涉及一种基于三维光顺交叉标架场的六面体网格生成方法。

背景技术

六面体网格在许多领域有着广泛的应用，比如几何设计，有限元模拟，B样条重建，三维纹理贴图等。相比起四面体网格，六面体网格可以更精确的获取模型的对称性等特征，从而使设计师的工作更加方便直观；但自动生成满足特征约束，尺寸以及规整度控制的六面体网格仍然具有很大的挑战。

目前六面体网格生成方法可以分为半自动方法和自动方法，半自动方法主要有下面两种：

1.设计人员采用人工或半自动的方法将三维模型分解成简单的子模型，然后将这个简单的子模型分解为六面体网格，这种方法对体网格的拓扑几何形状要求较高[可参考POCA，X.2009 Paving the path towards automatic hexahedralmeshgeneration.PhD thesis，Ph.D.Program in Applied Mathematics，UPC，Barcelona]。

2.多块方法，这种方法依赖于用户输入来生成一个非常粗糙的六面体网格，然后再细分成小的六面体[可参考POINTWISE，2009，Gridgen-reliable cfdmeshing.http://www.pointwise.com/gridgen/]。

自动生成六面体的方法有网格方法，分解方法等，但这些方法并不能保证生成高质量，保持边界特征的六面体网格[可参考TAUTGES，T.J.2001.Thegeneration of hexahedral meshes for assembly geometry：survey and progress.International Journal for Numerical Methods in Engineering 50，pp.2617-2642]。

最近，网格四边形化工作吸引了许多研究人员的注意并取得了很大的进展，其中利用二维对称方向场控制四边形网格的方向的方法已被研究人员广泛采用。但对三维方向场(标架场)的研究工作并没有将其中的对称性考虑进去。

发明内容

针对背景技术的不足，本发明的目的在于提出了一种基于三维交叉标架场的六面体网格生成方法。

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

基于三维交叉标架场的六面体网格生成方法包括如下步骤：

1)根据用户需求和模型特点，利用Spherical Harmonics表示生成满足相应表面对齐约束和光顺约束的三维交叉标架场；

2)根据生成的标架场和用户给定的边长，在体网格上计算三维驻波参数化；

3)从三维驻波参数化结果中抽取每个六面体单元，建立它们之间的拓扑连接关系得到最终的六面体网格。

所述的根据用户需求和模型特点，利用Spherical Harmonics表示生成满足相应表面对齐约束和光顺约束的三维交叉标架场的步骤为：

1)离散化组成模型的四面体网格的三维交叉标架场，并用SphericalHarmonics表示；

2)以模型表面网格的三维交叉标架场作为边界条件，对内部结点的三维交叉标架场梯度进行积分得到光顺能量；

3)根据光顺能量和特征约束能量优化得到光顺的三维交叉标架场。

所述的根据生成的三维标架场和用户给定的边长，在体网格上计算三维驻波参数化的步骤为：

对于体网格中的一个点p，它的三维驻波相位角用该点的标架所关联[d₁(p)，d₂(p)，d₃(p)]的ZYZ欧拉角Θ(p)＝[θ₁(p)，θ₂(p)，θ₃(p)]表示；对于它周围的相邻点q，顶点(p，q)之间的三维驻波振幅关系定义如下：

其中μ为全局波长，设定为用户输入边长的二倍，π为圆周率常数，k取值1，2，3表示x，y，z三个轴，<>为内积计算符。然后根据该振幅关系以及表面边界对齐约束能量在体网格内部构造三维驻波。

所述的从参数化结果中抽取六面体单元，建立拓扑连接关系生成六面体网格的步骤为：

1)找出三维驻波场中相位角为[n₁π，n₂π，n₃π]，n_i为整数的点作为六面体单元的顶点以及找出相位角为[(n₁+0.5)π，(n₂+0.5)π，(n₃+0.5)π]的点作为六面体单元的中心点；

2)根据找出的中心点和顶点构造满足质量要求的六面体单元并形成它们的拓扑连接关系；

3)对体网格内部未被六面体填充的区域，用棱柱和四面体填充；

4)细分棱柱和四面体形成完全的六面体网格细分得到最终的六面体网格。

与背景技术相比，本发明具有的有益效果是：

1)采用Spherical Harmonics表示三维交叉标架场，避免了整数规划问题而且计算更加高效；

2)可以对六面体网格的形状，密度，方向，特征对齐的需求进行灵活的控制；可以自动的生成高质量的六面体网格。

附图说明

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

图2是根据三维标架场和用户给定的边长，求解得到的满足特征对齐约束的三维驻波参数化结果。

图3是最终生成的六面体网格。

图4是生成的六面体网格的内部区域。

具体实施方式

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

如图1所示，基于三维交叉标架场的波动方程的六面体网格生成方法包括如下步骤：

1根据用户需求和模型特点，利用Spherical Harmonics表示生成满足相应表面对齐约束和光顺约束的三维光顺标架场；

2根据三维标架场和用户给定的边长，在体网格上计算三维驻波参数化；

3从参数化结果中抽取六面体单元的顶点和中心点位置，建立它们之间的拓扑连接关系，用棱柱和四面体填充内部的丢失区域，细分这些非六面体区域得到最终的六面体网格；

所述的根据用户需求和模型特点，生成满足相应的表面对齐约束和光顺约束的三维光顺标架场，并以Spherical Harmonics进行表示的步骤为：

1)离散化表示组成模型的四面体网格的三维标架场，并用SphericalHarmonics进行表示；

2)以模型表面网格的标架场作为边界条件，对内部点的Spherical Harmonics系数梯度进行积分得到光顺能量；

3)根据光顺能量和特征约束能量优化得到光顺的三维标架场。

引导六面体方向的三维交叉标架场是必须满足下面三个性质：

1)立方对称性：沿标架的任意轴旋转π/2倍的角度得到的仍是同一个标架；

2)光顺性：整个标架场中的变化应尽量光顺；

3)法向对齐：边界上的标架的某个轴应该和该处的法向对齐，从而使六面体的网格和整个边界对齐。

根据以上要求，定义该三维标架场

R为应用旋转矩阵所获得的所有标架的等价类，其中r₁，r₂，r₃∈{0，1，2，3}，

为沿标架的*轴旋转r×π/2角后得到的旋转矩阵；对于两个不同的标架F(R_a)和F(R_b)，它们之间的差异性度量定义如下：

${\&Integral;}_{s^2}{(f(s;R_{\mathrm{\&alpha;}})-f(s;R_b))}^2\mathrm{ds}$

其中 $f(s;R)=t_x^2{t}_y^2+t_y^2{t}_z^2+t_z^2{t}_x^2,{[t_x,t_y,t_z]}^T=t=R^T$ 为F(R)定义在单位球上的描述符。

利用Spherical Harmonics对f(s；R)进行表达可以更高效的计算。对于以ZYZ欧拉角给定的旋转矩阵R(α，β，γ)＝Z(γ)Y(β)Z(α)，f(s；R)所对应的SphericalHarmonics系数为：

$f\left(R\right)=(\left(\hat{Z}\left(r\right)\hat{X}(-\mathrm{\&pi;}/2)\hat{Z}\left(\mathrm{\&beta;}\right)\hat{X}(\mathrm{\&pi;}/2)\hat{Z}\left(\mathrm{\&alpha;}\right)\right)\hat{f}\left(I\right)=\hat{R}(\mathrm{\&alpha;},\mathrm{\&beta;},\mathrm{\&gamma;})\hat{f}\left(I\right)$

其中

为f(s；I)的Spherical Harmonics系数，I为单位阵，

为旋转矩阵所对应的Spherical Harmonics系数的变换。

根据上面的定义，两个不同标架的差异度量可以用Spherical Harmonics表示为：

${\&Integral;}_{s^2}{(f(s;R_a)-f(s;R_b))}^2\mathrm{ds}={\left|\right|\widehat{R_a}\hat{f}\left(I\right)-\widehat{R_b}\hat{f}\left(I\right)\left|\right|}^2$

f(s；I)可以无损的投射到Spherical Harmonics在band 0和band 4的3个spherical harmonics基上：

$f(s;I)=s_x^2{s}_v^2+s_v^2{s}_z^2+s_z^2{s}_x^2=-\frac{2\sqrt{\mathrm{\&pi;}}}{15}(Y_4^0+\sqrt{\frac{5}{7}{Y}_4^4+16\sqrt{\mathrm{\&pi;}}{Y}_0^0})$

其中

为位于band上阶为m在Spherical Harmonics基。

由于band 0部分

总是常量对差异度量并没有任何贡献所以可以忽略，从而得到：

$f(s;I)=\sqrt{7}{Y}_4^0+\sqrt{5}{Y}_4^4$

此外，对于表面边界对齐问题，采用Spherical Harmonics可以将该问题简化为：

$\hat{f}\left(R_{n\&RightArrow;z}\right)\left[4\right]={\hat{R}}_{n\&RightArrow;z}\hat{f}\left(I\right)\left[4\right]=\sqrt{7}$

其中R_n→z表示将n旋转到z轴。

在离散化表示中，将三维交叉标架场定义在四面体TET_j每个面的中心点p，并用ZYZ欧拉角关联Θ(p)＝[θ₁(p)，θ₂(p)，θ₃(p)]。令表示第k个系数的梯度值，根据上面的离散化表示可以将光顺度量定义为：

$E_s\left({\mathrm{TET}}_j\right)=E_s\left(p\right)=\underset{k=0}{\overset{8}{\mathrm{\&Sigma;}}}\&dtri;\hat{f}{\left[k\right]}^2,$ p∈TET_j

对所有四面体计算上面的光顺能量得到总的光顺能量为：

$E_s=\underset{\mathrm{\&Omega;}}{\&Integral;}{E}_s\left(p\right)\mathrm{dp}=\underset{{\mathrm{TET}}_j}{\mathrm{\&Sigma;}}\mathrm{volume}\left({\mathrm{TET}}_j\right)E_s\left({\mathrm{TET}}_j\right)$

其中volume(TET_j)表示第j个四面体在体积。

在这种离散化表示下，表面边界对齐约束能量可以近似为：

$E_{\mathrm{\&alpha;}}=\underset{S}{\&Integral;}{\hat{R}}_{n\&RightArrow;z}\hat{f}\left(p\right)\left[4\right]-{\sqrt{7}}^2\mathrm{dp}\&ap;\underset{{\mathrm{TET}}_i}{\mathrm{\&Sigma;}}\mathrm{area}\left({\mathrm{TET}}_i\right)E_s\left({\mathrm{TET}}_i\right){\hat{R}}_{n_i\&RightArrow;z}\hat{f}\left(p_i\right)\left[4\right]-{\sqrt{7}}^2$

其中area(TET_j)表示第j个四面体在表面面积。

$E_f=\frac{E_s}{\mathrm{volume}\left(\mathrm{\&Omega;}\right)}+w_a\frac{E_a}{\mathrm{area}\left(S_s\right)}$

把上面的能量结合起来，优化如下的总能量，得到三维光顺交叉标架场：

其中w_a为控制表面边界对齐约束的相关权值，volume(Ω)和area(S)分别表示四面体网格在体积和表面积。

所述的根据三维标架场和用户给定的边长，在体网格上计算三维驻波参数化的步骤为：

对于体网格中的一个点p，它的三维驻波相位角可以用该点的标架[d₁(p)，d₂(p)，d₃(p)]所关联的ZYZ欧拉角Θ(p)＝[θ₁(p)，θ₂(p)，θ₃(p)]表示。对于它周围的相邻点q，顶点(p，q)之间的三维驻波振幅关系定义如下：

$f(p,q)={\mathrm{\&Pi;}}_{k=1}^3\cos ({\mathrm{\&theta;}}_k\left(p\right)+{\mathrm{\&phi;}}_k(p,q))$

其中μ为全局波长，设定为用户输入边长的二倍，π为圆周率常数，k取值1，2，3表示x，y，z三个轴，<>为内积计算符。然后根据该振幅关系以及表面边界对齐约束能量在体网格内部构造三维驻波。

根据上面的定义，可以得到每个四面体的驻波能量：

$E_w\left({\mathrm{TET}}_j\right)=\underset{q\&Element;\left\{p^i\right\}\&cup;\left\{p_i\right\}}{\mathrm{\&Sigma;}}\left(\underset{p_{a,},p_b\&Element;\left\{p_i\right\}}{\mathrm{\&Sigma;}}{(f(p_a,q)-f(p_b,q))}^2\right)$

以及总的驻波能量 $E_w=\underset{\mathrm{TET}}{\mathrm{\&Sigma;}}\mathrm{volume}\left({\mathrm{TET}}_j\right)E_w\left({\mathrm{TET}}_j\right)$

边界对齐要求生成的六面体的表面面片和给定模型的表面对齐，这实际上等价于要求：cos(θ_k(p))＝±1。

具体求解是不直接优化变量θ，而是采用Zhang，M.，HUANG，J.，LIU，X.，Bao，H.2010.A wave-based anisotropic quadrangulation method.ACM Trans.Graph.

29(July)，118：1-118：8中的方法，用表示三维驻波状态的变量

替代相位角，其中：

$\widehat{\mathrm{\&Theta;}}\left[1\right]=\mathrm{ccc},\widehat{\mathrm{\&Theta;}}\left[2\right]=\mathrm{ccs},\widehat{\mathrm{\&Theta;}}\left[3\right]=\csc ,\&CenterDot;\&CenterDot;\&CenterDot;,\widehat{\mathrm{\&Theta;}}\left[8\right]=\mathrm{sss},$ c＝cos(θ_k)，s＝sin(θ_k)

$E_{\mathrm{valid}}={({\widehat{\mathrm{\&Theta;}}}^2-1)}^2+\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}\underset{j=i+1}{\overset{4}{\mathrm{\&Sigma;}}}{\left(\widehat{\mathrm{\&Theta;}}\right[i\left]\widehat{\mathrm{\&Theta;}}\right[9-i]-\widehat{\mathrm{\&Theta;}}[j\left]\widehat{\mathrm{\&Theta;}}\right[9-j\left]\right)}^2$

此外为了保证求得的是有效的三维驻波，仍需要下面额外约束：

最终优化下面的能量求得驻波：

其中E_b为表面边界对齐约束能量，w_w，w_b，w_valid为各自能量系数。

所述的从驻波参数化结果中抽取每个六面体网格的顶点和中心点位置，并建立它们之间的拓扑连接关系的步骤为：

1)从驻波参数化结果中抽取六面体单元的顶点和中心点；

2)根据抽取的顶点和中心点生成具有拓扑约束的六面体单元；

3)用棱柱和四面体填充未被六面体单元填充的内部区域；

4)细分棱柱和四面体形成完全的六面体网格。

对于每个六面体单元的顶点来说，它应该位于相位角为[n₁π，n₂π，n₃π]，n_i为整数的位置，同样每个六面体单元的中心点应该定位在

[(n₁+0.5)π，(n₂+0.5)π，(n₃+0.5)π]的位置。但由于数值原因及波长选择对驻波参数化计算的影响，得到的结果会有误差。在这里采用一种基于mean shift[可参考CHENG，Y.1995.Mean shift，mode seeking，and clustering.IEEE Trans.Pattern Anal.Mach.Intell.17(August)，790-799.]的算法计算这些点的位置。

对于模型中的每个结点p，先按下面的方程得到它周围的一个点集Q，

Q_vertex(p)＝{q：||q-p||＜μ/2|and θ_k(p)+φ_k(p，q)＝n_kπ，k＝{1，2，3}}

然后迭代更新p点的位置以及合并位置非常相近的结点。具体上在第k次迭代式，按下面的公式更新

点位置到

${\hat{p}}_i^k=\left({\mathrm{\&Sigma;}}_{p_i^k\&Element;N(p_i^k,r_i^k)}{w}_i^k{p}_i^k\right)/\left({\mathrm{\&Sigma;}}_{p_i^k\&Element;N(p_i^k,r_i^k)}{w}_i^k\right)$

其中

表示

半径内的

点的所有邻居，

为计算得到的权值。

然后按下面的策略将位置非常近的点合并为一个。如果某个点p附近点满足

在抽取出所有的顶点和中心点之后，按下面的方法生成每个六面体单元：对每个中心点找到它所对应的八个最优点顶点形成一个候选的六面体单元，同时对每个六面体单元定义如下的质量度量：

$\mathrm{\&sigma;}=\frac{\mathrm{\&eta;}}{{(\mathrm{\&mu;}/2)}^3}+\frac{{(\mathrm{\&mu;}/2)}^3}{\mathrm{\&eta;}}+{\mathrm{\&xi;}}^3+{\mathrm{\&psi;}}^2$

其中η为这个六面体单元的体积，μ为用户输入的六面体边长，ξ为该六面体单元边长，ψ为六面体表面四边形角度；然后按质量从好到差点顺序生成每个六面体单元，如果加入某个候选单元后导致了不正确的拓扑，则移除这个候选单元选取下一个。

按照上面的策略形成的六面体网格内部经常会有一个未被六面体但与完全填充的区域。用棱柱或四面体填充这些部分，然后再将这些单元细分为完全的六面体网格。

实施例：

本发明提出一种基于三维光顺交叉标架场的六面体网格生成方法，流程图如图1所示，下面介绍处理一个实际模型的例子：

1)根据网格特征和用户需求，利用Spherical Harmonics离散化表示三维交叉标架场，通过优化光顺能量和表面边界对齐约束能量，生成三维光顺交叉标架场；

2)根据三维交叉标架场和用户输入边长，通过优化驻波能量方程，计算体网格内部的三维驻波参数化，见图2；

3)从驻波参数化结果中抽取每个六面体网格的顶点和中心点位置，建立它们之间的拓扑连接关系，细分棱柱和四面体得到最终的六面体网格见图3；生成的六面体网格的内部区域见图4。

Claims (1)

1.一种基于三维光顺交叉标架场的六面体网格生成方法，包含如下步骤：

1)根据用户需求和模型特点，利用Spherical Harmonics表示生成满足相应表面对齐约束和光顺约束的三维交叉标架场；

2)根据生成的标架场和用户给定的边长，在体网格上计算三维驻波参数化；

3)从三维驻波参数化结果中抽取每个六面体单元，建立它们之间的拓扑连接关系得到最终的六面体网格；其特征在于：

所述的根据用户需求和模型特点，利用Spherical Harmonics表示生成满足相应表面对齐约束和光顺约束的三维交叉标架场的步骤为：

1.1)离散化组成模型的四面体网格的三维交叉标架场，并用SphericalHarmonics表示；

1.2)以模型表面网格的三维交叉标架场作为边界条件，对内部结点的三维交叉标架场梯度进行积分得到光顺能量；

1.3)根据光顺能量和特征约束能量优化得到光顺的三维交叉标架场；

所述的根据生成的三维标架场和用户给定的边长，在体网格上计算三维驻波参数化的步骤为：

对于体网格中的一个点p，它的三维驻波相位角用该点的标架[d₁(p)，d₂(p)，d₃(p)]所关联的ZYZ欧拉角Θ(p)＝[θ₁(p)，θ₂(p)，θ₃(p)]表示；对于它周围的相邻点q，顶点(p，q)之间的三维驻波振幅关系定义如下：

其中μ为全局波长，设定为用户输入边长的二倍，π为圆周率常数，k取值1，2，3表示x，y，z三个轴，<>为内积计算符；然后根据该振幅关系以及表面边界对齐约束能量在体网格内部构造三维驻波；

所述的从三维驻波参数化结果中抽取六面体单元，建立拓扑连接关系生成六面体网格的步骤为：

2.1)找出三维驻波场中相位角为[n₁π，n₂π，n₃π]，n_i为整数的点作为六面体单元的顶点，以及找出相位角为[(n₁+0.5)π，(n₂+0.5)π，(n₃+0.5)π]的点作为六面体单元的中心点；

2.2)根据找出的中心点和顶点构造满足质量要求的六面体单元并形成它们的拓扑连接关系；

2.3)对体网格内部未被六面体填充的区域，用棱柱和四面体填充；

2.4)细分棱柱和四面体形成完全的六面体网格细分得到最终的六面体网格。

Images