一种基于共形几何的数模表面网格生成方法与流程

本发明涉及计算共形几何技术和三维图像处理领域，具体涉及一种三维图像处理中基于共形几何的数模表面网格生成方法。

背景技术：

随着计算机技术、计算机辅助设计(cad)技术和计算流体动力学(cfd)的迅猛发展，在工程制图、视频游戏、科学研究、文物复原和生物医学等领域中三维模型的数量呈爆炸式增长。在不同的领域都使用了大量数模处理软件和数模文件，这给三维模型提供了良好的发展机遇，同时航空航天等诸多应用领域实际问题的复杂性又对计算流体动力学提出了严峻的挑战，其中一个重要问题就是如何自动生成绕复杂外形的网格。据文献介绍，网格生成大致占整个计算周期的60％～70％的时间，而复杂外形的表面描述及表面网格生成一般要占网格生成时间的70％～80％。由此可见，复杂外形的表面描述及表面网格生成是阻碍三维模型持续纵向发展的“瓶颈”问题之一。

计算共形几何是一个交叉学科，根植于纯数学领域，主要研究曲面上的共形结构，在拓扑与几何之间寻找关联及平衡，挖掘曲面潜在且内蕴的特征，来更好地理解曲面，用于图形分析及处理。共形几何在数字几何处理中扮演了一个重要角色，尤其是奠定了曲面参数化的理论基础，并提供了严格的算法。计算共形几何的核心是共形结构，具有共形结构的曲面是黎曼面。共形结构的研究方法有很多种，复分析、代数几何和微分几何。近来，随着三维扫描技术的发展，计算能力的提升，以及数学理论的进一步发展，人们已经把计算共形几何理论和算法从平面区域推广到了具有任意拓扑的度量曲面。

网格生成算法的种类是多种多样的，根据算法的几何分析域的特点，大致将各种算法归纳为：映射法、前沿推进法(advancingfrontmethod)、delaunay三角化法、拓扑分解法和几何分解法。其中映射法出现于20世纪70年代，是最早的网格生成方法，它是根据形体边界的参数方程定义适当的映射函数，将目标区域映射到参数空间中形成规则参数域，对规则参数域进行网格剖分，将参数域的网格(二维是正方形，三维是立方体)逆映射回原来的非规则区域上，从而生成实际的物理空间网格，优点在于计算效率高，实现简单，网格疏密度可控性强，可直接生成四边形、六面体等高质量的单元网格，能够用于曲面网格的生成，但是映射法对于形状较为复杂的形体适应性差，需要将复杂形体事先分解成若干形状简单的子区域，子区域分解繁琐费时，人工交互多，自动化程度低，难以实现全自动化。前沿推进法用于平面区域三角形网格自动生成，已经取得了很好的效果并推广到了二维自适应网格的生成和任意形状区域的三维四面体网格生成中，前沿推进法只需要给定区域边界便可生成网格，生成的网格具有可重复性，此外，前沿推进法还有适应性广，单元质量高等优点，缺点主要是由于没有较强的理论基础，不能保证该算法的收敛性。delaunay三角剖分是目前最为流行的且通用的全自动网格生成方法之一，最大最小角特性和外接圆特性保证能生成较高质量的单元格，随着delaunay技术的不断发展，现如今已经出现了许多行之有效的算法，一般可分为以下三大类：以bowyer和green、sibsons为代表的计算voronoi图方法；以watson为代表的空外接圆法；以lawson为代表的对角线交换算法。直接计算voronoi图的方法通常是比较复杂的，所需内存大，而且计算效率也比较低，这类算法现在很少被采用；lawson算法特别适用于二维delaunay三角化，它不存在退化现象，对约束情况同样适用，计算效率高，但是不能直接推广到三维情况；watson算法概念简单，易于编程实现，也能够实现约束三角化，因此应用频度最广，但该法也有不足，即可能出现内角小的薄片单元。拓扑分解法最初应用于二维平面问题，现如今已经推广到三维空间曲面中，拓扑分解法实质上就是在保证拓扑结构连续性的基础上，每次切割下一个三角形单元，使得三角形单元尽可能地接近等边三角形，而不考虑形体的具体形状，当待剖分区域已经形成初始网格之后，再用对半分割、中心分割、对角线置换等手段对单元格进行逐步的细化。拓扑分解法原理简单，实现起来比较方便，但是它只从形体的拓扑因素入手，过分依赖于形体的拓扑结构，而不考虑几何因素，因此很难获得质量良好的单元格，也难以对包含曲面的三维形体进行处理。几何分解法是一种递归对半分割法，几何分解法较多地考虑了形体的几何特征，是节点和单元格同时产生的方法。几何分解法的优点是从待剖分区域的几何因素着手划分单元，每次都可以得到形状良好的局部结果，所以，该方法通常可形成质量很高的单元格，它的缺点是对复杂区域和不规则形体首先要分割成简单的子区域，人工干预比较多，执行起来效率低，除此之外，逐层循环分割法只能做到每次取出的单元形状良好，但不能保证最终剩余的部分是合理的，即难以把空白地带解决好。

针对上述问题，如何结合共形几何的技术建立一套针对任意复杂数模的通用的表面网格生成方法，就显得尤为重要。

技术实现要素：

本发明提出了一种基于共形几何的数模表面网格生成方法，通过使用共形几何技术来解决三维模型表面网格的生成问题。

本发明采用的技术方案为：一种基于共形几何的数模表面网格生成方法，包括以下步骤：

(1)根据输入的三维模型数据判断三维模型表面是否可定向，计算曲面的亏格数目和边界数目；

(2)根据计算得到的拓扑结构选择合适的共形映射函数，根据曲面单值化定理，将三维数模表面共形映射到球面，欧式平面和二维双曲空间中的一种；

(3)根据所选的共形映射函数将三维模型共形映射到二维参数域上；

(4)根据共形映射的结果，在参数域上进行有限元网格剖分，按需生成三角形单元网格或四边形单元网格；

(5)根据共形映射过程中保持的角度，将重新剖分后的网格逆映射回三维空间域中，生成带表面网格的三维模型。

在步骤2中，将三维模型理解成连通的曲面进行处理，曲面的拓扑结构由是否可定向，环柄数目(亏格)和边界数目所决定。

在步骤3中，根据步骤2中判断的拓扑结构，将分割好的曲面共形地变换到三种标准空间中的一种：球面，欧式平面和二维双曲空间。

在步骤4中，网格剖分有两个重要的步骤：在待剖分区域内合理分布网格节点，即布点；将网格节点有效地连接起来，形成三角形或四边形网格单元，即单元格生成，其中三角网格生成采用delaunay三角化法。

在步骤5中，按照参数域到空间域的映射函数将参数域的网格映射成空间域的网格，并保持节点间的拓扑关系不变，共形变换将三维变成二维，角度保持，但是局部面积发生变化。局部面元的变化率和曲面的平均曲率完整地保持了原来曲面所有的几何信息，根据面元变化率和平均曲率来完全重建原始曲面。

所述三维模型均可视作连通的曲面，对三维模型进行参数化映射可视作对曲面进行参数化映射。整个三维模型表面网格生成过程在数学上等价于求从曲面空间域到平面参数域的一个光滑双射。

本发明与现有技术相比的优点在于：相比直接在三维数模表面进行网格剖分，参数化区域上进行的网格剖分操作更简单，计算复杂性更低，算法更易于实现。通过上述步骤，将三维数模表面网格生成问题简化成求三维空间域到二维参数域的光滑双射和二维参数域网格生成问题。通过该步骤生成所生成的三维模型表面网格质量高，精度好。

附图说明

图1是基于共形几何的数模表面网格生成具体实现示意图；

图2是根据本发明实现一个具体实例示意图；

图3是基于有限元网格剖分过程示意图。

具体实施方式

下面结合附图和具体实施方式对本发明进行描述。其中图1描述了基于共形几何的数模表面网格生成具体过程。图3描述了有限元网格剖分具体实现过程。

具体的实现步骤：

(1)根据所给三维模型，将其转换成可存入半边数据结构的.off或者.obj格式。首先，在半边数据结构中对模型进行边界识别，在遍历的过程中把只有一个邻接面的边当作边界边，如果边界边之间有闭环，则可当作存在一条边界；然后对曲面进行亏格识别，把能在曲面上画出而又不把曲面分割开的互不相交简单闭曲线的最多个数当作该模型曲面亏格的个数。统计曲面的边界个数和亏格个数，存为曲面的拓扑信息进行下一步。

对于亏格为0的曲面，可以由调和映射来计算其共形结构。假设两个0亏格网格曲面m1，m2,h:m1→m2表示曲面的度为1的映射，接下来最小化调和能量e(h)，

h的拉普拉斯方程可以简化为：

δ^plh＝(δ^plh0,δ^plh1,δ^plh2)

如果h是调和的，那么δ^plh的切向分量为0。定义投影算子：

其中为张量积，i为单位矩阵。当且仅当h满足：

时h是调和的，其中n是m2的法线。

为了确保过程收敛到唯一解必须添加额外约束。这里迫使表面的质心位于原点：

是m1的面积。这个约束能够保证在旋转过程中解的唯一性。在这个约束条件下构建偏微分方程：

此时，h的稳态解就是m1到m2的共形映射。该方程能够有效地被迭代求解。

对于亏格为1的曲面，可以通过全纯微分来计算，该方法基于hodge理论，并且对于带有任意亏格的闭合或者非闭合曲面都是行之有效的。该方法的核心思想就是找到一个全纯1形式的完备基，算法大致思想如下：

输入网格曲面m。

输出全纯微分基形式

计算同源群基b＝{e1,e2,…,e2g}；

计算对偶上同调群基ω＝{w1,w2,…,w2g}；

把每个wi分散成调和1形式ζi；

计算每个ζi的共轭，表示为*ζi。构造全纯1形式

对于任意亏格的曲面，曲面的共形结构可以通过离散ricci流来计算，具体步骤如下：

假设给定带度量的曲面间的光滑映射(s1,g1)→(s2,g2)，局部上，非线性的映射可以用其一阶近似，切平面间的线性映射来逼近，所谓的导映射：

会带来各种畸变，这些畸变的测量需要用到黎曼度量。线性映射的分类可以利用上述线性映射分类的想法，例如：如果对于任意一点p∈s1，都是等距(或者保角，保面积，拟共形)，则整体上是等距(或者保角，保面积，拟共形)。

共形度量的计算和嵌入过程如下：

通过黎曼度量可以计算长度、角度、面积和曲面s上的微分算子。在参数化过程中，黎曼度量被共形改变着，除角度外的其他相关量也相应的改变着。直到度量g变化为那么此时高斯曲率将变为：

其中δ为原始度量g的拉普拉斯算子；此时的测地曲率将变为：

根据gauss-bonnet公式，总曲率和黎曼度量的选取无关，可以证明，存在唯一的共形度量使得目标曲率是常数而测地曲率为零。

通过曲率控制来进行参数化的过程，其实就是寻找合适的共形度量e^2ug，使之诱导的曲率与预设之目标曲率一致即可。

在光滑情形下，共形形变将无穷小圆映射到无穷小圆并保持圆与圆之间的交角，这样做可以保证共形映射的保形性和保角性。离散情形下，通过circlepacking方法的引入在网格的每个顶点处设置一个半径为γi的圆，不同的圆互相相交。

在欧氏背景中，共形形变在保证圆与圆交角不变的情况下，通过改变半径达到参数化的目的。离散的ricci曲率流有着和光滑情形ricci流方程类似的形式：

在球面或双曲几何背景下，使用ricci曲率流来计算球面度量或者双曲度量时，过程和欧氏情形下一样，只要在计算几何量时，使用相应的球面几何或者双曲几何即可。

假设双曲圆(c,r)以c为圆心以r为半径是欧氏圆(c,r)。满足：

其中u＝(e^r-1)/(e^r+1)。

接下来要先将第一个面[v0,v1,v2]嵌入到双曲空间中，参数坐标分别为：

把与第一个面相邻的面fijk嵌入其中。假设顶点vi,vj已经在双曲平面中了，那么第三个顶点vk即为以vi,vj为顶点的两个圆的交点。以同样的方式，将其他与fijk相邻的面都压入堆栈，按顺序嵌入到双曲平面中。双曲空间中的嵌入方式，和欧氏空间的嵌入方式基本一样。同样，在球面几何下嵌入的方法，也和欧氏空间及双曲空间下一样。

直接变分法基本思路如下：用分片线性映射来表示从三角网格到平面区域的映射，m→d。固定一个三角形，将其铺到平面上，用vi,vj,vk来表示顶点的平面坐标，用来表示顶点像的平面坐标，则线性映射矩阵可以直接写出：

如果eijk为0，则限制在此三角形上的线性映射为保角变换。整个映射和等距变换的差距定义为：

这里sijk代表相应三角形的面积。然后用非线性优化方法来优化这个能量，并将所得结果作为全局保角变换的近似：

如果目的是寻求保面积变换，则：

eijk:＝(λ1λ2-1)²＝(det(aijk)-1)²

如果目的在于等距变换，定义局部能量：

eijk:＝(λ1-1)²+(λ2-1)²＝(λ1+λ2)²-2(λ1+λ2)-2λ1λ2+2

等价地：

eijk:＝tr(aijk)²-2tr(aijk)-2det(aijk)+2

下面是离散ricci曲率流算法的具体步骤：

1.通过最小化能量方程：

为网格的边及顶点分别赋初始边长及半径，可赋初始值：

2.使用余弦定理计算当前网格的黎曼度量即边长：

3.计算网格的所有内角：

4.为每个顶点计算离散的高斯曲率ki：

5.使用如下公式更新顶点半径γi：

其中ε为步长；

6.单位化网格顶点处的半径，使其乘积为1：

7.检查当前曲率ki与目标曲率之差，选取其中差距最大的：

如果误差已经小于阈值，那么就可以停止计算，否则调回第二步。整个收敛过程用时很短。选取一个顶点vi，误差曲线为一条指数曲线。

(2)网格剖分的两个主要步骤分别是布点和单元格生成。网格的节点可以分为边界节点和内部节点两类，以下为两种节点的布点方法：

1.边界节点的生成方法

如果给定区域边界是以直线、曲线的方程形式给出，设定三角形单元的剖分尺寸为d，以d来等分区域边界，从而可以计算出边界节点的坐标；如果区域带有内部约束，其内部特征约束可以作如下处理：如果是点，直接把点当成边界节点计入节点数组中；如果是直线段、圆弧、曲线等组成的不封闭的环，可将它们看成是有向曲线，并抽象成面积为零的封闭环，当作内边界环来处理。

2.内部节点的生成方法

对于内部节点的生成，在传统平行线布点的基础上，采用两组等距平行线交叉成60°角的布点方法，此方法所布点之间一般可形成正则三角形，可以极大提高网格的质量。

具体布点方法如下：

先求出区域的横坐标和纵坐标的最大值xmax、ymax及最小值xmin、ymin；再在平面区域的最小矩阵包容盒内布点，y＝ymax和y＝ymin这两条线是极限扫描线，在这两条线之间可以生成等距的横向扫描线，条数为：

m＝int[(ymax-ymin)/(length×sin60°)+0.5]

扫描线的方程为：

y＝yi＝ymin+(ymax-ymin)×i/mi＝0,1,2…m

扫描线上点的个数为：

n＝int[(xmax-xmin)/length+0.5]

扫描线上点x的坐标为：

其中i＝0,1,2,…m；j＝0,1,2,…n；k为自然数。

扫描线上相邻两点间的距离是length，相邻扫描线间的距离是length×60°，length/2为纵向相邻两点错开的距离，这样可使得内部节点容易形成等边三角形。

按上述方法布点之后，要删掉外边界以外和内边界以内的节点，还要去掉距离边界较近的节点以保证网格剖分的质量，剩余的节点就是内部节点。

平面区域三角网格生成采用delaunay三角剖分法。假定s为e^d空间的有限点集，点集s的delaunay三角化定义为满足以下条件的单纯形单元复形d(s)：

复形d(s)的0-单纯形组成的集合(即d(s)的所有顶点的集合)是s的子集；

复形的底空间是点集s的凸包ch(s)；

任意一个d-单纯形δt∈d(s)，|t|＝k+1，满足：任意q∈s-t，q在δt的外接球外。

在二维参数域中生成三角网格后，再按照参数域到空间域的映射函数将参数域的网格共形映射成空间域的网格，并保持节点间拓扑关系不变，得到的空间域的网格就是所求的曲面三角网格。