基于Laplace调和方程的网格重构算法

陈华伟　伍权　徐卫平　余泽云

摘要：针对三角形网格向四边形网格的转化问题，基于调和方程构建模型梯度场，追踪表面流线，实现了参数化网格重构。首先，建立了基于离散Laplace方程的三角网格梯度场理论模型、数据结构模型和稀疏矩阵求解方案;其次，提出了局部坐标变换和参数方程相结合求解流线节点的统一算法，并针对流线跟踪无交点、有多个交点等特殊情况，提出了梯度收敛、最短距离和参数极值等优选策略;最后，通过模型实验验证了算法。结果表明，流线网格具有等参、闭合特点，复杂模型网格划分没有歧义，而且网格质量随网格密度增加而提高。因此，相对于传统几何重构算法，数学方法对网格重构表达具有鲁棒性和唯一性，且应用场景更广泛。

关键词：算法理论;网格重构;参数化网格;调和方程;流线追踪

中图分类号：TP3016文献标志码：A

CHEN Huawei，WU Quan，XU Weiping， et al.Mesh reconstruction algorithm based on Laplace harmonic equation[J].Journal of Hebei University of Science and Technology，2019，40（3）：199-207.Mesh reconstruction algorithm based on Laplace harmonic equation

CHEN Huawei1， WU Quan1， XU Weiping1， YU Zeyun2

（1.School of Mechanical and Electrical Engineering， Guizhou Normal University， Guiyang， Guizhou 550025， China;2.Department of Computer Science， University of Wisconsin-Milwaukee， Milwaukee 53206， USA）

Abstract：Upon the problem of transforming triangular mesh into quadrilateral mesh， gradient field via harmonic equation is created， integral flow is tracked， and parameterized mesh is reconstructed. First， gradient field construction theory， data structure model and solution scheme of sparse matrix are constructed. Then， one uniform algorithm for solving flow line node by integration of local coordinate transformation and parameterized equation is advanced， and schemes like gradient convergence， shortest distance and extreme parameter are optimally occupied upon special cases of no intersection or multiple intersections when tracing flow line. Finally， the algorithm is verified via case studies. The result shows that flow lines are characterized as iso-parameterized and closed， mesh reconstruction of complicate models has no bifurcation， and mesh quality is increased with its density. The result proves that mathematical method has robustness and uniqueness for mesh reconstruction representation compared with traditional geometry method， and it will have more application scenarios.

Keywords：theory of algorithms; mesh reconstruction; parameterized mesh; harmonic equation; integral flow tracking

非結构化网格没有规则的拓扑结构，网格划分灵活性强，但是占用内存大，后续应用的数值计算量大;结构化网格有更好的网格质量，边界拟合性好，数值计算性能高、扩散性小。网格转化、重构及其衍生问题是学者持续关注的课题。以Q-Morph算法[1-2]为代表的几何推演迭代方法存在计算量大，优质网格生成还需要相当的后置处理等问题。其他的如方向场法[3]、边界简化法[4]、最少奇异点模板法[5]、体素法[6]等，可以通过参数优化生成自适应四边形网格，但都存在单元方向和单元密度难以控制等问题。在数学上，离散网格模型即为分段流形。基于Morse理论，为流形曲面选择合适的Morse函数，就能保证Morse-Smale复形将该曲面划分为纯四边形结构[7-8]。该方法很容易实现模型特征约束和单元密度控制，但是由于模型的复杂性和标量场数据的多样性，很难鲁棒地构建Morse-Smale复形。为此，越来越多的学者开始采用调和方程进行等参数化网格的重构研究[9-11]。本文以三角网格模型为研究对象，展开离散Laplace调和方程创建和求解、流线跟踪及流线网格重构方面的研究。该研究可扩展至有限元网格划分，用于材料、流体力学等综合性能分析;可扩展至图像处理领域，用于医学分层图像、地理等高图像重构和网格分析;还可扩展至仿生结构设计，用于多孔结构、微观分形结构建模。

1三角网格模型的场求解

在拓扑几何上，三角网格模型为节点V和三角面F集合构成的分段线性流形，即M=（V， F）。根据应用需求，为每个节点vi∈V赋标量值ui，即构造M上的标量场u：V。

河北科技大学学报2019年第3期陈华伟，等：基于Laplace调和方程的网格重构算法1.1理论分析

对标量场u构造Laplace调和方程[12-13]：Δu=0。    （1）给定第一类边界条件，即狄利克雷（Dirichlet）条件：u|C=φ，   （2）其中：C为边界;φ为边界函数。

对式（1）和式（2）构成的约束Laplace方程进行离散化，有：Δui=∑j∈Niωij（uj-ui）=0，          （3）

ui=ci，ci∈C 。   （4）图1三角形上的标量和矢量

Fig.1Scalar and vector on triangle

Ni表示节点i的一阶邻域，三角网格中Ni为i的一环（1-ring）邻域。三角网格离散化，一般使用cotangent余切值计算Laplace系数[14]：ωij=[cot（αij）+cot（βij）]/2 ，   （5）式（3）—式（5）构成线性方程组，通过方程组的求解，得标量场ui∈u。

进一步地，可求解其对应的梯度场（矢量场）。对于分段线性流形，在每个分段内梯度矢量为常量。如图1所示，已知三角形T0（v0，v1，v2）或（e0，e1，e2） ，对应的标量参数u0，u1，u2，三角形面的单位法矢为N，则T0的梯度G0可通过以下线性系统求解：[e0，e1，N]T[G0]=[v1-v0，v2-v1，n]T[G0]=[u1-u0，u2-u1，0]T 。     （6）观察式（6）易知，G0是边e0和e1的线性组合，即：G0=ae0+be1=[e0，e1][a，b]T 。                      （7）将式（7）代入式（6），则有：[e0，e1]T[G0]=[e0，e1]T[e0，e1][a，b]T=[v20，v0v1;v1v0，v21][a，b]T=[u1-u0，u2-u1]T 。（8）至此，系統降为二维线性系统，可直接通过手工计算求解参数a，b。令Δ=v20×v21-（v0v1）2，则有：a=[（u1-u0）×v21-（u2-u1）×v0v1]/Δ ，

b=[（u2-u1）×v20-（u1-u0）×v0v1]/Δ 。把参数a，b代入式（7）即可得G0，对应地，可求正交梯度：G⊥=n×G0 。                                    （9）1.2数据准备

面向Laplace方程的构造，三角网格模型的数据准备工作主要有2个方面，即节点数据结构的创建和节点类型划分。

1）节点数据结构的创建

节点数据结构应同时反映邻域点和邻域三角形，因此，最简单的方法是记录节点的一环邻域边，本文节点结构体NodeDetail中主要记录节点的对边：

struct NodeDetail

{

long v[2];//记录对边两点的序号

long idxF;//记录三角形的序号

long posF;//posF=0， 1， 2，记录当前点在该三角形三点中的顺序号

};

显然1个NodeDetail对象代表了节点的1个邻接三角形，如果节点i有Ni个邻接三角形，则应有Ni个NodeDetail对象。

2）节点分类

为了保证Laplace方程有非零解，还需要设置边界约束条件。对应于Laplace方程，将三角形网格中事先指定的约束点和网格边界点都归为约束节点，其他节点（均为内部节点）全部归为非约束节点，如图2所示。如果模型为闭合曲面，则认为模型无边界点，只有事先指定点为约束节点。

1.3Laplace方程的构造与求解

式（3）中，每个节点只与其临近的一阶邻域节点存在系数关系，与其他非邻接节点无系数关系，具有局部影响性。因此，面向全局网格节点构建Laplace方程，其系数矩阵是稀疏矩阵[15]，对此本文采用了SuperLU稀疏矩阵求解方案，其时间和空间复杂度分别为T（n）=O（n1.5），S（n）=O（nlogn）。

为了减少数据存贮量，需在程序设计中直接构造稀疏矩阵。建立LaplaceData数据结构，存贮与当前节点邻接节点的索引和权重系数，并使用SortLaplaceData函数按节点索引对邻接节点进行排序。

struct LaplaceData

{

long idx;//节点索引

float data;//cotangent计算值

};

static int SortLaplaceData（const void* e1， const void* e2）

{

LaplaceData* e11 = （LaplaceData*）e1;

LaplaceData* e22 = （LaplaceData*）e2;

return e11->idx - e22->idx;

}

直观上，邻域节点的LaplaceData数据按行存贮。但SuperLU采用Harwell-Boeing矩阵存贮格式[16]，即列压缩格式，需要使用dCompRow_to_CompCol函数将行存贮参数转换为列存贮参数，再使用dCreate_CompCol_Matrix函数构造列压缩存贮的稀疏矩阵。最后，使用dgssv函数解矩阵方程，求解结果就是对应于每个节点的参数ui∈u。进一步地，可求解三角网格上的梯度场。

2流线跟踪

交互输入种子点，即可通过该点沿梯度方向追踪一条唯一的流线。设种子三角形Ts（v0，v1，v2），T0内一点ptSeed为种子点，Ts内的梯度矢量为Gs，Ts中梯度线交点为ptInt0和ptInt1，对应的共边三角形为T0和T1，如图3所示。

设最小极点和最大极点分别为vmin和vmax，则有如下流线跟踪的总体过程：

1）在Ts内，过点ptSeed和方向Gs做直线，与Ts的2个交点ptInt0和ptInt1，即为初始节点;

2）在共边三角形T0内，从ptInt0出发，沿T0的负梯度方向搜索下一个交点，记为ptInt2;

3）重复步骤2），直到ptInt2=vmin;

4）在共边三角形T1内，从ptInt1出发，沿T1的梯度方向搜索下一个交点，记为ptInt2;

5）重复步骤4），直到ptInt2=vmax;

6）以上交点按vmin至vmax的正梯度方向整理输出为流线。

2.1初始梯度线的计算

该问题可描述为求过ptSeed且方向为Gs的直线L与Ts各边的2个交点ptInt0和ptInt1。为了便于叙述，假设网格内三角形方向DT为逆时针，法向为N。令直线L与Ts的三条边分别求交，最终解出2个交点。

本文使用局部坐标变换和参数方程法统一求解。首先，以三角形一个角点为坐标原点，以该点的下一条有序边（DT方向）为X轴，使用Y=N×X，即X轴沿DT方向旋转90°构造Y轴，形成三角面上的局部坐标系XOY。假设当前角点为v0，则X轴与v0v1边重合，方向一致。接下来，将原三维坐标系参数转换至XOY二維坐标系，分别构造直线L和v0v1的参数方程，并联立求解即可求得直线L和v0v1的交点，如图4所示。

v0v1直线方程为

y=0，0≤x≤v1.x。

直线L的参数方程为

x=ptSeed.x+t×Gsx，y=ptSeed.y+t×Gsy，0≤t≤1。

联立求解，可得交点参数t，记为t0：

t0=-ptSeed.y/ptSeed.x 。

相应地：

x0=ptSeed.x+t0×Gs.x 。

接下来，验证x0的取值范围，如果满足0≤x0≤v1.x，则求得一个交点（x0，0） ，否则直线L与三角形边v0v1无交点。

如果有交点，则需要如下的后置处理：

1）通过逆变换将二维点 （x0， 0） 转换至原三维坐标系，求得交点的三维坐标;

2）搜索与Ts共v0v1边的三角形T0或T1，用于梯度线的连续搜索。

依次将原点轮换至另外2个三角形角点，按上述过程继续求解直线L与三角形另外两条边v1v2，v2v0的交点。三次求解结束后，最终可得2个交点ptInt0和ptInt1。

2.2从边出发追踪交点

该问题可描述如下：已知三角形T0（v0，v1，v2），T0内的梯度矢量为G0，其中一边（假设为v0v1）上有一点ptSeed，过ptSeed且方向为G0的直线L与T0的另外两边有一个交点ptInt，求该交点，如图5所示。

1）L与v1v2求交

构造图5所示的局部坐标系XOY。显然，XOY坐标系中，ptSeed.y=0。直线L的参数方程同2.1节，直线v1v2采用斜截式方程：

x=ky+b 。

不难得出，斜率k=（v2.x-v1.x）/v2.y，X轴截距b=v1.x。以上各式联立求解，得交点坐标：

（x0，y0）=（ptSeed.x+t0×G0.x，t0×G0.y）。

接下来，判断交点有效性。构造直线v1v2的参数方程：

x=v1x+s×（v2.x-v1.x），y=s×v2.y，0≤s≤1。

将y0代入，可解参数s：

s=y0/v2.y 。

如果满足0≤s≤1，则说明交点（x0，y0）在线段v1v2之内，为有效交点，否则认为直线L与v1v2边无交点。

有交点的情况下，按前述方法做后置处理，提前结束搜索;否则，继续求解L与v2v0的交点。

2）L与v2v0求交

参数方程可表示为

x=s×v2.x，y=s×v2.y，0≤s≤1。

同理求解交点坐标，并判断参数s的有效性。

上述计算中，需要针对以下2种情况做特殊处理：

1）当s=0或1，说明交点与节点重合。此时，无法输出共边三角形，只能按2.3节的方法继续搜索。

2）如果均不满足0 ≤s≤1条件，即有向直线L与线段v1v2和v2v0均无交点，此时，将会出现梯度收敛现象。如图6所示，当前三角形T0内，过点ptSeed沿G0方向，无交点。但是，综合T0及其共边三角形T2可知，梯度矢量G0和G2是向节点v1收敛的。可见，该情况下节点v1应为所求点。

2.3从节点出发追踪交点

当搜索起点为三角形节点时，本文按以下方法处理：遍历当前节点的所有邻接三角形，在每个三角形内分别求交点，并判断交点有效性。

该问题可描述如下：已知节点v0及其一阶邻域三角形Ti，i=0，1，…，m，m为邻域三角形数，对应梯度矢量为Gi，对每个三角形Ti，求过v0且方向为Gi的直线L与v0对边（v1v2）的交点，最终在所有有效交点中输出一个最优点，如图7所示。

图7Case3：从节点出发追踪

Fig.7Case3： Track from node如2.2节所述，只需将搜索起点置于原点，即ptSeed=v0，即可同理求解。

接下来根据参数t0和s判断交点有效性，如表1所示。

序号条件交点1t>0， s ≤ 0v12t>0， s ≥ 1v23t>0， 0 < s < 1ptInt4t0≤0无交点

如果出现多个有效节点，则需要进行冲突消解处理。可根据实验效果，实施以下2种消解策略。

1）最短距离策略

对所有有效交点，计算至v0的距离，将最短距离的交点作为输出。

2）s极值策略

对所有有效交点，比较对应的参数s，s<0时，令s0=|s|，s>1时，令s0=|s-1|，将s0极小值对应的交点作为输出。

3流线网格生成

按上述方法，可在三角网格表面追踪一条过种子点ptSeed，且两端分别抵达极点vmin，vmax的流线，记为种子线lnSeed。接下来基于种子线，重建模型表面的四边形网格。设梯度方向及其正交方向网格密度为nu×nv，则有以下重建流程：

1）分段累计种子线lnSeed的长度，等分lnSeed为nu份，得到（nu-1）个等分点ptSegi，i=0，1，…，nu-2;

2）以等分点ptSegi为新的种子点ptSeedi，从该点出发，沿G⊥方向跟踪正交流线（等参线）;

3）在（nu-1）条G⊥向流线中，选择最长的一条作为G⊥向种子线lnSeed⊥;

4）等分lnSeed⊥为nv份，再从这些等分点出发沿G向追踪，得到（nv-1）条G向流线。

4实验

在Windows系统下，以VC和OpenGL为开发平台，集成SuperLU稀疏矩阵求解方案，实现了上述算法。

图8使用球体网格模型（节点数N=2 562，三角形数T=5 120）展示了算法各主要步骤的结果。图8 b）中，极点附近梯度矢量从四面八方指向极点，这保证了所有G向流线都能最终收敛至极点，如图8 f）图8球体模型（N=2 562， T=5 120）

Fig.8Sphere model（N=2 562， T=5 120）

所示;图8 d）中，G⊥向流线即等参线均为闭合曲线，从lnSeed上的每个等分点均可追踪出一条等参线，该等分点是对应等参线的起止点。

图9使用人体网格模型展示了不同种子点对流线网格的影响。图9 a）和图9 c）分别在人体正面和背面选定不同的种子点ptSeed，可分别在人体正面和背面追踪种子线;图9 b）和图9 d）为2种情况下的流线网格，对比来看流线网格并未发生变化。理论上，只要极点vmin，vmax一定，所求的梯度场G及正交梯度场G⊥就是确定的，对应地，模型表面流线和等参线也是确定的。

图10使用股骨网格模型展示了不同网格密度下的流线网格。图10 b）—d）分别是网格密度nu×nv为10×10，100×20，400×100的结果。图10 c）中外侧髁部位（图示右上部）没有G向流线，图10 d）中，在网格加密后，同一部位划分出了G向流线，对该部位及整体模型的表达也更准确。

以角度偏斜度（Skewness）为网格质量评价参数，Skewness=max[（θmax-θe）/（180-θe），（θe-θmin）/θe]，其中θmax和θmin为单元中的最大角和最小角，θe为单元均衡角，本文网格在极点处为三角形网格（图8 f）），因此θe有60（三角形）和90（四边形）2个取值。对图10的网格质量进行分析，10×10，100×20，400×100网格的平均偏斜度分别为0.68，0.33，0.25，分别达到良、优、优的网格质量等级，3种情况下，均未出现偏斜度大于090的网格。

与传统几何方法相比，如文献＼[2＼]的Q-Morph方法，调和方程法是以底层网格的梯度场为指导重新划分网格，这一点对复杂网格同样实用，而且，网格精度和质量很容易通过网格密度进行调整。这些特点也体现在实验结果之中：1）分段流形的梯度场是流线追踪的依据，相对于传统几何方法，数学方法对网格重构表达的鲁棒性更高;2）给定模型、约束条件及网格密度等條件，可确定唯一的梯度场和流线网格;3）网格密度越高，流线网格越能准确表达原始模型，但网格密度过高，将延长流线追踪时间，增加存贮空间。

此外，借助稀疏矩阵方案，以上实验模型的梯度场求解均没有时间和空间压力;而流线追踪过程主要是利用线性方程求解交点，也不会存在性能问题。

5结语

将调和理论应用至分段流形网格重构可有效解决几何重构法中存在的算法不稳定、迭代时间长等问题。论文针对三角化网格重构问题，使用离散Laplace调和方程求解模型梯度场，进而通过流线追踪重构等参数化网格，并通过算法设计、开发和模型实验验证了该过程。相对传统几何算法，调和理论算法更具鲁棒性和网格表达的唯一性，在参数化网格重构研究领域具有扩展应用的能力，具体如下。

1）双极点约束模型很好地验证了前文的理论和方法，而且每个网格内都有确定的梯度方向，因此能有效规避零曲率曲面或平面区域网格再划分问题。但是，双极点法只能表达一个方向的等参线（G⊥向的闭合流线），流线网格并非等参网格。而且，流线网格也非纯四边形网格，在极点处网格退化为三角形。通过改变约束条件可有效解决以上问题，例如，文献＼[17—18＼]使用8点约束在模型表面构造六面体模型，从而很好地重构出纯四边形的等参数化网格。

2）对多亏格（genius-n）模型＼[19＼]，本文的等参线追踪方法需要改进。如图9，人体模型两腿之间有一个孔洞，为genius-1模型，在双腿部位，只在一侧追踪了等参线。genius-n模型中，在同一高度参数下可能追踪出等参线，这一点可通过完善流线追踪方法实现。此外，也可以采用区域分割法将genius-n模型“降维”为无亏格（genius-0）模型进行处理。在分支特征（图9手臂部位）或尖锐特征（图10外侧髁部位）处，即使加密了整体网格，流线还是难以到达、局部网格加密困难等问题，此时可以采用区域分割法＼[20＼]对这些特征区域单独划分网格。

参考文献/References：

[1]汪攀， 张见明， 韩磊， 等. 基于带约束前沿推进的四边形网格生成方法[J]. 湖南大学学报（自然科学版）， 2017， 44（8）： 29-34.

WANG Pan， ZHANG Jianming， HAN Lei， et al. An advancing front quadrilateral mesh generation method with constraint[J]. Journal of Hunan University（Natural Sciences）， 2017， 44（8）： 29-34.

[2]吴禄慎， 王伟杰， 陈华伟， 等. 基于区域划分的各向异性四边形网格重建方法[J]. 机械设计与制造， 2015（3）： 212-216.

WU Lushen， WANG Weijie， CHEN Huawei， et al. Anisotropic quadrilateral remeshing based on region division[J]. Machinery Design & Manufacture， 2015（3）： 212-216.

[3]李天华. 基于方向场的四边形重网格化技术研究[D]. 长春：吉林大学， 2017.

LI Tianhua. Research on Quad Remeshing based on Vector Field [D]. Changchun： Jilin University， 2017.

[4]徐岗， 舒来新， 朱亚光， 等. 边界简化与多目标优化相结合的高质量四边形网格生成[J]. 中国图象图形学报， 2018， 23（1）： 61-73.

XU Gang， SHU Laixin， ZHU Yaguang， et al. High-quality quadrilateral mesh generation by combining boundary simplification and multi-objective optimization[J]. Journal of Image and Graphics， 2018， 23（1）：61-73.

[5]VERMA C S， SURESH K. Alpha-MST： A robust unified algorithm for quadrilateral mesh adaptation[J]. Computer-Aided Design， 2018， 103（SI）： 47-60.

[6]郑永川， 关柏良， 林淑金， 等. 针对密集点云的快速自适应四边形网格生成算法[J]. 计算机辅助设计与图形学学报， 2019， 31（1）： 39-46.

ZHENG Yongchuan， GUAN Boliang， LIN Shujin， et al. Rapid and adaptive quadrilateral mesh generation algorithm from dense point cloud [J]. Journal of Computer-Aided Design & Computer Graphics， 2019， 31（1）： 39-46.

[7]PATANE G. Mesh-based and meshless design and approximation of scalar functions[J]. Computer Aided Geometric Design， 2017， 57： 23-43.

[8]袁洁， 周明全， 耿国华， 等.  基于Morse-Smale拓扑特征的文物碎片拼接算法[J]. 自动化学报， 2018， 44（8）： 1486-1495.

YUAN Jie， ZHOU Mingquan， GENG Guohua， et al. Automatic reassembly of fractured fragments using Morse topological features[J]. ACTA Automatica Sinica， 2018， 44（8）： 1486-1495.

[9]齊学义， 杨帆. 双调和方程生成正交曲线坐标网格技术的研究[J]. 水力发电学报， 2000（3）： 84-89.

QI Xueyi， YANG Fan. Study on technique for generating boundary-orthogonal curvilinear coordinate grid using the biharmonic equations [J]. Journal of Hydroelectric Engineering， 2000（3）： 84-89.

[10]LING R， HUANG J， SUN F， et al. Spectral quadrangulation with feature curve alignment and element size control[J]. ACM Transactions on Graphics， 2014， 34（1）： 1-11.

[11]MARCHANDISE E， WIART C， VOS W G， et al. High-quality surface remeshing using harmonic maps（part II）： Surfaces with high genus and of large aspect ratio[J]. International Journal for Numerical Methods in Engineering， 2011， 86（11）： 1303-1321.

[12]ABDULHADI Z， MUHANNA Y， PONNUSAMY S. Dirichlet problem， univalency and schwarz lemma for biharmonic mappings[J]. Mediterranean Journal of Mathematics， 2018， 15（4）： 187.

[13]BEZRODNYKH S， VALSOV V. On a problem of the constructive theory of harmonic mappings[J]. Journal of Mathematical Sciences， 2014， 201（6）： 705-732.

[14]GLICKENSTEIN D. Geometric triangulations and discrete Laplacians on manifolds[J]. Mathematics， 2005：1-43.

[15]紀国良， 丁勇， 周曼， 等. 工程计算中大型稀疏矩阵存储方法研究[J]. 数值计算与计算机应用， 2018， 39（3）： 217-230.

JI Guoliang， DING Yong， ZHOU Man， et al. Research on the storage method of large-scale sparse matrix in engineering calculation[J]. Journal on Numerical Methods and Computer Applications， 2018， 39（3）： 217-230.

[16]张永杰， 孙秦. 稀疏矩阵存储技术[J]. 长春理工大学学报， 2006， 29（3）： 38-41.

ZHANG Yongjie， SUN Qin. Sparse storage technique for sparse matrix[J]. Journal of Changchun University of Science and Technology， 2006， 29（3）： 38-41.

[17]CHEN Huawei， GUO Ye， ROSTAMI R， et al. Porous structure design using parameterized hexahedral meshes and triply periodic minimal surfaces[C]// Computer Graphics International 2018. Bintan Island：[s.n.]， 2018.

[18]MARTIN T， COHEN E， KIRBY R M. Volumetric parameterization and trivariate B-spline fitting using harmonic functions[J]. Computer Aided Geometric Design， 2009（26）： 648-664.

[19]钱坤， 张家玲， 李映华， 等. 高亏格曲面共形参数化方法[J]. 计算机辅助设计与图形学学报， 2017， 29（12）： 2225-2234.

QIAN Kun， ZHANG Jialing， LI Yinghua， et al. Conformal parameterization for high genus surfaces[J]. Journal of Computer-Aided Design & Computer Graphics， 2017， 29（12）： 2225-2234.

[20]LI Bo， QIN Hong. Feature-aware reconstruction of volume data via trivariate splines [C]// Pacific Graphics 2011. Kaohsiung：[s.n.]， 2011.第40卷第3期河北科技大学学报Vol.40，No.3

2019年6月Journal of Hebei University of Science and TechnologyJune 2019