本发明属于几何处理领域中技术，提出了一种基于莫尔斯-斯梅尔复形和参数化的混合四边形化方法，在用户输入标架场导引下，本发明方法能够生成与输入标架场和特征线对齐的的高质量纯四边形网格。
背景技术：
四边形网格作为一种半规整网格具有相当广泛的应用范围,学术界对此进行了几十年的研究,取得了长足的进展。而用户对网格质量的要求也随着技术的发展逐渐提高。将三角形网格重网格化为四边形网格之前也有一些方法可以实现，但是有其局限性。
到目前为止，具有特征和任意引导框架区域的表面的全自动四边形仍然是一个挑战：快速且可靠地生成能与用户输入的标架场紧密对齐的纯四边形网格。传统的基于参数化的方法，依赖于输入标架场的奇异点的拓扑结构，需要利用混合整数优化得到参数化，然后通过提取参数化中的等高线来得到四边形网格。虽然这种一般方法很高效，但不能保证得到纯四边形网格。另一类方法基于莫尔斯-斯梅尔复形的四边形重网格化，该方法具有理论保证。然而，该方法需要足够密的三角网格作为输入，使其很难运用于实际应用中。
技术实现要素：
针对背景技术的不足，本发明的目的在于提供一种基于莫尔斯-斯梅尔复形和参数化的混合四边形化方法，将参数化和莫尔斯-斯梅尔复形相结合到网格四边形处理中，在保证得到四边形网格的同时计算效率高。
本发明首先根据用户输入的流形三角网格、标架场和特征线，优化得到周期性四维向量场，然后利用该周期性四维向量场，构建全局参数化，并从该全局参数化中抽取四边形网格，对于该全局参数化中不能直接抽取四边形网格的部分，使用莫尔斯-斯梅尔复形生成四边形网格，最终将两部分网格合成为最终的四边形网格。
为实现上述的目的，本发明采用的技术方案具体是如下步骤：
1)对于给定的一个流形三角网格M，所有顶点集合表示为V，三角面片集合表示为T，用户输入定义在三角面片t上的标架为Ft，三角面片t∈T，所有三角面片上的标架集合构成标架场F，对于流形三角网格M中的每个顶点v∈V，顶点v处的标架定义为Fv，从与顶点v相邻的三角面片中任意选取一个三角面片标架作为Fv，特征线和边界的所有顶点的集合为Vb，特征线角点和边界角点的集合为Vc；
所述的特征线是用户输入的边集和流形三角网格M的所有边界边的并集，并集中的每个元素为特征边。特征线角点是指多条特征边相交处的顶点集合和特征线的端点集合的并集。
2)在流形三角网格上定义周期性四维向量场Ψ(p)和以下向量场能量ε(Ψ)的公式：
其中，Ψ是所有顶点p的周期性四维向量Ψp的集合，即周期性四维向量场，Ψ＝{Ψp|p∈V}；是定义在顶点p处的四维向量，分别表示周期性四维向量Ψp的第零、第一、第二和第三分量，向量场能量ε(Ψ)代表周期性四维向量场Ψ的光滑程度；
然后求解能量公式得到光滑并且满足所有要求的周期性四维向量场Ψ，求解首先是通过求解特征值问题得到初值，然后从初值出发进行带惩罚项的非线性优化计算得到周期性四维向量场；
3)从优化得到的周期性四维向量场中，对每一个顶点p计算得到与之相容的局部参数化坐标
4)利用生成树将所有局部参数化坐标组合构成全局参数化坐标；将全局参数化坐标中满足无缝与特征对齐条件的区域作为平凡区域，从平凡区域抽取获得四边形网格QR；
其中，无缝条件指的是四边形网格QR中相邻的四边形的边能够重合，特征对齐条件指的是在用户定义的特征线处有对应的四边形边。
抽取获得平凡区域的四边形网格QR具体是从所有全局参数化坐标中选取坐标值为整数的点构成四边形集合QR。
5)取流形三角网格M中未被平凡区域的四边形网格QR覆盖的区域作为奇异区域，然后在奇异区域使用莫尔斯-斯梅尔复形抽取获得奇异区域的四边形网格QS；
6)通过合并平凡区域的四边形网格QR和奇异区域的四边形网格QS获得最终的四边形网格Q，完成混合四边形化。
所述步骤2)中向量场能量ε(Ψ)定义为：
在上式中，Nb-f表示与特征线和边界相连的边的集合，pq是该集合中的一个元素，epq是边pq的向量表示，μ用于调节边界和特征线处对齐程度的权重，|t|是三角面片t的面积，εpq表示边pq上的能量，Wsum表示权重归一化因子；
边pq上的能量εpq通过以下公式计算得到：
cc(a,b)＝cos(πa)cos(πb)，sc(a,b)＝sin(πa)cos(πb)，cs(a,b)＝cos(πa)sin(πb)，ss(a,b)＝sin(πa)sin(πb)
其中，W表示四维向量的旋转矩阵，kpt、kqt分别为顶点p和顶点q在三角面片t上的旋转系数，为旋转矩阵，是一个二维向量，定义为(a,b)T，a,b分别表示二维向量的第一分量和第二分量；
表示旋转矩阵的kqt次幂，表示旋转矩阵的kpt次幂，旋转矩阵采用以下公式计算：
其中，顶点p在三角面片t上的旋转系数kpt和顶点q在三角面片t上的旋转系数kqt采用以下公式计算：
kqt＝argmink||FqJ4-k-Ft||2,k∈{1,2,3,4}
kpt＝argmink||FpJ4-k-Ft||2,k∈{1,2,3,4}
其中，Fp表示顶点p的标架，Ft表示三角面片t的标架，标架实际上是建立在三角面片t上的2×2矩阵，k表示待优化系数，J表示逆时针旋转π/2的旋转矩阵。
所述步骤2)的求解过程中，首先求解以下公式的特征值问题：
HΨ*＝λminΨ*
其中，H是能量ε(Ψ)的Hessian矩阵，Ψ*表示周期性四维向量场Ψ的初值；上述特征值问题的方程求解中线性的特征对齐约束通过变量替换消除。
通过求解上式得到周期性四维向量场Ψ的初值Ψ*，然后利用带惩罚项的非线性优化采用以下公式求解得到最终的周期性四维向量场Ψ作为四维向量场：
其中，wr是控制惩罚项的权重，|V|表示顶点集合V中的顶点总数。
上式中，是惩罚项。
所述步骤3)中每一个顶点p的局部参数化坐标计算公式如下：
其中，分别表示周期性四维向量Ψp的第零、第一、第二和第三分量。
所述步骤4)中，利用生成树将所有局部参数化坐标组合构成全局参数化坐标，具体为：
4.1)以三角面片作为树的节点、以流形三角网格M的对偶边作为树的边构造生成树，对偶边是连接相邻三角面片的共同边，将生成树中的根节点三角面片记为tmin具体实施中随机选取一个三角面片为根节点；
4.2)按照以下方式从根节点出发构造全局参数化坐标，全局参数化坐标为一个二维向量：
4.2.1)根节点tmin的全局参数化坐标计算
设置根节点tmin处的第一个顶点p1的全局参数化坐标为然后对于根节点tmin处的其余两个顶点p2和p3，按照根节点tmin的三角面片逆时针序的采用以下公式计算两个顶点的全局参数化坐标：
其中，表示偶数集合，分别表示顶点pi的局部参数化坐标到根节点t全局参数化坐标的两个分量对应的整数平移量，表示三角面片tmin标架场的逆矩阵，表示以顶点pi为起点且以顶点pj为终点的边，表示根节点处顶点pi的旋转类型，使用下式计算：
4.2.2)从根节点tmin开始广度优先搜索方法遍历树的边，如果当前的树边是从父节点的三角面片s到子节点的三角面片t，其中顶点p和q是三角面片s和t的公共顶点，r是三角面片t上除p和q之外的第三个顶点，根据以下方程求解得到r的全局参数化坐标
其中，表示三角面片t中顶点p的全局参数化坐标，表示三角面片t中顶点q的全局参数化坐标，表示三角面片t中顶点r的全局参数化坐标；表示顶点r的局部参数化坐标，是三角面片t上的形变梯度，为三角面片t上顶点r的整数平移量，表示三角面片t的顶点r上从局部参数化坐标到全局参数化坐标对应的整数平移量的两个分量；是旋转矩阵，kt表示三角面片t上的旋转类型，根节点处的ks表示三角面片s上的旋转类型，遍历当前树边时该值已经被计算得到，kts表示从三角面片t到三角面片s的旋转类型，krt表示从顶点r到三角面片t的旋转类型，krt表示从顶点r到面片t的旋转类型，kt＝kts+ks；
上述kts和krt表示为：
kts＝argmink||FtJ4-k-Fs||2,k∈{1,2,3,4}
krt＝argmink||FrJ4-k-Ft||2,k∈{1,2,3,4}，
其中，Fr是顶点r的标架，Ft是三角面片t的标架，Fs是三角面片s上的标架。
6.根据权利要求1所述的一种基于莫尔斯-斯梅尔复形和参数化的混合四边形化方法，其特征在于：所述步骤5)中：
5.1)取流形三角网格M中未被平凡区域MR的四边形网格QR覆盖的区域作为奇异区域MS，即MS＝M/QR；
5.2)在奇异区域MS，利用周期性四维向量场的第零分量构建莫尔斯-斯梅尔复形，即抽取所有极大值点、极小值点和鞍点，从每个鞍点出发搜索最近的两个极大值点和两个极小值点构成一个四边形，最终形成奇异区域MS的四边形网格QS。
所述步骤6)具体是通过合并平凡区域和奇异区域的四边形网格QS∪QR获得最终的四边形网格Q。
本发明的有益效果是：
本发明方法能根据输入的流形三角网格，快速且鲁棒地生成与用户输入标架场相匹配、与用户输入特征线严格对齐的纯四边形网格。本发明以任意的非正交且各向异性的标架场作为输入，并能严格对齐用户给定的复杂特征线约束。本发明方法既有莫尔斯-斯梅尔复形方法的理论保证的特点，又具备全局参数化方法的效率优势，所以能从理论上保证得到纯四边形网格，且计算效率高。
附图说明
图1是本发明的流程图。
图2是实施例输入流形三角形网格的示意图。
图3是实施例输入的标架场示意图。
图4是实施例输入的用户约束示意图。
图5是实施例由全局参数化生成的平凡区域的四边形网格和奇异区域的莫尔斯-斯梅尔复形的结果示意图。
图6是实施例最终得到的四边形网格结果示意图。
具体实施方式
下面结合附图和实施例对本发明作进一步说明。
如图1所示，按照本发明发明内容完整方法实施的实施例如下：
1)对于给定的一个流形三角网格M，所有顶点集合表示为V，三角面片集合表示为T(如图2所示)，用户输入定义在三角面片t上的标架为Ft，三角面片t∈T，这些面片上的标架集合构成标架场F，对于M中的每个顶点v∈V，顶点v处的标架定义为Fv，从与顶点v相邻三角面片中任意选取一个三角面片标架作为Fv，特征线和边界的所有顶点的集合为Vb，特征线角点和边界角点的集合为Vc(标架场如图3所示，用户约束如图4所示)。本例中使用的模型是一个汽车外壳模型，具有42465个三角网格。
2)在流形三角网格上定义周期性四维向量场Ψ(p)，建立向量场能量ε(Ψ)的公式：
其中，向量场能量ε(Ψ)代表周期性四维向量场Ψ的光滑程度，ε(Ψ)定义为：
在本例中μ＝20，|t|是三角面片t的面积。
然后求解能量公式得到光滑并且满足所有要求的周期性四维向量场Ψ，求解首先是通过求解特征值问题得到初值，然后利用初值进行带惩罚项的非线性优化得到周期向量场。
首先求解以下公式的特征值问题：
HΨ*＝λminΨ*
其中，H是能量ε(Ψ)的Hessian矩阵，Ψ*表示周期性四维向量场Ψ的初值；上述特征值问题的方程求解中线性的特征对齐约束通过变量替换消除。
通过求解上式得到周期性四维向量场Ψ的初值Ψ*，然后利用带惩罚项的非线性优化采用以下公式求解得到最终的周期性四维向量场Ψ作为周期向量场：
计算中，wr是控制惩罚项的权重wr＝1，在本例中，|V|表示顶点的个数。
3)从优化得到的周期向量场中，对每一个顶点p计算得到与之相容的局部参数化坐标
利用生成树将所有局部参数化坐标组合构成全局参数化坐标；从生成树生成网格的部分将满足无缝与边界条件的区域抽取获得平凡区域的四边形网格QR，流形三角网格M中被QR覆盖的区域为平凡区域；
4)利用生成树将所有局部参数化坐标组合构成全局参数化坐标
4.1)以三角面片作为树的节点、以流形三角网格M的对偶边作为树的边构造生成树，将生成树中的根节点三角面片记为tmin，具体实施中随机选取一个三角面片为根节点；
4.2)按照以下方式从根节点出发构造全局参数化坐标，全局参数化坐标为一个二维向量：
4.2.1)根节点tmin的全局参数化坐标计算
设置根节点tmin处的第一个顶点p1的全局参数化坐标为然后对于根节点tmin处的其余两个顶点p2和p3，按照如下公式计算两个顶点的全局参数化坐标：
4.2.2)从根节点tmin开始，通过遍历生成树的节点求出全局参数化坐标。
由生成树生成网格的部分，满足无缝与特征对齐条件的区域，称为平凡区域MR。在MR中，通过抽取半整数点(坐标值为0.5的整数倍)和整数点(坐标值都是整数)得到对应的四边形网格。如果某个半整数点周围的四个整数点不能通过等值线相连接，则忽略该半整数点。对于每一个满足条件的半整数点，得到一个对应的四边形，构成平凡区域的四边形集合QR(如图5中黑色网格线的部分)。
5)5.1)取流形三角网格M中未被平凡区域MR的四边形网格QR覆盖的区域作为奇异区域MS，即MS＝M/QR；
5.2)在奇异区域，利用周期性四维向量场的第零分量构建莫尔斯-斯梅尔复形，即抽取所有极大值点、极小值点和鞍点，从每个鞍点出发搜索最近的两个极大值点和两个极小值点构成一个四边形，最终形成奇异区域的四边形网格QS(如图5中填充有黑白渐变图案的部分)。
6)通过合并平凡区域和奇异区域的网格MS∪MR，可以获得最终的四边形网格Q输出(如图6中的四边形网格)。
7)可以看到，在生成的四边形网格中，在保证纯四边形网格的前提下，对于边界的部分，生成的网格很好地保留了模型的边界，在非边界的区域，生成的网格很好地体现了用户输入的标架场要求。
上述具体实施方式用来解释说明本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明作出的任何修改和改变，都落入本发明的保护范围。