0引言计算无翻转和低扭曲的体映射是计算机图形与几何处理中的重要问题。无翻转体映射在网格变形、重新网格化、网格优化以及形状分析等诸多领域有着广泛应用。研究鲁棒的无翻转体映射生成方法具有重要的现实意义。现实世界中的物体包含内部结构。体映射不仅需要考虑边界，还需要处理内部复杂的几何和拓扑结构。由于现实中的物体不存在负的体积，所以要求体映射必须是无翻转的，即体映射的雅可比矩阵行列式处处大于0。生成无翻转的体映射主要有以下困难：1)输入映射可能存在翻转。由于无翻转约束是非线性和非凸的，所以将这些翻转去除是困难的。在保证体映射是无翻转的同时，还需要尽可能降低它的扭曲。2)体映射需要满足用户指定的位置约束。位置约束和无翻转约束是相互制约的，同时满足这两个约束比较困难。为了满足上述要求，已经提出了很多计算无翻转体映射的方法。体映射的计算通常视为一个带约束的优化问题来求解，一般要求在无翻转的条件下尽量减少映射的扭曲。一类方法是从一个无翻转的初始映射开始优化扭曲，通过使用线搜索的方式使得无翻转约束一直不被违反(Aigerman和Lipman，2013；Kovalsky等，2016；Zhu等，2018)。然而对于3维网格，给定严格的位置约束时，计算无翻转的初始体映射是非常困难的。另一类方法不需要无翻转的初始映射，尝试通过消除翻转的方式生成无翻转的体映射。这类方法主要分为扭曲有界的映射方法(Aigerman和Lipman，2013；Kovalsky等，2014, 2015)、基于表示的方法(Paillé等，2015；Fu等，2015)、基于惩罚的方法(Liu等，2016)等。这类方法有一些限制因素，一般需要额外的输入信息或复杂的计算。扭曲有界的映射方法需要输入精确的扭曲界线。基于惩罚的方法需要输入额外的旋转场，若输入的数据不准确，则会造成生成的体映射扭曲过高，甚至无法消除翻转。基于表示的方法需要计算二面角和仿射变换等数据，计算成本大幅增加。而基于映射的方法(practical foldover-free volumetric mapping construction，FF)(Su等，2019)能够在严格位置约束下生成无翻转的体映射，无需额外的输入信息，使用扭曲有界的投影技术即可消除翻转，但不能保证完全消除翻转。事实上，经过大量实验发现，对一些复杂模型，使用基于投影的方法处理后，仍然会有部分翻转存在。本文提出了一种新的计算体映射的方法，能够满足上面提出的几点要求。图 1是对本文方法的基本介绍。对有翻转的初始映射(图 1(b))，基于映射的方法(FF)不能完全消除翻转(图 1(c))，而本文方法可以完全消除翻转(图 1(d))。本文算法的核心是一种雅可比指导的变形算法。基于投影的方法虽然不能保证完全消除翻转，但是可以消除大部分翻转，并可以保证映射严格满足位置约束。因此本文将基于映射的方法得到的雅可比矩阵为引导，对原始网格进行变形操作，变形过程中通过线搜索保证体映射不发生翻转。此外，添加了位置能量，通过不断增大位置能量的权重保证变形结果最终能够满足位置约束。最后，在变形优化收敛后，再固定边界优化内部网格的扭曲，最终映射的扭曲达到较低水平。 图1 本文方法简介 Introduction to the proposed methodFig 1((a) input mesh; (b) initial mapping; (c) FF; (d) ours) 本文提出的雅可比引导的无翻转体映射生成方法的优点和贡献主要如下：1)不需要任何额外的输入，能够处理现有的无翻转体映射生成算法失败的复杂模型，算法鲁棒性高。2)对初始有翻转的体映射，可以生成无翻转的低扭曲体映射。3)提出一种新的雅可比引导的变形方法，能够满足用户指定的位置约束。本文测试了多种复杂模型，结果表明本文方法均能生成无翻转的体映射。如图 2中高亏格的复杂模型以及输入与输出差别较大的模型，本文方法都有较好的效果。此外，在保证满足位置约束的基础上，降低了体映射的扭曲。与现有的方法相比，本文算法具有更强的鲁棒性。 图2 复杂模型与差异较大模型 Complex model and model with big differencesFig 2((a) high-genus model; (b) human model) 图 3展示了本文算法的详细流程，其中红色四面体表示翻转的四面体。图 3(a)—(c)分别表示输入的初始网格、输入模型初步映射到多立方体上和基于映射方法的结果，图 3(d)—(f)表示本文算法的主要流程，分别为雅可比引导的变形优化、位置约束以及降低扭曲。 图3 本文算法的整体流程 The pipeline of the algorithmFig 3((a) input mesh; (b) initial mapping; (c) FF; (d) deformation optimization; (e) position constraints; (f) post-optimization processing) 1相关工作1.1基于维护的方法基于维护的方法通常从一个无翻转的初始出发，通过优化目标能量方程来降低映射的扭曲(Aigerman和Lipman，2013；Jiang等，2017)。这类方法在优化过程中通过检测雅可比矩阵的行列式来判断映射是否发生翻转，当映射的雅可比矩阵的行列式趋于0时，映射的扭曲能量将趋于正无穷。并且这类方法可以通过线搜索来控制下降步长，使得在优化能量过程中始终不产生翻转。对表面映射而言，最早的双射方法是Tutte嵌入方法(Tutte，1963)，这是唯一有理论保证的可以保证无翻转的映射方法。因此基于维护的方法通常应用于与表面映射关系密切的参数化中(Claici等，2017；Liu等，2018)。然而对于体积映射来说，想要得到一个无翻转的初始映射是困难的。例如，可以将一个四面体网格映射到一个立方体或球上，同时保证该映射是一个双射(Campen等，2016)。但是如果将原始模型映射到任意的边界形状上，就很难保证该映射不发生翻转。因此基于维护的方法很难拓展到体映射的计算中。类似基于维护的方法，本文方法同样使用能量优化来防止翻转并降低扭曲，通过变形生成无翻转的体映射。但是本文方法可以处理具有翻转的初始映射，对3维中具有复杂边界的翻转模型仍然能产生较好的结果。基于维护的方法优化的能量通常是非线性的，很难直接求解这些能量方程。通常做法是通过迭代产生一系列的近似值，直到收敛至最优值。为此，大部分方法使用目标函数的局部二次逼近或代理来求解该优化问题，根据代理的不同，基于维护的方法可以大致分为3种类型：1)一阶方法(Fu等，2016；Rabinovich等，2017)，构造比较简单，只考虑能量的一阶导数信息，可以达到不错结果，但往往收敛较慢。2)拟牛顿法(Zhu等，2018)，仅使用梯度差来模拟二阶导数，收敛速度高于一阶方法。3)二阶方法，收敛速率最高，其中牛顿法是最传统的二阶方法。该方法一般适用于凸函数，对非凸函数通常需要修正目标函数的海森矩阵(Shtengel等，2017；Golla等，2018)，使其成为半正定矩阵才能保证能量优化最终能够收敛。本文方法同样使用二阶方法(Smith等，2019)，通过一些不变量系统推导出等距能量的特征系统，能够将能量方程的海森矩阵解析地投影到半正定矩阵，可以加快优化算法的收敛速度，提高算法效率。1.2基于翻转消除的方法很多方法尝试通过消除翻转来生成无翻转的体映射。这类方法对初始映射没有要求，一般可以处理初始较差的模型，鲁棒性较好。在连续情况下，准保形映射方法(Weber等，2012)能够保证局部单射，并且符合给定的位置约束。但是该方法的单射保证并不能很好地转换为离散逼近，并且在凹角处经常容易发生翻转(Weber和Zorin，2014)。扭曲有界的映射方法(Aigerman和Lipman，2013；Kovalsky等，2014, 2015)将映射结果投影到扭曲有界的空间范围内来保证无翻转，但需要用户输入额外的界线信息。然而如何定义扭曲的界线仍然是一个困难问题，不准确的输入会导致该方法不能完全消除翻转。另外，基于惩罚的方法(Liu等，2016)通过添加一些凸函数和一些二次旋转弹性能量来惩罚翻转元素的存在，从而消除映射中的翻转。然而这种方法仍需要额外的旋转场作为输入信息。还有一些基于表示的方法(Paillé等，2015；Fu等，2015)不再将顶点作为变量，而是选择二面角或者仿射变换等数据来表示映射。这种方法虽然提高了表示的准确性，但大幅增加了计算的复杂度，体映射生成效率比较慢。基于面积的方法(Xu等，2011)通过优化无符号面积来保证映射无翻转，但这类方法具有不能排除退化三角形、导数不连续和梯度消失的缺点。Du等人(2020)将网格提升到更高的维度进行测量解决了上述问题，但也不能保证得到全局的最小值，即不能保证映射不会产生翻转。基于映射的方法(Su等，2019)使用的是有界扭曲的投影技术，不需要额外输入信息，并且具有鲁棒高效的优点，但对一些复杂模型仍然不能完全消除翻转。2雅可比矩阵引导的无翻转体映射生成算法2.1概述 2.1.1输入首先输入一个3维空间中的四面体网格，记为$\mathit{\boldsymbol{M}}$。网格的顶点数和四面体数分别为$N_v$和$N_t$。网格顶点集合$\mathit{\boldsymbol{V}} = \left\{ {{v_i}\mid i = 1, \cdots, {N_v}} \right\}$，网格四面体集合$\mathit{\boldsymbol{T}} = \left\{ {{t_i}\mid i = 1, \cdots, {N_t}} \right\}$。然后输入一个定义在四面体网格$\mathit{\boldsymbol{M}}$上的初始体映射，记为$f$。对于输入网格的$\mathit{\boldsymbol{V}}$和$\mathit{\boldsymbol{T}}$，将通过$f$映射后的结果分别记为$\mathit{\boldsymbol{\tilde V}}$和$\mathit{\boldsymbol{\tilde T}}$。即$\mathit{\boldsymbol{\tilde {V}}} = \left\{ {{{\tilde {\mathit{\boldsymbol{v}}}}_i}\mid i = } \right.\left. {1, \cdots, {N_v}} \right\}, \mathit{\boldsymbol{\tilde{ T}}} = \left\{ {{{\tilde{ \mathit{\boldsymbol{t}}}}_i}\mid i = 1, \cdots, {N_t}} \right\}$。在离散网格上定义的体映射$f$可以由定义的在每个四面体$t_i $上的分片线性映射$f_i$组成，即${f_i}(\mathit{\boldsymbol{x}}) = {\mathit{\boldsymbol{J}}_i}\mathit{\boldsymbol{x}} + {\mathit{\boldsymbol{b}}_i}$。其中$\mathit{\boldsymbol{J}}_i$表示$f_i$的雅可比矩阵，它是一个3×3的仿射矩阵，通常用来检验映射是否翻转或者用来定义映射的扭曲。对于四面体${t_i} = \Delta {v_{i1}}{v_{i2}}{v_{i3}}{v_{i4}}, {v_i} \in \mathit{\boldsymbol{V}}$，它的像为$f\left({{t_i}} \right) = \Delta {\tilde v_{i1}}{\tilde v_{i2}}{\tilde v_{ij}}{\tilde v_{i4}}, {\bar v_i} \in \mathit{\boldsymbol{\tilde V}}$。基于上述定义，可以得到雅可比矩阵${\mathit{\boldsymbol{J}}_i}$的表示形式，具体为1$\begin{array}{l}{\mathit{\boldsymbol{J}}_i} = \left[ {{{\mathit{\boldsymbol{\widetilde v}}}_{i2}} - {{\mathit{\boldsymbol{\widetilde v}}}_{i1}}, \cdots, {{\mathit{\boldsymbol{\widetilde v}}}_{i4}} - {{\mathit{\boldsymbol{\widetilde v}}}_{i1}}} \right] \cdot \\\;\;\;\;{\left[ {{\mathit{\boldsymbol{v}}_{i2}} - {\mathit{\boldsymbol{v}}_{i1}}, \cdots, {\mathit{\boldsymbol{v}}_{i4}} - {\mathit{\boldsymbol{v}}_{i1}}} \right]^{ - 1}}\end{array}$ 式中，${{\mathit{\boldsymbol{v}}_{i}}}$表示顶点$v_i$相应的位置坐标。最后输入用户给定的位置约束，即$\mathit{\boldsymbol{M}}$中的某些顶点需要映射到指定的位置。记位置约束中的顶点集合为${\mathit{\boldsymbol{V}}_p} = \left\{ {v_i^p|i = 1, \cdots, {N_p}} \right\}$，其中$N_p$表示位置约束顶点的数目。2.1.2目标及方法介绍若雅可比矩阵的行列式${\rm{det}}\;{\mathit{\boldsymbol{J}}_i} 0$，则称四面体$t_i$是无翻转的。仅当所有的四面体都不翻转时，一个映射才称为是无翻转的。本文的目标是得到一个无翻转的映射，并且满足给定的位置约束。为了完成给定的目标，现有的无翻转映射生成方法一般都将其作为一个优化问题来求解。因为输入的初始映射往往是翻转的，这些方法都是在满足位置约束的前提下尝试去除所有的翻转，但是往往没有理论保证能成功。因此，本文提出一种雅可比矩阵引导的变形算法用于计算无翻转的体映射。本文算法放松了位置约束，不依赖无翻转的初始映射。具体来讲，映射首先设置为恒等映射，然后进行两步优化: 1)使用雅可比矩阵指导映射优化，使得映射逐渐满足给定的位置约束，并通过线搜索的方式防止在优化过程中产生翻转；2)通过优化位置能量使得最终结果满足位置约束。2.1.3挑战及解决方案使用上述算法生成体映射时，主要面临以下两个挑战：1)给定的初始映射是有翻转的，因此找到合适的雅可比矩阵作为指导比较困难。2)在保证体映射没有翻转的同时，满足用户给出的位置约束且尽量降低映射的扭曲是非平凡的。针对上述挑战，分别提出以下解决方案：1)对于输入的有翻转的初始映射，首先使用现有的消除翻转的方法进行处理。但由于算法的限制，处理后的映射可能仍然存在部分翻转。本文直接使用处理后映射无翻转部分的雅可比矩阵作为本文变形算法的引导，对翻转部分使用单位矩阵作为引导。这样就得到了无翻转的引导雅可比矩阵。2)经过上述雅可比矩阵引导优化后，位置约束已经接近满足。然后添加一项位置约束能量，通过不断增大位置约束能量的权重，使得网格最终满足位置约束要求。最后固定位置约束继续优化网格的扭曲，得到一个无翻转的、满足位置约束的且扭曲较低的体映射。2.2雅可比矩阵引导的变形优化 2.2.1引导雅可比矩阵生成输入的初始映射$f$存在大量翻转，不能直接使用上述雅可比矩阵作为本文变形算法的引导。对此，采用如下方法找到合适的引导雅可比矩阵：对初始映射存在大量的翻转，首先使用基于映射的方法(Su等，2019)消除初始映射中的翻转。基于映射的方法使用扭曲有界的投影技术消除翻转，不仅可以严格地保证位置约束，并且鲁棒性较好。但与其他消除翻转方法相同的是，基于映射的方法也无法保证能够完全消除翻转。所以还需对上述映射的雅可比矩阵${\mathit{\boldsymbol{J}}_i}$进行处理，具体为2${\mathit{\boldsymbol{J}}_i} = \left\{ \begin{array}{l}{\mathit{\boldsymbol{J}}_i}\;\;\;\det {\mathit{\boldsymbol{J}}_i} 0\\\mathit{\boldsymbol{I}}\;\;\;\;其他\end{array} \right.$ 式中，$\mathit{\boldsymbol{I}} $表示单位矩阵。使用单位矩阵的原因是可以减少变形时翻转部分的扭曲。式(2)中矩阵对应的映射不存在翻转，且具备部分基于映射方法的良好性质，因此将式(2)中的矩阵作为变形算法的引导雅可比矩阵。2.2.2变形能量优化本文变形的目标是将映射变形到满足给定的位置约束。将算法中映射的雅可比矩阵记为$\mathit{\boldsymbol{J}}'$，四面体$t_i$上对应的雅可比矩阵记为$\mathit{\boldsymbol{J}}_i^'$。为了满足位置约束的目标，需要将$\mathit{\boldsymbol{J}}{'_i}$优化到接近引导雅可比矩阵${\mathit{\boldsymbol{J}}_i}$。因为$\mathit{\boldsymbol{J}}_i^\prime (\mathit{\boldsymbol{x}})$是网格顶点坐标$\mathit{\boldsymbol{x}}$的线性组合，所以本文变形算法可以使用的能量为3$\mathop {\min }\limits_\mathit{\boldsymbol{x}} E(\mathit{\boldsymbol{x}}) = \sum\limits_{{t_i} \in \mathit{\boldsymbol{T}}} {{A_i}} {\cal D}\left({\mathit{\boldsymbol{J}}_i^\prime (\mathit{\boldsymbol{x}})} \right)$ 式中，${\cal D}(\cdot)$用来衡量两个雅可比矩阵之间的差异，$A_i $表示四面体$t_i$的体积。为了定义${\cal D}(\cdot)$，首先定义矩阵${\mathit{\boldsymbol{\hat J}}_i}$，具体为4${\mathit{\boldsymbol{\hat J}}_i} \buildrel \Delta \over = \mathit{\boldsymbol{J}}_i^\prime \mathit{\boldsymbol{J}}_i^{ - 1}$ 当${\mathit{\boldsymbol{\hat J}}_i}$是一个单位矩阵时，说明此时$\mathit{\boldsymbol{J}}_i^\prime $和${\mathit{\boldsymbol{J}}_i}$完全一致，即变形结果与想要的最终结果是旋转不变的。为了保证在优化过程中不出现翻转，使用如下能量定义${\cal D}(\cdot)$，即5$\begin{array}{l}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\cal D}\left({\mathit{\boldsymbol{J}}_i^\prime (\mathit{\boldsymbol{x}})} \right) = \\\left\{ \begin{array}{l}\left({{{\cal D}_{{\rm{sym}}}}\left({{{\mathit{\boldsymbol{\hat J}}}_\mathit{i}}(\mathit{\boldsymbol{x}})} \right) + \omega {{\cal D}_{12}}\left({{{\mathit{\boldsymbol{\hat J}}}_\mathit{i}}(\mathit{\boldsymbol{x}})} \right)\;\;\det {{\mathit{\boldsymbol{\hat J}}}_\mathit{i}}(\mathit{\boldsymbol{x}}) \ge 0} \right.\\\infty \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\det {{\mathit{\boldsymbol{\hat J}}}_\mathit{i}}(\mathit{\boldsymbol{x}}) 0\end{array} \right.\end{array}$ 式中，$\omega $表示一个正的权值，${{{\cal D}_{{\rm{sym}}}}\left(\cdot \right)}$是对称狄利克雷度量(Schreiner等，2004；Smith和Schaefer，2015)，表示为6${D_{{\rm{sym}}}}\left({{{\mathit{\boldsymbol{\hat J}}}_\mathit{i}}(\mathit{\boldsymbol{x}})} \right) = \left\| {{{\mathit{\boldsymbol{\hat J}}}_\mathit{i}}(\mathit{\boldsymbol{x}})} \right\|_{\rm{F}}^2 + \left\| {\mathit{\boldsymbol{\hat J}}_i^{ - 1}(\mathit{\boldsymbol{x}})} \right\|_{\rm{F}}^2$ 式(5)中，${{\cal D}_{{\rm{L}}2}}(\cdot)$表示两个雅可比矩阵之间的${\rm{F}}$范数，具体为7${D_{{\rm{L2}}}}\left({{{\mathit{\boldsymbol{\hat {J}}}}_\mathit{i}}(\mathit{\boldsymbol{x}})} \right) = \left\| {\mathit{\boldsymbol{J}}_i^{'}(\mathit{\boldsymbol{x}}) - {\mathit{\boldsymbol{J}}_i}(\mathit{\boldsymbol{x}})} \right\|_{\rm{F}}^2$ 当${{{\mathit{\boldsymbol{\hat J}}}_\mathit{i}}}$退化时，${{\cal D}_{{\rm{sym}}}}\left(\cdot \right)$能量会变成无限大，从而可以防止在优化过程中发生翻转。在实验中发现，当$\omega $取值很小时，本文算法仍能得到较好结果。本文取$\omega $的值为${10^{ - 4}}$。2.2.3求解器本文使用牛顿法优化式(3)定义的能量。牛顿法是二阶优化方法，相较于其他一阶的优化方法，牛顿法迭代次数更少，收敛速率更快。但由于本文使用的狄利克雷能量是非凸的，所以它的海森矩阵并不能保证是半正定的，此时传统的牛顿法往往会发生停滞甚至发散。为了解决这一问题，使用解析的特征系统方法(Smith等，2019)来修正对称狄利克雷能量的海森矩阵，从而满足牛顿法中矩阵半正定性要求，使得本文的变形算法能够快速收敛。2.2.4线搜索方法在变形过程中，为了防止产生翻转的四面体, 使用线搜索方法(Smith和Schaefer，2015)确定下降步长$t$，从而保证每一步迭代都不会产生翻转。考虑每个未翻转的四面体$t_i$，记四面体上的4个顶点坐标分别为${\mathit{\boldsymbol{v}}_1}、{\mathit{\boldsymbol{v}}_2}、{\mathit{\boldsymbol{v}}_3}$和${\mathit{\boldsymbol{v}}_4}$，它们对应的搜索方向分别是${\mathit{\boldsymbol{p}}_1}、{\mathit{\boldsymbol{p}}_2}、{\mathit{\boldsymbol{p}}_3}$和${\mathit{\boldsymbol{p}}_4}$。当四面体发生退化时，四面体的有向体积会变为0。因此，可以通过求解以下方程求得该四面体的下降步长$k_i$8$\det \left(\begin{array}{l}\left({{\mathit{\boldsymbol{v}}_2} + {\mathit{\boldsymbol{p}}_2}k} \right) - \left({{\mathit{\boldsymbol{v}}_1} + {\mathit{\boldsymbol{p}}_1}k} \right)\\\left({{\mathit{\boldsymbol{v}}_3} + {\mathit{\boldsymbol{p}}_3}k} \right) - \left({{\mathit{\boldsymbol{v}}_1} + {\mathit{\boldsymbol{p}}_1}k} \right)\\\left({{\mathit{\boldsymbol{v}}_4} + {\mathit{\boldsymbol{p}}_4}k} \right) - \left({{\mathit{\boldsymbol{v}}_1} + {\mathit{\boldsymbol{p}}_1}k} \right)\end{array} \right) = 0$ 式(8)是关于下降步长$k$的三次方程，求解比较简单。本文只使用正下降方向上的线搜索，当下降步长$k_i$接近所求的最小正根时，相应的四面体则会接近退化。因此所求的最小的正根即是该四面体的最大下降步长。对所有的四面体计算式(8)中$k$的最小正根，取其中的最小值作为最大下降步长，记为$k_{\rm{max}}$。当下降步长$k \in \left({0, {\mathit{k}_{{\rm{max}}}}} \right)$时，变形后的网格不会发生翻转。2.2.5终止条件与其他消除翻转方法不同，本文算法始终不会产生翻转的四面体，所以不能以是否存在翻转作为终止条件。为了尽量满足位置约束，便于后面算法的进行，要求变形算法的能量尽可能地优化到收敛。所以终止条件设定为9$\frac{{\left| {{E^n} - {E^{n - 1}}} \right|}}{{{E^{n - 1}}}} {10^{{\rm{ - }}4}}$ 式中，${{E^n}}$表示定义的能量在第$n$次迭代的结果。当扭曲能量满足上述条件时终止变形算法。对于复杂的输入，算法可能迭代多次仍然不能达到上述收敛条件，因此还设定了一个迭代次数上限作为终止条件。取迭代上限为250次迭代，当超过迭代次数上限时也停止变形算法，通过后续的处理最后仍然可以得到较好结果。由于本文给出的引导雅可比矩阵是修正过的，即使上述能量优化到最优值，即${\mathit{\boldsymbol{\hat J}}_i}$是一个单位矩阵，变形的结果可能仍然不能满足给出的位置约束。图 4展示了实验得到的变形结果，通过图 4(b)初始映射和图 4(d)变形结果进行对比，发现变形结果与位置约束相差不大。 图4 雅可比矩阵引导的变形优化 Deformation optimization guided by Jacobian matrixFig 4((a) input mesh; (b) initial mapping; (c) FF; (d) deformation result) 2.3位置约束虽然本文给出的修正雅可比引导是不准确的，但由于基于映射方法处理后的四面体网格翻转较少，所以在变形优化收敛后，映射基本满足给定的位置约束。为了尽量地满足位置约束，将位置能量加入到定义的变形能量当中。由于此时的变形网格接近满足位置约束，因此优化起来仍然比较简单。添加位置约束后的能量表达式为10$\begin{array}{l}\mathop {\min }\limits_x {E_t}(\mathit{\boldsymbol{x}}) = {\lambda _D}\sum\limits_{{t_i} \in \mathit{\boldsymbol{T}}} {{A_i}} {\cal D}\left({\mathit{\boldsymbol{J}}_i^\prime (\mathit{\boldsymbol{x}})} \right) + \\\;\;\;\;\;\;\;\;\;\;\;\;{\lambda _S}\sum\limits_i^{{N_p}} {{{\left\| {{\mathit{\boldsymbol{v}}_i} - {{\mathit{\boldsymbol{\widetilde v}}}_i}} \right\|}^2}} \end{array}$ 式中，${\lambda _D}$和${\lambda _S}$是给定的两个参数，用来控制这两部分能量的大小。本文设定${\lambda _D}$为1。由于变形刚收敛的时候，映射与指定的位置约束相差较大，因此将${\lambda _S}$的初值设定得比较小，本文将其初始值设置为0.001。继续使用现有的优化方法优化该能量，为了保证能够满足位置约束，在能量收敛后增大${\lambda _S}$的值继续优化，直到第2项位置约束的能量值达到一个较低的水平。改变参数${\lambda _S}$的方法为11$\lambda _S^{{\rm{new}}} = \alpha {\lambda _S}$ 式中，$\lambda _S^{{\rm{new}}}$是新一轮优化使用的位置约束能量的参数，即令$ {\lambda _S} = \lambda _S^{{\rm{new}}}$。而$\alpha $是一个增长因子，本文将其设定为10 000。通过增大参数${\lambda _S}$，变形网格逐渐满足给出的位置约束。当位置能量小于给定的阈值$\theta $时，优化终止。实验所取的阈值$\theta $大小为${10^{ - 6}}$。图 5展示了位置约束优化的流程，其中图 5(e)表示给出的位置约束，图 5所示模型迭代3次就能满足位置约束要求。 图5 位置约束优化过程 Position constraint optimization processFig 5 ((a) input mesh; (b) deformation result; (c) the first iteration; (d) the third iteration; (e) position constraints) 2.4后优化处理本文算法还做了一步后处理操作，用来降低最终映射的扭曲。在带有位置约束的能量优化结束后，此时映射已经满足位置约束的要求。为了进一步降低映射的扭曲，固定位置约束，通过优化以下能量降低扭曲。12$\mathop {\min }\limits_\mathit{\boldsymbol{x}} {E_d}(\mathit{\boldsymbol{x}}) = \sum\limits_{{t_i} \in \mathit{\boldsymbol{T}}} {{A_i}} {\cal D}'\left({{\mathit{\boldsymbol{J}}_i}(\mathit{\boldsymbol{x}})} \right)$ 式中13$\begin{array}{l}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\cal D}'\left({{\mathit{\boldsymbol{J}}_i}(\mathit{\boldsymbol{x}})} \right) = \\\left\{ {\begin{array}{*{20}{l}}{\left\| {{\mathit{\boldsymbol{J}}_i}(\mathit{\boldsymbol{x}})} \right\|_{\rm{F}}^2 + \left\| {\mathit{\boldsymbol{J}}_i^{ - 1}(\mathit{\boldsymbol{x}})} \right\|_{\rm{F}}^2}&{\det {\mathit{\boldsymbol{J}}_i}(\mathit{\boldsymbol{x}}) \ge 0}\\\infty &{\det {\mathit{\boldsymbol{J}}_i}(\mathit{\boldsymbol{x}}) 0}\end{array}} \right.\end{array}$ 在优化处理中，仍然使用变形算法中的能量求解方法。因为降低扭曲时固定了位置约束，所以能够在满足位置约束的条件下降低映射的扭曲，并且在优化过程中不会产生新的翻转四面体。如图 6所示，红色代表扭曲，红色越深扭曲越大。经过后处理优化，本文算法明显降低了映射的扭曲。 图6 优化处理中的网格扭曲能量 Mesh distortion energy in optimization processFig 6((a) input mesh; (b) before optimization process; (c) after optimization process) 3实验结果与分析本文提供了一种能够鲁棒地生成无翻转和低扭曲体映射的算法。针对复杂模型进行大量实验，并通过与现有方法对比，验证本文算法的鲁棒性。本文算法基于C+ +实现，所有实验都在Core i7-8700处理器、16 GB内存的台式电脑上完成。优化能量时使用解析的特征系统方法来修正海森矩阵，并使用英特尔的数学内核库求解大型的线性方程组(Schenk等，2007；Petra等，2014a, b)。3.1质量测量标准为了评估体映射的质量，使用如下测量标准：1) 是否存在翻转。无翻转是体映射的最基本要求，通过记录四面体网格中翻转的四面体数量来衡量体映射的质量。本文所有图例中红色的四面体表示翻转的四面体。2) 体积扭曲度量。对每个雅可比矩阵${{\mathit{\boldsymbol{J}}_i}}$，将扭曲度量定义为$\delta \left({{\mathit{\boldsymbol{J}}_i}} \right) = \left\| {{\mathit{\boldsymbol{J}}_i}} \right\|_{\rm{F}}^2 + \left\| {\mathit{\boldsymbol{J}}_i^{ - 1}} \right\|_{\rm{F}}^2$。实验中首先得到网格中每个四面体的扭曲度量，然后将它们各自的体积占总体积的比例作为权重，计算网格总体的扭曲，记为${\delta _t}$。3) 位置约束度量。使用变形网格与位置约束顶点的距离来衡量位置约束满足程度。对每个约束顶点定义位置度量为$\mathit{\rho }({\mathit{\boldsymbol{v}}_i}, {\mathit{\boldsymbol{\tilde v}}_i}) = {\left\| {{\mathit{\boldsymbol{v}}_i} - {{\mathit{\boldsymbol{\tilde v}}}_i}} \right\|^2}$。将网格整体的位置约束能量记为$\mathit{\rho }$，并给定一个很小的阈值$\theta $作为是否满足位置约束的判断标准。3.2结果分析 3.2.1不同网格质量本文测试了不同质量的网格输入。如图 7所示，输入两种内部质量不同、形状相同的原始网格，位置约束是它们对应的多立方体。经过测试，对于不同质量的网格，本文算法均能生成满足测量标准的体映射。 图7 不同网格质量 Different mesh qualityFig 7((a)input mesh 1; (b)output mesh 1;(c)input mesh 2;(d)output mesh 2) 3.2.2不同网格类型基于映射的方法(Su等，2019)将不同类型的网格作为映射目标进行大量实验。本文算法对这些不同类型的网格也进行大量测试。使用现有算法生成的多立方体(Fu等，2016)、球(Hu等，2018)以及自由边界(Yang等，2019)作为边界点的位置约束。实验结果表明，对于基于映射方法失败的模型，本文方法同样能够生成良好的体映射。图 8展示了对dancer-children、uu-memento、ramses、sphinx、airplane和santa等不同目标的实验结果。 图8 不同网格类型及更多结果展示 Different mesh types and more resultsFig 8((a) input original meshs; (b) output results) 3.2.3不同位置约束及初始化一个合格的体映射生成算法能够适应不同的位置约束，并且能够处理不同的初始映射。本文将同一个模型映射到不同的目标区域进行实验，如图 9所示，对同一个网格构建目标形状不同的体映射。结果表明，对于不同的目标区域，本文算法均能生成较好的体映射。 图9 不同位置约束 Different position constraintsFig 9((a) input mesh 1; (b) position constraint 1;(c) position constraint 2;(d) position constraint 3;(e) input mesh 2;(f) output mesh 1; (g) output mesh 2;(h) output mesh 3) 此外，实验中本文使用了不同的体映射初始化，包括3种方式：1)谐波映射；2)将网格的内部顶点放置在同一位置；3)随机放置网格内部的顶点位置。图 10展示了同一个模型3种不同的初始化。可以看出，本文算法得到的结果不存在翻转的四面体，能够满足体映射测量标准的要求。 图10 不同初始映射 Different initial mappingsFig 10((a) input mesh; (b) harmonic mapping; (c) same position; (d) random position; (e) position constraint; (f) output mesh 1; (g) output mesh 2;(h) output mesh 3) 3.2.4鲁棒性分析本文算法以基于映射的方法为基础。基于映射的方法对大量数据进行测试，鲁棒性较好。即使对于不能完全消除翻转的模型，该方法仍能消除大部分翻转，保证了本文的雅可比引导是比较合理的。通过雅可比矩阵引导的变形算法处理，能够使映射尽可能地逼近位置约束，便于后面包含位置约束的能量优化。并且本文算法在变形过程中还可以使用线搜索的方法防止产生翻转。因此，对于生成无翻转的体映射，本文方法具有更好的鲁棒性。3.2.5效率分析与现有的体映射方法相比，本文算法在消除翻转方法的基础上添加了雅可比引导的变形算法，因此耗费时间更长。最耗费时间的地方主要在优化变形网格以满足位置约束这一步。为了使优化顺利进行，并尽可能满足位置约束，本文算法将位置约束能量的权重逐渐增大并进行多次优化。这些步骤保证了本文算法能够生成高质量的无翻转的体映射，同时也通过对变形算法进行加速，使算法的效率在一个合理的范围内。3.3与其他方法对比 3.3.1与扭曲有界方法的对比扭曲有界方法是消除翻转这类方法中一种比较有效的方法。扭曲有界的方法不依赖初始映射，通过将映射结果投影到扭曲有界的空间范围内消除翻转，实验中主要与这类方法中的大规模有界扭曲映射(large-scale bounded distortion mappings，LBDM)(Kovalsky等，2015)方法进行比较，对比结果如图 11所示。 图11 本文算法与扭曲有界方法对比 Comparison between our method and LBDMFig 11((a) input mesh; (b) LBDM; (c) ours) 由于扭曲有界的方法比较依赖扭曲界线的选取，一个比较差的界线可能导致LBDM方法产生振荡。而本文算法以基于映射的方法作为基础，使用基于映射方法中的最大扭曲作界线。为了更好地消除翻转，对LBDM方法进行最大1 000次迭代的实验或能够完全消除翻转时停止迭代。通过实验发现，对较复杂的模型，LBDM方法进行1 000次迭代不能完全消除翻转，仍有少部分区域存在翻转的四面体(图 11(b)上)。而对于映射前后差异较大的网格，LBDM方法处理后留存的翻转的四面体数目还有很多，仍有183个翻转的四面体(图 11(b)下)。相比于LBDM方法，本文算法不仅迭代次数更少，并且保证不会产生翻转(图 11(c))。3.3.2与基于映射方法的对比基于映射的方法不依赖无翻转的初始映射，具有较高的鲁棒性，并且也通过后处理降低了体映射的扭曲，与本文目标相似。但通过实验发现，对于某些复杂模型，基于映射的方法往往需要多步迭代才能收敛，并且有些模型在经过多步迭代后仍不能完全消除翻转。而使用本文算法对大量复杂模型进行实验，均能生成无翻转的体映射。本文算法与基于映射的方法主要从3个方面进行对比。1)无翻转是体映射的最基本要求，生成无翻转的体映射是本文的核心工作。如图 12所示，基于映射的方法未能完全消除翻转，在进行了多次迭代之后仍存有少量的翻转四面体(图 12(b))。本文算法通过使用雅可比引导的变形方法，保证变形过程中不产生翻转，因此最终结果不会有翻转的四面体(图 12(c))。2)基于映射的方法可以保证严格的位置约束。本文算法虽然放松了位置约束，但通过后续位置约束的优化，可以将定义的位置约束能量控制在很小的值。3)在时间复杂度上，本文算法需要使用基于映射方法的结果作为引导，因此需要更多的时间。但是对于复杂模型，基于映射的方法也需要多次迭代来消除翻转。由于本文算法不需要完全无翻转的雅可比矩阵引导，因此可以提前终止基于映射方法的迭代，从而减少算法的时间。此外，基于映射的方法也使用了类似的后优化处理，两种方法均能生成扭曲较低的映射。 图12 本文算法与基于映射的方法对比 Comparison between our method and FFFig 12((a) input mesh; (b) FF; (c) ours) 本文算法在本文各图示模型的顶点和四面体数目、初始映射的翻转数目以及每个模型的扭曲能量${\delta _t}$和位置能量$\rho $如表 1所示。可以看出，本文算法不仅不会产生翻转，而且能够满足用户能量和位置的要求。表1 模型数据统计 模型 顶点 四面体 输入翻转 δt ρ 图 1(d) 18 189 69 048 394 7.33 5.51E-18 图 2(a) 16 509 60 600 741 7.20 1.80E-12 图 2(b) 14 653 56 550 253 26.58 5.41E-8 图 3(d) 21 599 81 114 348 7.59 45.82 图 3(f) 21 599 81 114 348 6.83 1.05E-12 图 4(d) 24 519 94 441 19 768 10.01 2.63E-14 图 5(b) 6 888 24 483 228 9.34 38.17 图 5(c) 6 888 24 483 228 9.27 1.67E-4 图 5(d) 6 888 24 483 228 9.01 4.95E-12 图 6(b) 18 051 69 115 265 50.30 9.88E-8 图 6(d) 18 051 69 115 265 30.69 9.88E-8 图 7(b) 18 428 71 512 293 6.61 1.45E-12 图 7(d) 23 398 96 948 287 6.60 1.86E-13 图 8(b)dancer-children 16 433 59 769 2 662 7.46 6.42E-19 图 8(b)uu-memento 15 293 54 886 1 268 6.83 3.48E-17 图 8(b)ramses 19 729 76 850 8 454 14.56 2.12E-18 图 8(b)sphinx 15 534 61 186 4 092 9.35 7.42E-13 图 8(b)airplane 27 610 106 240 33 117 8.45 1.55E-11 图 8(b)santa 26 800 102 629 31 114 6.15 6.69E-14 图 9(f) 35 122 141 171 5 937 7.83 2.67E-22 图 9(g) 35 122 141 171 6 983 8.78 6.85E-18 图 9(h) 35 122 141 171 24 250 11.56 1.29E-12 图 10(f) 30 493 123 576 853 6.86 5.55E-13 图 10(g) 30 493 123 576 68 433 6.86 1.08E-11 图 10(h) 30 493 123 576 6 938 6.86 5.42E-13 图 11(c)elephant 13 188 48 623 2 579 7.46 2.27E-19 图 11(c)ant 12 898 46 172 9 908 24.05 1.05E-31 图 12(c)bear 27 626 107 460 1 055 7.83 4.25E-12 图 12(c)hand 21 416 92 369 50 405 13.98 6.72E-11 Model data statisticsTable 14结论本文提出了一种新的计算无翻转体映射的方法，不依赖无翻转的初始映射，对不同质量、不同目标的网格均有较好表现，并且能够处理现有方法处理失败的复杂网格。本文算法的优势在于，使用雅可比矩阵作为引导进行变形操作，通过构建合理的雅可比矩阵，可以保证变形后的网格不存在翻转的四面体，并且使变形网格尽量接近目标形状。同时通过控制位置能量来满足给定的位置约束，并通过后优化处理进一步降低映射的扭曲。但是，本文算法有一定的局限性。首先，本文算法依赖基于映射的方法，且需要进行多步的优化处理，因此算法时间复杂度较高。其次，本文算法可以理论上保证映射是无翻转的，大量实验也表明本文算法生成的映射都能满足位置约束的要求，但目前没有理论保证映射一定会满足给定的位置约束。