限制失真的网格参数化方法

蔡兴泉 孙辰 葛亚坤

摘 要：针对当前网格参数化效率较低、映射失真较严重的问题，提出一种限制失真的网格参数化方法。首先，预处理原始网格模型。输入原3D网格模型，采用Half-Edge数据结构来重新组织网格并切割网格模型产生相应的切缝;构建Tutte映射把3D网格映射到一个2D凸多边形域，即构建2D网格模型。然后，进行限制失真的网格参数化计算。将Tutte映射后的2D网格模型作为限制失真计算的初始数据，建立相对于原3D模型网格的失真度量函数;求得该度量函数的最小值点，即为映射后的网格坐标集合;将映射后的网格作为限制失真映射的输入网格，设定迭代终止条件，循环迭代直至迭代结束，得到收敛的最优网格坐标;在计算映射失真度时，针对等距映射失真采用Dirichlet能量函数度量，针对共形映射失真采用尽可能等距（MIPS）能量函数度量;在求解映射失真度量函数的最小值點时采用代理函数法结合组合牛顿法的最优解方法。最终，实现了该方法并开发了一个原型系统。在原型系统中，分别设计了限制等距失真和限制共形失真的网格参数化实验，对程序执行时间和失真能量下降情况进行了统计和对比，提供了相应的纹理映射效果展示。实验数据表明，所提出的方法执行效率高、映射失真能量下降快，最优值收敛质量稳定;纹理映射时纹理着色均匀、布局紧致、线条均匀，符合实际应用的标准。

关键词：网格参数化;限制失真;等距映射;共形映射;能量函数;最优坐标点

中图分类号： TP391.9

文献标志码：A

Abstract：Aiming at the low efficiency and serious mapping distortion of current mesh parameterization， a mesh parameterization method with limiting distortion was proposed. Firstly， the original mesh model was pre-processed. After inputting the original 3D mesh model， the Half-Edge data structure was used to reorganize the mesh and the corresponding seams were generated by cutting the mesh model. The Tutte mapping was constructed to map the 3D mesh to a 2D convex polygon domain， that is to construct the 2D mesh model. Then， the mesh parameterization calculation with limiting distortion was performed. The Tutte-mapped 2D mesh model was used as the initial data for limiting distortion calculation， and the distortion metric function relative to the original 3D model mesh was established. The minimum value points of the metric function were obtained， which form the mapped mesh coordinate set. The mapped mesh was used as the input mesh to limit the distortion mapping， and the iteration termination condition was set. The iteration was performed cyclically until the termination  condition was satisfied， and the convergent optimal mesh coordinates were obtained. In calculating the mapping distortion， the Dirichlet energy function was used to measure the isometric mapping distortion， and the Most Isometric Parameterizations （MIPS） energy function was used for the conformal mapping distortion. The minimum of the mapping distortion energy function was solved by proxy function combining assembly-Newton method. Finally， this method was implemented and a prototype system was developed. In the prototype system， mesh parameterization experiments for limiting isometric distortion and limiting conformal distortion were designed respectively， statistics and comparisons of program execution time and distortion energy falling were performed， and the corresponding texture mapping effects were displayed. Experimental  results show that the proposed method has high execution efficiency， fast falling speed of mapping distortion energy and stable quality of optimal value convergence. When texture mapping is performed， the texture is evenly colored， close laid and uniformly lined， which meets the practical application standards.Key words：  mesh parameterization; limiting distortion; isometric mapping; conformal mapping; energy function; optimal coordinate point

0 引言

网格参数化是数字几何处理中的热点和难点[1]。网格参数化的任务是将三维空间的网格映射到二维空间，称为平面参数化。网格参数化可用在纹理映射、网格变形、网格编辑等领域。

最初，数学家Tutte[2]证明了网格模型任何一个顶点的坐标可表示为其邻接顶点坐标的加权组合，提出了边界域凸组合映射，为网格参数化提供了理论基础。在映射计算时，容易产生失真，如三角网格的边长可能会变，三角网格的内角可能会变。三角网格的边长在映射前后发生变化称为等距失真;三角网格的内角在映射前后发生变化称为共形失真。但是，Tutte算法在映射的过程中并没有对映射失真进行约束。Hormann等[3]提出了保證无网格翻转的映射，以保证网格不发生翻转并降低映射产生的失真。Lipman[4]于2012年提出了限定失真的映射，设定了失真阈值，从而进一步减少映射带来的失真，并拓宽了网格参数化的应用范围。受该方法启发，本文主要研究限制失真的网格参数化方法。

在限制失真时，需要构建映射失真度量函数。关于等距失真映射，Hormann提出了尽可能等距（Most Isometric Parameterizations， MIPS）能量函数，Alexa等[5]提出尽可能刚体（As Rigid As Possible， ARAP）能量函数，Smith等[6]于2015年提出了对称狄利克雷（Symmetric Dirichlet， Dirichlet）能量函数。关于共形失真映射，Sorkine等[7]提出了比值能量函数，随后又在文献[8]中提出了尽可能相似（As Similar As Possible， ASAP）能量函数。综合考虑能量函数的度量精确性和易优化性，本文选用狄利克雷能量函数作为等距映射的失真度量函数，选用尽可能等距能量函数作为共形映射的失真度量函数。

映射失真函数的最优坐标解依赖于数值最优化方法。文献[9]中常用的数值优化方法有梯度下降法、牛顿法、拟牛顿法等。失真度量函数通常是非线性非凸的，有些能量函数带有约束项。因此，需要改进这些最优化方法以解决最优坐标快速求解问题。对于梯度下降法，主要是根据原函数的一阶导数确定下降方向。Kovalsky等[10]于2016年提出利用网格的Laplacian矩阵确定能量函数初始的下降方向。该方法在迭代前期失真下降较快，但是由于一阶导数缺少自由度之间的耦合信息，导致在迭代后期收敛质量和效率较差。Rabinovich等[11]于2018年提出了代理函数方法，将难以优化求解的度量函数转化为易优化的简单函数，利用局部全局（Local-Global）的思想结合代理函数的一阶导数信息确定下降方向。该方法在最优值求解迭代初期，能量下降速度很快。对于牛顿法，主要利用能量函数的Hessian矩阵确定原函数下降方向，要保持Hessian的正定。为此，Teran等[12]提出了投影牛顿法，对Hessian矩阵进行特征值分解，将非正的特征值设定为较小的正值，保持了矩阵的正定，从而保持了最优解的稳定性。该方法在实际求解过程中其计算过于复杂，时间复杂度表现较差。Shtengel等[13]提出组合牛顿法，对失真函数进行局部凸化，把非凸函数转化为凸函数，并保持Hessian正定，保证了最优解的收敛稳定性和求解的效率。对于拟牛顿法，则介于梯度下降法和牛顿法之间，主要思想是利用一阶导数信息来近似二阶导数确定下降方向。Nocedal提出了低内存占用的拟牛顿法，此方法解决与初始解的选取无关的优化问题的能力较弱。为了在最优坐标解的求解中充分利用各种优化方法的优势，当初始解离最优解距离较远时，用代理函数法求解最优坐标;在当前解离最优解距离较近时，用组合牛顿法求解最优坐标。

当前网格参数化算法效率较低、映射失真较严重，本文主要研究限制失真的网格参数化方法。首先预处理原始网格模型，然后进行限制失真的网格参数化计算。

1 预处理原始网格模型

在进行网格参数化运算前，需要预处理原始网格模型。首先，输入的是原网格模型OriginalModel。该网格模型是以顶点集和面片集为存储格式的3D网格模型。这种3D网格模型是网格建模常用的存储结构;然后，采用OpenMesh网格处理库中高效的Half-Edge数据结构来重新组织网格数据;接着，切割网格模型，产生Seam，即在模型上产生相应的切缝;最后，利用切割好的模型构建Tutte映射，把3D网格映射到一个2D凸多边形域，即2D网格模型。

1.1 切割网格模型

输入原网格模型，利用Half-Edge数据结构来重新组织网格数据后，对原3D网格模型进行切割，产生切缝。目的是将亏格大于1、有多个洞的模型转化为亏格为0、无洞的片状拓扑模型。

对原3D网格模型进行切割时，有两种方法，即Seamster网格切割法和几何图像法。Seamster网格切割法根据3D模型的高曲率点和模型中不可视区域确定切线，采用最小Steiner树算法尽量缩短切缝长度;几何图像法采用网格参数化与网格切割交替迭代，将参数化失真严重的点作为下次切割的切点。

在实现时，本文用这两种方法完成网格模型切割。这两种方法都是运用代数拓扑理论保证切线能把任意亏格的模型转化为亏格为0的开网格。

1.2 构建Tutte映射

由于Tutte映射既可以保证网格不发生翻转，又可以保证映射后的网格不发生交叉，即可保证双射。因此，本文在完成切割网格模型后，对模型构建Tutte映射，为后续参数化计算做好准备。

对模型构建Tutte映射时，首先将网格模型的边界顶点集序列组成凸多边形的边界顶点集序列。边界顶点较多，因此所围成的凸多边形边界近似一个圆。本文在实现时，将凸多边形边界顶点集序列依次平均分布在一个单位圆上。

然后，解重心坐标映射方程（1）得到此凸多边形区域所有内点的坐标。

这些内点的坐标即构建Tutte映射后的2D网格顶点的坐标。

重心坐标映射方程把内点坐标的求解运算转化为线性系统的计算，也就是求解两个线性系统Au=和Av=。在方程（1）中，向量u和向量v代表凸多边形域所有内点的（u， v）坐标，向量和向量代表凸多边形域所有边界点的（u， v）坐标。

2 限制失真的网格参数化计算

网格参数化计算最基本的要求是：保证映射计算后的2D网格模型无网格翻转，即映射前后两个阶段的2D网格的三个顶点顺序均按逆时针顺序排列。在映射计算时，容易产生失真，包括等距失真和共形失真。在网格参数化计算的过程中，严格消除等距失真或共形失真是难以做到的，但是，采用限制失真的方法在映射计算中尽量减少失真是可行的。

本文采用限制失真的网格参数化计算，具体步骤如下：

Step1：将Tutte映射计算后的2D网格作为限制失真计算的初始数据，建立相对于原3D模型网格OriginalModel的失真度量函数;

Step2：利用失真度量函数求得该度量函数的最小值点，即为映射后的网格坐标集合;

Step3：将映射后的网格作为限制失真映射的输入网格，设定迭代终止条件，循环迭代Step1～2，直至迭代结束，得到收敛的最优网格坐标。

2.1 构建失真度量函数

网格参数化的处理过程实际上是每个网格的映射过程。该映射类型为分段线性映射。该映射定义在每个三角网格上，表达式如式（3）所示：

其中：J是映射函数的Jacobian二阶矩阵，决定了映射的性质;b是位移平移量，为一个常量向量，在实际使用时一般用零向量代替;x、 y是映射前的网格顶点坐标;u和v是映射后的网格顶点坐标。无网格翻转的约束条件为J的行列式值大于0。

计算映射失真度需要建立失真刻画精确、计算简单、容易求得最小值的失真度量函数。鉴于此，采用Dirichlet能量函数度量等距映射失真，采用MIPS能量函数度量共形映射失真。

2.1.1 基于Dirichlet能量函数度量的等距映射

等距映射保持三角网格的边长在映射前后不变。由于Jacobian矩阵决定了变换的几何性质，J的特征值决定了J的性质，故能量函数主要以J的特征值为变量。Dirichlet能量函数描述如式（4）所示：

共形映射保证映射前后网格的内角不变，即保证相似性。MIPS能量函数最初用来约束等距映射失真，但是其约束目标主要是σ1=σ2，更倾向于共形映射的约束目标，对等距失真的约束效果不如Dirichlet函数，所以用它作为共形映射的失真能量函数。此函数同样具有计算量小、度量精确、鲁棒性强等特点，MIPS能量函数描述如式（6）所示：

2.2 基于能量函数求解最优坐标点

映射失真度量函數的最小值点是映射后网格的最优坐标点。考虑到数值最优化方法各有优势，采用代理函数法结合组合牛顿法的最优解方法。为了加速失真能量的下降，在前期失真能量值高的迭代步引入了参考网格，用参考网格代替原3D模型网格建立能量函数;接着利用此最优解方法求解映射后最优坐标，优化方法采用全局优化。通过数学理论分析结合实验数据统计得出实验的控制参数，参数具有普遍意义，对于不同的数据集有良好的鲁棒性。

当失真能量下降率大于0.1时采用代理函数法;当能量下降率小于0.1大于0.01时采用组合牛顿法;当能量下降率小于0.001时，就不再建立参考网格，直接用组合牛顿法迭代得出最优解。对于Dirichlet能量函数，当能量梯度范数小于1×10-6时迭代停止;对于MIPS能量函数，当能量梯度范数小于1×10-4时迭代停止，最大迭代步数为100。

2.2.1 基于梯度下降的代理函数最优坐标解法

代理函数法的主要思想是通过优化比能量函数更易优化的代理能量函数来间接优化能量函数，而代理能量的优化可以看作尽可能刚体失真能量的局部全局优化方法的推广。该方法的代理函数P（x）写成式（8）的形式：

其中：矩阵U和V是当前Jacobian矩阵做特征值分解后的两个旋转矩阵;k指当前第k次迭代;λ取1×10-4;A表示网格的面积; f表示每个三角网格;‖·‖表示求矩阵的欧氏范数;SJ是Jacobian矩阵的特征值矩阵;D（SJ）是失真能量函数;

Λ为一个常量矩阵，根据失真能量函数确定。对于Dirichlet能量函数，旋转矩阵的特征值矩阵SΛ两个值均为1;对于MIPS能量函数，SΛ两个值均为σ1 σ2。

最后利用函数的一阶导数求得代理函数的下降方向。

2.2.2 组合牛顿最优坐标解法

组合牛顿法的主要思想是把非凸的失真能量函数做局部凸化转化为凸函数，再用牛顿法得出函数的下降方向。失真能量函数可写成复合函数，如式（11）所示：

最后利用所求得的Hessian矩阵确定能量函数的下降方向。

为求解能量函数的下降步长，首先根据能量函数的下降方向确定所有网格发生翻转时的步长，并取所有步长的最小值作为最大步长;接着根据最大步长用基于Armijo准则的回溯算法求得失真函数的实际下降步长;最后结合能量函数的下降方向和下降步长得出映射后的网格坐标，即本次迭代后的网格坐标。

2.2.3 建立参考网格

递进参数化的思想主要是建立参考网格代替原模型网格。本文在实现递进参数化的计算时，具体步骤如下：

Step1：首先设定一个失真阈值，然后根据前一个迭代步结束后生成的模型网格，再结合原模型网格建立能量函数，最后根据能量函数确定映射函数Jacobian矩阵的特征值，再结合原模型网格计算得到参考网格;

Step2：首先由前一个迭代步产生的模型网格结合当前的参考网格建立能量函数，然后求解能量函数的最优坐标值得出当前迭代步的模型网格;

Step3：设定迭代终止条件，循环迭代Step1和Step2，直至迭代结束，得到收敛的最优网格坐标。

递进参数化的原理是，较高的失真能量会影响最优解的求解速度。采用建立参考网格的方法，每次优化一个较小的失真，可以提高优化的效率。特别在选用了产生高失真能量值的能量函数，或迭代初始点离最优解距离较远时使用效果更佳。但此方法会增加内存的开销，故迭代次数不宜过多，当失真能量降到一个较低的水平时，此方法没有明显优势。对Dirichlet能量约束的等距映射，首先采用指数插值法让失真能量值等于200，方程如式（15）所示：

然后对每个网格用牛顿求根法计算t值。接着对所有网格取t的最小值，以最小值作为参数求得新的Jacobian矩阵特征值，最后根据特征值求出每个参考网格的顶点。对MIPS能量约束的共形映射，直接采用线性插值法根据约束目标迭代让特征值相互接近，失真能量值小于3时迭代停止，取当前的特征值，以此求得参考网格的顶点坐标。

3 实验与结果分析

本文设计并实现了所提出的限制失真的网格参数化方法，并开发了一个验证系统。本文实验中所采用的输入模型选自法国INRIA的Unwarpped Meshes网格模型库和中国科技大学的Disk-topology Meshes网格模型库。实验中的牛、兔、手掌模型取自Unwarpped Meshes库，鱼、人、玩偶模型取自Disk-topology Meshes网格模型库。牛模型网格数量为69451个，兔模型网格数量为34504个，手掌模型网格数量为72598个，鱼模型网格数量为16428个，人模型网格数量为19012个，玩偶模型网格数量为25082个。实验硬件环境为Intel Core i7 7700 HQ，2.8GHz，四核心八线程CPU，Samsung 2400MHz 8GB内存，NVIDIA 1050 Ti 4GB显卡，Windows 10 操作系统的笔记本电脑。在实现验证系统时，采用C++作为编程语言，采用Visual Studio 2017作为开发工具。在进行相关稀疏矩阵和稀疏向量的计算时，利用Eigen计算库完成。

网格参数化计算非常注重运行效率和映射失真的下降速度。失真分为等距失真和共形失真，因此，本文设计了限制等距失真映射实验和限制共形失真映射实验。两个实验均进行了较为客观的实验数据统计和较为主观的纹理映射效果展示。实验数据统计包括程序执行时间对比、失真能量下降情况对比等。在统计程序执行时间数据时，分别对牛、兔、手掌、鱼、人、玩偶等模型进行1000次网格参数化计算，统计时间，最终求得平均值。在统计失真能量下降情况数据时，对每个模型迭代计算，统计相应的失真能量值。在纹理映射的实验中，利用OpenGL纹理映射函数绑定参数化后的网格坐标，即原模型的纹理坐标，从而生成OriginalModel模型的纹理贴图。

3.1 限制等距失真的网格参数化实验

限制等距失真的网格参数化执行时间如表1所示。

从表1可看出：在处理相同模型时，本文算法程序执行时间明显比代理函数算法和组合牛顿算法短。

等距映射失真能量下降对比情况如图1所示，取前10次迭代。通过对比可以看出，在失真能量下降速度方面，本文算法明显比代理函数法和组合牛顿法下降速度快，收敛质量优于代理函数法。

3.2 限制共形失真的网格参数化实验

限制共形失真的网格参数化执行时间如表2所示。从表2可看出：在处理相同模型时，本文算法的程序执行时间明显比组合牛顿法的短，但比代理函数法稍长。共形映射失真能量下降对比情况如图3所示，取前10次迭代。通过对比可以看出，在失真能量下降速度方面，本文算法比代理函数法和组合牛顿法下降速度快，收敛质量优于代理函数法。模型共形映射后的结果和纹理贴图效果如图4所示。可以看出各区域纹理着色均匀，不同区域的纹理图案具有相似性。

3.3 实验结果分析

对映射实验的统计数据进行分析。对于等距映射，能量函数每一项均为平方项，所以在映射失真时所产生的能量值较大，递进式参数化策略在迭代前期能充分地发挥其优势。采用代理函数法结合组合牛顿的方法使得求解效率得到提高，既能发挥基于梯度下降的代理函数方法在前期迭代的快速下降的优势，又能发挥出牛顿法的二次收敛性和收敛稳定性。对于共形映射，由于所选用的失真能量函数每一项是比值项，所以产生的失真能量值相对较小，建立参考网格需要一定的程序时间和内存空间，所以最优坐标点的求解效率比等距映射低，但是相对于不采用参考网格的方法在程序时间上和优化结果上仍有优势。值得注意的是，在迭代前期单独使用代理函数法优化效率较高，但随着迭代次数增加，当前解逐渐靠近最优解时收敛速度很慢。

从纹理映射效果上直观分析，对于等距映射，映射产生的纹理效果布局紧致，纹理条纹呈现均匀，但在局部区域出现角度失真。对于共形映射，纹理布局在局部区域呈现松弛，出现距离失真，但是不同区域的相似性得到了体现，各区域纹理着色均匀，对于角度的失真约束更加可观。在实际应用中，在对距离失真和纹理的均匀性要求严格的情况下，建议采用等距映射的参数化方法，在对角度失真和纹理的相似性要求严格的情况下，建议采用共形映射参数化。经实验分析，该算法的鲁棒性强，参数化结果对参数的调整很稳定。

4 结语

针对当前网格参数化效率较低、映射失真较严重的问题，本文提出了一种限制失真的网格参数化方法。首先，预处理原始网格模型。然后，进行限制失真的网格参数化计算。在计算映射失真度时，采用Dirichlet能量函数度量等距映射失真;采用MIPS能量函数度量共形映射失真。在求解映射失真度量函数的最小值点时，采用代理函数法结合组合牛顿法的最优解方法。最终设计并实现了该限制失真的网格参數化方法，开发了一个验证系统，分别设计了限制等距失真和限制共形失真的网格参数化实验。实验数据表明，本文方法执行效率高，映射失真能量下降快，最优值收敛质量稳定;纹理映射时，纹理着色均匀、布局紧致、线条均匀、符合实际应用的标准。该方法可用在纹理映射、网格变形、网格编辑等领域。

下一步研究工作包括將本文方法扩展，应用于二维网格模型变形、三维网格模型变形、网格编辑、形状插值、网格质量提升等工作中。

参考文献（References）

[1] 郭凤华， 张彩明， 焦文江. 网格参数化研究进展[J]. 软件学报， 2016， 27（1）： 112-135（GUO F H， ZHANG C M， JIAO W J. Research progress on mesh parameterization[J]. Journal of Software， 2016， 27（1）： 112-135）.

[2] TUTTE W T. How to draw a graph[J]. Proceedings of the London Mathematical Society， 1963， 3（1）： 743-767.

[3] HORMANN K， GREINER G. MIPS： an efficient global parameterization method[EB/OL]. [2019-01-10]. https：//pdfs.semanticscholar.org/02e4/f09c9a6d0d770d31c9289d30b7b4e9b5d974.pdf.

[4] LIPMAN Y. Bounded distortion mapping spaces for triangular meshes[J]. ACM Transactions on Graphics， 2012， 31（4）： 108.

[5] ALEXA M， COHEN-OR D， LEVIN D. As-rigid-as-possible shape interpolation[C]// Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques. New York： ACM， 2000： 157-164.

[6] SMITH J， SCHAEFER S. Bijective parameterization with free boundaries[J]. ACM Transations on Graphics， 2015， 34（4）： 70.

[7] SORKINE O， COHEN-OR D， GOLDENTHAL R. Bounded-distortion piecewise mesh parameterization[C]// Proceedings of the 2002 Conference on Visualization. Piscataway： IEEE， 2002： 355-362.

[8] SORKINE O， ALEXA M. As-rigid-as-possible surface modeling[C]// Proceedings of the 5th Eurographics Symposium on Geometry Processing. Aire-la-Ville， Switzerland： Eurographics Association， 2007： 109-116.

[9] NOCEDAL J， WRINGT S J. Numerical Optimization[M]. Berlin： Springer， 1999： 1-10.

[10] KOVALSKY S Z， GALUN M， LIPMAN Y. Accelerated quadratic proxy for geometric optimization[J]. ACM Transactions on Graphics， 2016， 35（4）： 134.

[11] RABINOVICH M， PORANNE R， PANOZZO D， et al. Scalable locally injective maps[J]. ACM Transactions on Graphics， 2017， 36（4）： 37a.

[12] TERAN J， SIFAKIS E， IRVING G， et al. Robust quasi static finite elements and mesh simulation[C]// Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation. New York： ACM， 2005： 181-190.

[13] SHETENGEL A， PORANNE R， SORKINE-HORNUNG O， et al. Geometric optimization via composite majorization[J]. ACM Transactions on Graphics， 2017， 36（4）： 38.

[14] SHEFFER A， HART J C. Seamster： inconspicuous low-distortion texture seam layout[C]// Proceedings of the 2002 Conference on Visualization. Piscataway： IEEE， 2002： 291-298.

[15] GU X F， STEVEN J G， HOPPE H. Geometry images[J]. ACM Transactions on Graphics， 2002， 21（3）： 355-361.

[16] FLOATER M S. One-to-one piecewise linear mappings over triangulations[J]. Mathematics of Computation， 2003， 72（242）： 685-696.

[17] LIU L， YE C， NI R， et al. Progressive parameterizations[J]. ACM Transations on Graphics， 2018， 37（4）： 41.

[18] AIGERMAN N， LIPMAN Y. Injective and bounded distortion mappings in 3D[J]. ACM Transactions on Graphics， 2013， 32（4）： 106.

[19] FU X， LIU Y， GUO B. Computing locally injective mappings by advanced MIPS[J]. ACM Transactions on Graphics， 2015， 34（4）： 71.

[20] WEBER O， ZORIN D. Locally injective parameterization with arbitrary fixed boundaries[J]. ACM Transactions on Graphics， 2014， 33（4）： 75.

[21] ZHU Y， BRIDSON R， KAUFMAN D M. Blended cured quasi-newton for distortion optimization[J]. ACM Transactions on Graphics， 2018， 37（4）： 40.

[22] RODOL E， LHNER Z， BRONSTEIN A M， et al. Functional maps representation on product manifolds[J]. Comouter Graphics Forum， 2019， 38（1）： 678-689.