research-article

Open access

Geometric Optimisation Via Spectral Shifting

Authors:

Roman Poya,

Rogelio Ortigosa,

Theodore KimAuthors Info & Claims

ACM Transactions on Graphics, Volume 42, Issue 3

Article No.: 29, Pages 1 - 15

https://doi.org/10.1145/3585003

Published: 07 April 2023 Publication History

All formats PDF

Abstract

We present a geometric optimisation framework that can recover fold-over free maps from non-injective initial states using popular flip-preventing distortion energies. Since flip-preventing energies are infinite for folded configurations, we propose a new regularisation scheme that shifts the singular values of the deformation gradient. This allow us to re-use many existing algorithms, especially locally injective methods for initially folded maps. Our regularisation is suitable for both singular value- and invariant-based formulations, and systematically contributes multiple stabilisers to the Hessian. In contrast to proxy-based techniques, we maintain second-order convergence. Compact expressions for the energy eigensystems can be obtained for our extended stretch invariants, enabling the use of fast projected Newton solvers. Although spectral shifting in general has no theoretical guarantees that the global minimum is an injection, extensive experiments show that our framework is fast and extremely robust in practice, and capable of generating high-quality maps from severely distorted, degenerate and folded initialisations.

1 Introduction

Generating fold-over free maps is a common goal in geometry processing and physical simulation. Nonlinear techniques for generating locally injective maps often require a valid starting point (such as a Tutte [1963] embedding) as many popular distortion energies become infinite if the initial domain is not locally injective. Such energies have been studied in continuum mechanics since the foundational work of Ball [1976, 1981, 1983] on invertibility theorems.

Coercivity (i.e., barriers) is a key concept in flip-preventing energies which dictates that the energy should go to infinity when the map inverts [Hartmann and Neff 2003]. A series of recent works in computer graphics have proposed such barrier energies [Fu et al. 2015; Rabinovich et al. 2017; Schüller et al. 2013; Smith and Schaefer 2015], while in computational physics, the topic of recovering local injectivity has been studied since the work of Winslow [1966] and its variational formulation proposed later in Brackbill and Saltzman [1982]. We refer the reader to Garanzha et al. [2021a] for a historical overview and Fu et al. [2021] for a detailed survey.

In this paper, we propose an optimisation framework for generating fold-over free mappings using popular rotation invariant flip-preventing energies, even when the initial domain is non-injective or degenerate. Our work uses a continuum mechanics approach to introduce the following contributions:

•

A regularisation technique based on shifted singular values that yields a simple, consistent and general formulation.

•

Compact eigensystems for regularised energy Hessians which can be embedded into fast second-order solvers, using extended, regularised stretch invariants [Smith et al. 2019] with a tailor-made optimisation scheme.

•

Between \({\bf 1.43\times }\) to \({\bf 12.9\times }\) mean speedups over existing methods on a large-scale benchmark.

More specifically, we solve for the regularisation parameter ahead of minimisation (in an outer loop; see Algorithm 2) and run a series of Projected Newton (PN) steps that make use of energy eigensystems. Once the map is repaired, our workflow reduces to well-known locally injective methods. This property then enables the re-use of established techniques for tackling global injectivity, while still starting from folded configurations (see Figure 4). While there are no theoretical guarantees that the global minimum of the energy corresponds to an inversion-free map, our framework passes large-scale datasets with up to 100% success rate. To aid reproducibility, we provide kernel computation code in the supplement.

2 Related Work

2.1 Problem Statement

To facilitate our discussion in the next section, we will first articulate the problem from a mechanics perspective [Bonet et al. 2016b]. Given a mesh of nodes \(\mathcal {V}\) and elements \(\mathcal {T}\), potentially containing fold-overs, our goal is to minimise energies with barriers such that¹

\begin{align} \Pi (\boldsymbol {x}^{\ast }) & = \min _{\boldsymbol {x}} \int _{\Omega _0} \mathcal {D}(\boldsymbol {F}) \; \mathrm{d} \Omega \end{align}

(1a)

\begin{align} J & = \text{det}\boldsymbol {F} \gt 0\qquad \end{align}

(1b)

where \(\Pi\) is the total potential energy of the system, \(\boldsymbol {x}\) is the vertex positions in the mesh, \(\boldsymbol {x}^\ast\) the optimal solution, \(\boldsymbol {F}\) the deformation gradient tensor, \(\mathcal {D}(\boldsymbol {F})\) the distortion energy density function, and J the determinant (also known as the Jacobian). Second-order methods require both first and second derivatives of the integral in (1a):

\begin{align} D\Pi [\delta \boldsymbol {u}] & = \int _{\Omega _0} \frac{\partial \mathcal {D}(\boldsymbol {F})}{\partial \boldsymbol {F}} : \boldsymbol {\nabla _w} \delta \boldsymbol {u} \; \mathrm{d} \Omega \end{align}

(2a)

\begin{align} D^2\Pi [\delta \boldsymbol {u}, \Delta \boldsymbol {u}] & = \int _{\Omega _0} \boldsymbol {\nabla _w} \Delta \boldsymbol {u} : \frac{\partial \mathcal {D}^2(\boldsymbol {F})}{\partial \boldsymbol {F} \partial \boldsymbol {F}} : \boldsymbol {\nabla _w} \delta \boldsymbol {u} \; \mathrm{d} \Omega \end{align}

(2b)

where \(:\) denotes the inner product of matrices (\(\boldsymbol {A} : \boldsymbol {B} = \mathrm{tr}(\boldsymbol {A}^T\boldsymbol {B})\)), \(\boldsymbol {\nabla _w}\) is the gradient operator with respect to a desired domain (such as the initial mesh) and \(\delta \boldsymbol {u}\) and \(\Delta \boldsymbol {u}\) are finite element trial and test functions. The first term in integrand (2a) is the gradient of energy, also called the first Piola-Kirchhoff stress tensor (PK1), which we express as

\begin{equation} \boldsymbol {P} = \frac{\partial \mathcal {D}(\boldsymbol {F})}{\partial \boldsymbol {F}}. \end{equation}

(3)

The middle term in integrand (2b) is the Hessian of energy, which we denote as

(6)

2.2 Inversion-free Mappings

In graphics, methods for generating inversion-free mappings include Tutte [1963], Bounded Distortion (BD and LBD) [Aigerman and Lipman 2013; Kovalsky et al. 2015], extremal quasi-conformal mapping [Weber et al. 2012], Simplex Assembly (SA) [Fu and Liu 2016], Adaptive Block Coordinate Descent (ABCD) [Naitsat et al. 2020], Total Lifted Content (TLC) [Du et al. 2020a] and Fold-over Free Maps (FFM) [Garanzha et al. 2021a; Su et al. 2019]. Progressive embeddings and intrinsic triangulations have also been successfully applied in 2D (parametrisation) contexts [Campen et al. 2021; Gillespie et al. 2021; Liu et al. 2018; Shen et al. 2019], which relax the requirements on the input triangulation to guarantee finding a discrete, locally injective map. To our knowledge, TLC and FFM are the only frameworks that pass the whole local injectivity benchmark of Du et al. [2020b] in 2D and 3D.

To guarantee inversion-free maps, Equation (1b) must be satisfied, so one common approach is to perform an unconstrained minimisation of barrier energies that implicitly encode this constraint. However, the barrier prevents these energies from being applied when the initialisation contains fold-overs (\(J \lt 0\)). Consequently, state-of-the-art algorithms either lift the space of simplices to a logical space by minimising some simple field [Du et al. 2020a], or consider a metric that bounds the amount of distortion [Aigerman and Lipman 2013; Fu and Liu 2016; Kovalsky et al. 2015]. Most of these approaches change the metric (transformation matrix) by employing an auxiliary strictly non-degenerate map.

Minimising barrier energies with fold-overs without changing the triangulation/tetrahedralisation relies on either filtered folded elements [Naitsat et al. 2020] or regularised Jacobians [Garanzha et al. 2021a]. Filtering then requires adding energies for inverted elements that can unnecessarily compete with the primary distortion energy. Alternatively, Jacobian regularisation allows the same energy to be used, but leads to non-convex energies and inconsistent treatment of the map’s metric, which implies that the regularised Jacobian cannot be retrieved from the product of singular values of the deformation gradient. This complicates the use of state-of-the-art second-order approaches [Shtengel et al. 2017; Smith et al. 2019], which exploit the structure of the deformation gradient.

Invertible distortion energies have also been studied in simulation [Irving et al. 2004; Smith et al. 2018; Stomakhin et al. 2012], but these formulations tackle a different problem. They are neither designed to solve inversion-free mappings nor are they a substitute for locally injective methods (further discussion in Section 5.8).

2.3 Energy Minimisation

Minimising barrier energies typically require stabilisation of the Hessian to make them semi-definite positive. We use the term stabilisation here as opposed to regularisation because Hessian stabilisation is distinct from Jacobian regularisation. The former is necessary to obtain a descent direction while the latter makes a barrier energy finite under fold-overs. To overcome the problem of Hessian indefiniteness, one can clamp the negative eigenvalues of the Hessian [Fu and Liu 2016]. Alternatively, a proxy Hessian such as a Laplacian [Gargallo-Peiró et al. 2015], or a linearised version of the same Hessian [Poya et al. 2016] can be used/added, or a semidefinite-indefinite decomposition can be performed [Garanzha et al. 2021a].

Projected Newton (PN) [Stomakhin et al. 2012; Teran et al. 2005] is not a good candidate for certain cases of untangling, as it can lead to damped convergence, and coupled with a non-filtered line search, can even stall. Over 10,000 rounds of L-BFGS or Projected Newton are commonly needed [Du et al. 2020a; Garanzha et al. 2021a]. Alternatively, Quasi-Newton (QN) methods can smooth away barrier terms and cause the optimisation to converge to the incorrect (still folded) minima. As reported in Du et al. [2020a], the QN method in TLC can fail to recover local injectivity in moderately complex problems, even though the energy is guaranteed to achieve its global minimum at a locally injective state.

3 Our Regularisation Approach

We will present our overall approach by starting with continuum mechanics preliminaries (Section 3.1), and present our core regularisation scheme (Section 3.2). From this, we will show how to build a set of regularised stretch invariants (Section 3.3), which allow us to derive compact expressions for their Hessians. After describing a strategy for setting a key regularisation parameter (Β) in Section 3.4, we will apply our regularisation scheme to several common energies (Section 3.5), and end with a brief discussion (Section 3.6).

3.1 Continuum Mechanics Preliminaries

In the context of mechanics, the mapping between a mesh with inverted elements and an inversion-free mesh [Aigerman and Lipman 2013; Du et al. 2020a; Fu and Liu 2016; Garanzha et al. 2021a; Kovalsky et al. 2015; Su et al. 2019] can be formalised as the following Arbitrary-Lagrangian-Eulerian (ALE) scheme.

Consider an element in the auxiliary (“ideal”) configuration \(\Omega _D\), with position vectors \(\boldsymbol {w}\) and a desired shape and size. A corresponding element \(\Omega _0\) in the physical configuration has position vectors \(\boldsymbol {y}\). This element undergoes motion to arrive at the final physical configuration \(\Omega\), denoted with position vectors \(\boldsymbol {x}\) [Bonet et al. 2016b].

To measure distortion we consider the affine mapping \(\boldsymbol {\phi }\), from \(\Omega _D\) to \(\Omega\) (see Figure 2). This can be decomposed into \(\boldsymbol {\phi }_D : \Omega _N \rightarrow \Omega _D\) from the natural element \(\Omega _N\) to the auxiliary element \(\Omega _D\), and \(\boldsymbol {\phi }_P : \Omega _N \rightarrow \Omega _0\), from the natural element \(\Omega _N\) to the physical element \(\Omega _0\). Thus, \(\boldsymbol {\phi }\) can be determined by the composition \(\boldsymbol {\phi } = \boldsymbol {\phi }_{P} \circ \boldsymbol {\phi }^{-1}_D\).

Following nonlinear mechanics [Bonet et al. 2016b], the fibre map transforms an edge between an auxiliary and physical element:

\begin{align} \boldsymbol {F} = \boldsymbol {\nabla _w}\boldsymbol {x} = (\boldsymbol {\nabla _\xi }\boldsymbol {x}) (\boldsymbol {\nabla _\xi }\boldsymbol {w})^{-T}. \end{align}

(5)

The domain \(\Omega _D\) can be a single element, a polygon soup, or a fully conformal mesh. In the case of a conformal mesh, \(\boldsymbol {F}\) is the deformation gradient tensor \(\boldsymbol {F} = \boldsymbol {F_y}\), where \(\boldsymbol {w}=\boldsymbol {y}\) is the position vector of the initial configuration.

Following Bonet et al. [2015, 2016a] the three fundamental distortion measures are \(\boldsymbol {F}\) (edge map), its cofactor \(\boldsymbol {H}\) (area map) and its determinant (volume map) such that (see Figure 3)

\begin{align} \boldsymbol {F} = \boldsymbol {\nabla _w}\boldsymbol {x}, \;\;\; \boldsymbol {H} = \mathrm{cof}\boldsymbol {F} = J\boldsymbol {F}^{-T} = \frac{1}{2} \boldsymbol {F} \times \boldsymbol {F}, \;\;\; J = \mathrm{det}\boldsymbol {F}, \end{align}

(6)

where \(\times\) is the tensor cross product between two second order tensors expressed as \([\boldsymbol {A} \times \boldsymbol {A}]_{iI} = \mathcal {E}_{ijk}\mathcal {E}_{IJK} A_{jJ}B_{kK}\) where \(\mathcal {E}\) is the third order Levi-Civita tensor [de Boer 1982]. Furthermore, \(\boldsymbol {F}\) admits a left polar decomposition

\begin{align} \boldsymbol {F} = \boldsymbol {U} \boldsymbol {\Sigma } \boldsymbol {V}^T = \boldsymbol {RS}, \end{align}

(7)

into a rotation tensor \(\boldsymbol {R} = {\boldsymbol {U}} {\boldsymbol {V}}^{T}\) and a stretch tensor \(\boldsymbol {S} = {\boldsymbol {V}} \boldsymbol {\Sigma }{\boldsymbol {V}}^{T}\). Here, \({\boldsymbol {U}}\), \(\boldsymbol {\Sigma }\) and \({\boldsymbol {V}}\) are the tensors from the signed singular value decomposition (SSVD) of \(\boldsymbol {F}\) [Irving et al. 2004; Kovalsky et al. 2015].

3.2 Regularisation of F

We will now regularise the deformation gradient tensor \(\boldsymbol {F}\), which in turn will allow us to formulate regularised invariants which we can then use to regularise arbitrary barrier energies. Recall that a tensor invariant (e.g. the determinant) is a scalar function of a tensor that remains unchanged under a basis change or a rotation of the target frame. We use the subscript r to denote a regularised quantity, e.g. \(\boldsymbol {F}_r\) or \(\sigma _{r_i}\). We regularise \(\boldsymbol {F}\) by spectrally shifting its signed singular values:

\begin{align} \boldsymbol {F}_r & = {\boldsymbol {U}} (\boldsymbol {\Sigma } + \boldsymbol {\Sigma }_s) {\boldsymbol {V}}^T = {\boldsymbol {U}} \boldsymbol {\Sigma }_r {\boldsymbol {V}}^T \end{align}

(8a)

\begin{align} \mathrm{with} \qquad \sigma _{r_i} & = \sigma _i + \beta \quad i=1,2,3\qquad\qquad\qquad\ \ \ \ \end{align}

(8b)

\begin{align} & \mathrm{s.t.\; det}\boldsymbol {F}_r \gt 0. \end{align}

(8c)

Above, \(\boldsymbol {\Sigma }_s = \beta \boldsymbol {I}\) is a strictly positive definite regularisation tensor. The tensor \(\boldsymbol {\Sigma }_r\) then has a smallest signed singular value of at least \(\sigma _i + \beta \ge \epsilon\) where \(\epsilon\) is a small positive number that we iteratively decrease over the course of optimisation (see \(\epsilon _k\) in Algorithm 2). We set \(\epsilon = 10^{-6}\) as the minimum value.

We consider Β a constant during minimisation; this assumption is discussed extensively in Section 3.4. We could alternatively regularise only one signed singular value using the convention (\(\mathrm{sign}(\sigma _3) \lt 0\)) [Irving et al. 2004], but such an approach prohibits the derivation of generic regularised energies. The complexity of the expressions becomes unwieldy [Stomakhin et al. 2012], as regularisation in one principal direction introduces anisotropy into isotropic energies. Geometrically, regularising \(\boldsymbol {F}\) scales-up the deformation metric i.e., increases the Jacobian (volume ratio), edge length, and area ratios (see Figure 3). However, as soon as fold-overs are repaired, we fall back to the original metric.

3.3 Regularisation of Stretch Invariants

We can now use the regularised deformation gradient \(\boldsymbol {F}_r\) to formulate regularised stretch invariants. We first express an arbitrary distortion energy \(\mathcal {D}\) in terms of an extended set of stretch tensor invariants:

\begin{align} \boxed{ \mathcal {D}(\boldsymbol {F}) = \hat{\mathcal {D}} (I_F, II_F, I_H, II_H, J) } \end{align}

(9)

with the following definitions for invariants in 3D and 2D

(10)

Two additional cofactor invariants, \(I_H\) and \(II_H\), are introduced on top of those from Smith et al. [2019], and are “cross” terms that model area distortion in 3D (see Figure 3). They are redundant in 2D, as \(I_H=I_F=\sigma _1 + \sigma _2\) and \(II_H=II_F=\sigma _1^2 + \sigma _2^2\), but in 3D they can be expressed as

\begin{align} \ \ I_H & = \frac{1}{2} (I^2_F - II_F) & \text{area-rotation in 3D}\quad\ \end{align}

(11a)

\begin{align} II_H & = I^2_H - 2I_F\,J & \text{area-preservation in 3D} \end{align}

(11b)

We discuss these two cofactor invariants further in Section 3.6.

Using Equation (8), we can now formulate a set of regularised invariants,

(12)

where \(III_F = \sigma _1^3 + \sigma _2^3 + \sigma _3^3\) is an auxiliary invariant.

These invariants can be plugged in to any existing distortion energy to yield a regularised counterpart. The gradients and Hessian of these regularised invariants, as well as the exact stabilisers that they contribute to the Hessian, are given in the supplement. Critically, our choice of invariants allows us to derive compact eigensystems for the regularised energy Hessians.

3.4 Finding the Regulariser, Β

The decision to regularise using shifted singular values can seem empirical. However, we will show that all state-of-the-art approaches, including TLC [Du et al. 2020a], FFM [Garanzha et al. 2021a] and ABCD [Naitsat et al. 2020] can be understood as similar regularisers. Our scheme will then combine the merits of these approaches.

3.4.1 State-of-the-art Approaches.

We start with TLC [Du et al. 2020a], where the value of regulariser emerges from a strictly non-degenerate (lifted) transformation. The regulariser is part of the TLC energy itself and cannot be applied to other energies.

In contrast, FFM [Garanzha et al. 2021a] performs a penalty-style regularisation of J such that

\begin{align} J_r(J) = \frac{1}{2}\left(J + \sqrt {J^2 + \delta ^2}\right) \end{align}

(13)

for every element for which \(J \lt 0\). The \(\delta\) is a positive constant that typically encodes the minimum Jacobian value over the entire mesh. This approach allows the same energy to be used for folded elements, but introduces metric inconsistency.² This inconsistency implies that: (a) \(J_r\) does not correspond to the product of singular values of any metric (\(\boldsymbol {F}\) or \(\boldsymbol {F}_r\)) and (b) results in non-convex energies, even when applied to convex energies. Moreover, it becomes impossibles to formulate the problem in terms of singular values and obtain Hessian eigensystems. Our regularised invariants do not suffer from these problems.

Finally, ABCD [Naitsat et al. 2020] considers a separate flip term to modify non-positive Jacobians, which is then used in conjunction with a desired distortion energy:

\begin{align} J_r(J) = \prod \limits _i^{d-1}\sigma _i \delta , \end{align}

(14)

where negative Jacobians are replaced with a positive constant \(\delta\).

3.4.2 Our Approach.

We can combine the strengths of the aforementioned approaches. In general, given a regularised Jacobian of, e.g. Equations (13) or (14), we can solve the last line in Equation (12) for Β. This allows for both Jacobian and singular-value regularisations to be embedded into a unified approach (see Algorithm 2).

In 2D, this amounts to solving a quadratic equation (\((\sigma _1 + \beta)(\sigma _2 + \beta)=J_r\)), and a cubic in 3D (\((\sigma _1 + \beta)(\sigma _2 + \beta)(\sigma _3 + \beta)=J_r\)). This is achieved by polynomial root finding on Line 12 of Algorithm 1.

Alternatively, Β can be computed for every element, since the last line of Equation (12) can be solved for Β analytically. In 2D, we obtain \(\beta = \frac{1}{2}(-I_F + \sqrt {I^2_F + 4 (J_r - J)})\), which is always positive and can be back-substituted to carry out the energy derivatives.

We have now fully parametrised over Β, where Β can originate from any kind of constraint, not just those described here. In contrast to Garanzha et al. [2021a], a constant Β in 2D amounts to a linear transform of a standard energy into a regularised energy (note that, \(I_F\) is linear in \(\boldsymbol {F}\) and has SPD Hessian), and linear transforms preserve convexity. With this set of regularised invariants, we have established a consistent way of expressing infinite energies over injective, non-injective and degenerate mappings.

3.5 Regularised Energies

We can now regularise several well-known distortion energies and give explicit expressions for their eigensystems. As mentioned, we can either pre-compute Β ahead of minimisation or re-substitute its analytical form and compute energy derivatives. The latter is particularly useful in 2D energies such as MIPS or Symmetric Dirichlet (SD) where Β does not explicitly appear in the energy (c.f. Table 1), and only the Jacobian regulariser \(\delta\) from Equations (13) or (14) are needed for a consistent formulation. In such cases, our approach becomes a generalisation of FFM. In the supplement we work out the MIPS energy with such a \(\delta\) regulariser. We also follow this setup whenever possible (back-substituion/per-element Β for 2D MIPS and SD and a global Β for the rest) for numerical experimentation.

Table 1.

For the case when Β is a constant, the 2D energies, gradients and analytic Hessian eigensystems are given in Table 1 for Regularised MIPS (RMIPS), Regularised Symmetric Dirichlet (RSD) and Regularised Symmetric ARAP (RSARAP). The Hessian expressions for the 3D cases are given in the supplement, as well as a script to generate them symbolically. We follow Smith et al. [2019] notation.

3.6 Discussion of Invariant-based Formulation

Distortion energies are typically expressed in terms of singular values [Xu et al. 2015], Cauchy-Green invariants [Bonet et al. 2016b], polyconvex invariants [Bonet et al. 2015, 2016a],³ or stretch tensor invariants [Smith et al. 2019].

Stretch tensor invariants and singular value formulations can both reveal the full eigensystem of the Hessian [Chen and Weber 2017; Smith et al. 2019; Zhu 2021], but invariant formulations generally have the advantage that they encode specific distortion metric such as length, area, volume. While these metrics can also be re-expressed in terms of singular-values, invariant-based formulations additionally allow convexity to be studied in closed-form, see Equation (15). By adopting an invariant-based formulation, we can re-use nearly a century of well-established literature [Ball 1976; Bonet et al. 2016b]. Without loss of generality, we chose to work with stretch tensor invariants because it allowed us to derive analytic eigensystems for our regularised energy Hessians.

4 Optimisation Approach

We use a numerical continuation technique, coupled with a backtracking scheme [Nocedal and Wright 2006] inspired by FFM, that runs a series of minimisation by decreasing an auxiliary penalty parameter \(\epsilon _k\) starting from \(\epsilon _0 = 1\) (see line 9 of 2nd procedure in Algorithm 2). This parameter is also called the finite untangling sequence parameter. We compute Β (Algorithm 1) in an outer loop so it remains constant for all elements and minimisation iterations. As shown in Algorithm 2, Β can be computed using a wide range of techniques such as singular-value regularisation [Liao et al. 2021; Naitsat et al. 2020], Jacobian regularisation using a well-known empirical formula (line 11 of 2nd procedure) [Gargallo-Peiró et al. 2015], or using the finite untangling sequence [Garanzha et al. 2021a]. We have found the latter two result in faster convergence. In general, we perform the following (in case of per-element/analytical Β, step 2 is not necessary and computing \(\delta\) in Algorithm 2 is sufficient; see supplement for C++ code):

(1)

At every round (outer loop in Algorithm 2) compute the worst Jacobian, and then compute \(\delta\) by setting \(\epsilon _0=1\).

(2)

Compute the maximum regulariser Β over the entire mesh.

(3)

Perform a full PN minimisation (inner loop in Algorithm 2):

(a)

For elements with \(J \le 0\) compute the regularised invariants, energy, PK1 and filtered Hessian.

(b)

For elements with \(J \gt 0\) use the standard (non-regularised) invariants and compute the energy, PK1 and filtered Hessian.

(4)

At the end of every minimisation, if \(J_{\mathrm{min}} \le 0\), recompute Β using \(\epsilon _k\) based on a sufficient energy decrease, i.e., \(\epsilon _k = \epsilon _{k-1}(1 - \mathcal {D}^{k}/\mathcal {D}^{k-1})\) and \(\delta\), and repeat.

4.0.1 Non-injective States.

In a non-injective state, it is not necessary to drive the minimisation to gradient- or energy-based convergence. We show a simple example in Figure 6. Two opposite boundary vertices are moved inwards and optimised for two minimisation rounds. We show convergence of the first minimisation round, where displacements drop quickly despite a high residual. To detect such cases, we first determine if the number of fold-overs have remained unchanged in three consecutive iterations, and whether the displacements have dropped below minimum characteristic length of the mesh. The characteristic length is computed by scaling minimum edge lengths with the ratio of energy decrease. Similar and better alternatives can be used here [Naitsat et al. 2020; Zhu et al. 2018]. Furthermore, when Β is assumed constant, the positivity of \(J_r\) may not be guaranteed since new negative Jacobian may appear during iteration that exceed the value of precomputed Β. In such cases, we assign the absolute minimum energy over the entire mesh to the corresponding flipped elements that exceed Β, and filter their contributions from the line search.

5 Numerical Experiments

In this section, we provide results obtained with our framework and compare against other algorithms. For most studies, the code for other frameworks were obtained either from the authors or official project page. For fairness, we also implemented a per-element, numerically projected Newton, similar to the PN solvers used in Simplex Assembly [Fu and Liu 2016] and TLC-PN [Du et al. 2020a].

For completeness, we also implemented the Hessian modified version of FFM which performs a semidefinite-indefinite decomposition as presented in Garanzha et al. [2021a] (see Appendix A.2 in that paper and Appendix B here). We call this version Modified-FFM. The paper presented but did not implement the approach, and their final version of FFM used an L-BFGS solver. A direct comparison with other code bases give us approximate performance comparisons, but a clearer picture emerges by comparing our framework’s results to our own implementation of other techniques such as Modified-FFM and numerical PN. Unless otherwise stated, all analyses were performed single-threaded on an Intel Xeon(R) E5-2650. We use MKL PARDISO, Eigen [Guennebaud and Jacob 2010] and Fastor [Poya et al. 2017] for tensor manipulations.

5.1 Surface Parametrisation

5.1.1 Comparing with FFM.

We begin by comparing against FFM. We include our Modified-FFM. We chose two boundary surface parametrisation problems [Liu et al. 2018] with a small number of folded elements. For comparison with FFM, we used the same termination tolerance for both our analytic PN solver and FFM’s L-BFGS solver (\(10^{-4}\)) with an outer loop tolerance of \(10^{-5}\). For fairness, we used the same \(\epsilon _k\)-decay formula (finite untangling sequence [Garanzha et al. 2021a]) and back-substituted the analytical Β into energy derivatives to ensure the differences were only solver related. The volumetric term in FFM was chosen as \(1/128\) which corresponded well with our \(\boldsymbol {F}\)-based MIPS. We confirmed that the initially computed Jacobians matched.

Figure 5 shows comparison of our framework with FFM and Modified-FFM. All three repaired fold-overs successfully, but our regularisation achieved the same quality almost 5\(\times\) faster. Modified-FFM almost degraded to a first order in the subsequent iteration of outer loop but performed better than FFM.

5.1.2 Comparison with TLC.

To compare against TLC, we used the David-to-H local injectivity benchmark from Du et al. [2020a] (Figure 7). We chose the regularised MIPS energy and symbolically generated a regularised Advanced MIPS (AMIPS) energy for quasi-isometric mapping. This adds a volumetric term to MIPS which can be written in terms of regularised stretch invariants as

\begin{align*} \hat{\mathcal {D}}_{\mathrm{RAMIPS}}(II_F, J) = \mu \frac{II_F}{2J_r} + \mu \left(1- \frac{J}{J_r}\right) + \lambda \big(J_r + J_r^{-1}\big) - (2\lambda +\mu) \end{align*}

where \(\mu\) and \(\lambda\) are, respectively, shape and area/volumetric coefficients.

We ran experiments for \(\lambda = 0\) (MIPS) and \(\lambda =1\), and set \(\mu =1\) in both cases. All solvers were stopped as soon as local injectivity was recovered. Both FFM (with \(\lambda =0\)) and TLC’s QN solvers failed to produce a valid “all” injective parametrisation for this model. With FFM, regardless of how we tuned \(\epsilon _k\), it was unable to solve the original MIPS energy. FFM with a volumetric term succeeded.

We further ran the Modified-FFM and a numerical PN solver for the case of quasi-isometric mappings (\(\lambda =1\)). In this example, Modified-FFM was faster than FFM, and while PN converged quickly, the cost of the per-element eigen-decomposition dominated. On the other hand, our framework recovered injectivity up to 10\(\times\) faster (Figure 7). While this decomposition is embarrassingly parallel, we ran our benchmarks serially for fairness of comparison.

TLC generated huge distortion with \(J_{\text{min}} = 1.2\times 10^{-20}\), which is \(10^{19}\) times worse than the ones obtained by our RAMIPS energy. Best running times and qualities are highlighted in Figure 7.

5.1.3 Comparison with SA.

Simplex Assembly (SA) [Fu and Liu 2016] can respect user-specified distortion bounds, but loses the robustness of TLC, FFM, or our framework, with respect to repairing fold-overs. SA minimises the same AMIPS energy as FFM, so it is natural to compare it against our regularised MIPS energy.

We chose to analyse quality with the example in Weber and Zorin [2014], where a convex polygon is morphed into a concave shape using a harmonic mapping, resulting in 28 fold-overs (see inset). We also compared the results of our framework using RSD and RSARAP. As seen at the bottom of Figure 8, the quality with respect to J obtained with SA is comparable to our RMIPS and FFM.

While SA creates a few elements with extremely low J, our framework pushes the elements away from this zone. SA’s conformal results are strictly within a bound and better than our RMIPS. Given that RSD and RSARAP create an isometric mapping, the isometric distortion of RSARAP and RSD is lower.

5.2 Locally Injective Mappings Benchmark

We conducted a parametrisation study with SA, ABCD, FFM, TLC, and our regularised energies on the most challenging examples from the recent benchmark of Du et al. [2020b]: the Lucy-to-G and Gargoyle-to-S embeddings (Figure 9). To cover the spectrum from conformal to isometric, we used RMIPS (and FFM with \(\lambda =1/128\), though quasi-isometric maps can also be obtained by increasing the value of \(\lambda\)), RSD and RSARAP energies. The initial embeddings for these models have \(\sim\)98K triangles and over 3K fold-overs, requiring the optimisation to move a large layer of elements into a concave shape. SA and ABCD (PN) were unable to create inversion-free domains, and respectively failed after \(\sim\)300 and \(\sim\)400 seconds. TLC (PN) succeeded, but yielded nearly collapsed elements (\(J\approx 0\)) due to stopping criterion. FFM succeeded but took over an hour on the Lucy model.

Our framework minimised the RMIPS energy almost 75\(\times\) faster than FFM on this model. It successfully minimised RMIPS, RSD, and RSARAP in a few seconds, and yielded high-quality, fixed-boundary, locally injective maps. The results of FFM and our RMIPS energy were visually identical, so only FFM is shown in Figure 9. However, the histograms should give a better representation of quality (J, conformal \(\frac{\sigma _1}{\sigma _2}\) and isometric \(II_F\)) distribution.

We ran our regularised energies over all the 2D and 3D benchmarks of Du et al. [2020b], as shown in Figure 10. Our framework passed the entire benchmark (similar to TLC and FFM). The benchmark covers a variety of 2D parameterisations, 3D polycube parameterisations and 3D deformations, and range from easy to extremely difficult. We chose \(\lambda =0.01\) in 2D (to be compatible with FFM although, this creates low quality elements) and \(\lambda =1\) in 3D for the RAMIPS energy. Unlike Du et al. [2020a]; Garanzha et al. [2021a], we used the significantly more challenging source mesh as the rest configuration, instead of an auxiliary “ideal” element.

We also ran a side-by-side comparison with FFM, Modified-FFM, and ABCD by tuning all solver parameters (e.g., number of steps, termination criteria, and \(\epsilon _k\) decay (for FFM and ours)) to be identical. Figure 11 shows our framework was statistically 12.9\(\times\) faster than FFM over 10,000 2D meshes and 2\(\times\) faster over 900 3D meshes. ABCD did not have 100% success rate (91.1% in 2D and 88.3% in 3D) but, it was the second fastest method in 2D for cases that it succeeded (within 1.8\(\times\) of our analytic PN). The performance of ABCD and our framework confirms that for 2D problems, PN solvers are preferred.

In many 3D cases, ABCD stalled at the repair phase and in general performed over 10\(\times\) slower (excluding failure times) compared to ours. For some 3D cases when the initialisation was close to the rest shape, FFM’s L-BFGS was faster than our framework. However, it had a more non-uniform convergence pattern and dispersed performance profile over the entire dataset. As expected, FFM’s performance was closer to ours in 3D since, unlike L-BFGS, the cost of the linear solver factors in for our framework. The second fastest method in 3D was our implementation of Modified-FFM, which was 43% slower than our analytic PN. This is because Modified-FFM, despite doing more iterations than analytic PN, does not require an SVD which dominates the assembly cost in 3D.

Overall, our framework was strictly faster than other methods on 72% (7728/10706) of the 2D cases and 65% (593/904) of the 3D cases. Such performance disparities between different optimisers have been well documented in the context of locally injective methods [Kovalsky et al. 2016; Rabinovich et al. 2017; Zhu et al. 2018]. Quality-wise, best performing methods were our RSD (see Figure 10 instead) and ABCD, although all frameworks yielded comparable results.

5.3 Volumetric Polycube Parametrisation

Fixed boundary volumetric polycube parametrisation is an ideal candidate for our regularised energies, since fold-overs are common. In Figure 12 we show several examples with many initially folded tetrahedra, which we then optimised with the RMIPS energy. Given the fixed boundaries, we did not optimise the individual faces, which can lead to higher distortion. The exception is the Max Planck model, where we processed the individual faces (the smaller cube in Figure 12) using 2D RMIPS. This yielded a more uniform distribution of conformal distortion along the surface, and improved the quality of initial parametrisation.

5.4 Free Boundary Surface Parametrisation

In Figure 13 we embed multiple 3D models onto star-shaped domains, resulting in \(\sim\)10K fold-overs. We used RSARAP to optimise our map and obtained results similar to previous works [Claici et al. 2017; Rabinovich et al. 2017; Shtengel et al. 2017; Smith et al. 2019] in approximately the same number of iterations as Shtengel et al. [2017] and Smith et al. [2019]. Our RSARAP energy is stabler, and converges quickly despite starting from an invalid Tutte embedding. Since the energy has a barrier, regularising J decreases the residual, which accelerates convergence. These problems each took 10–12 iterations to converge.

In Figure 1, we start from an invalid mapping and use our RSD energy. As described in Section 2, using regularised barrier energies and filtered Hessians within PN solvers can cause oscillatory convergence, (see Figure 13), while unfiltered Hessians lead to stalls. Over a wide range of examples, we observed that our solver did not have this issue.

Fig. 1.

Fig. 2.

Fig. 3.

5.5 Globally Injective Surface Parametrisation

The purpose of our framework is to generate inversion-free mappings; a first step towards local injectivity. However, many frameworks (including ours), cannot ensure a locally injective map due to overwinding issue and k-covering traps (see Figure 4). Overlaps can occur under free boundaries, which means the parametrisation may not be globally injective as well. The authors of TLC fixed this issue with the Smooth Excess Area (SEA) proposal [Du et al. 2021], while FFM proposed using phantom triangulation [Garanzha et al. 2021b].

Fig. 4.

Fig. 5.

Fig. 6.

Fig. 7.

Fig. 8.

Fig. 9.

Fig. 10.

Our regularisation scheme allows for the complete re-use of existing techniques for locally injective mappings. To circumvent both overwinding and boundary self-intersection, we can use a “curve-repulsive” energy on the boundaries, such as those proposed in Smith and Schaefer [2015] and Su et al. [2020]. In Figure 14 we combined a simplified implementation of Smith and Schaefer [2015] contact penalty with our RSD energy to obtain a globally injective mesh parametrisation. We used the triangle inequality enhancement of Su et al. [2020] but did not use a surrounding shell mesh. In some hard cases, overwinding can still occur around interior vertices [Weber and Zorin 2014]. Furthermore, such boundary energies do not work when the initialisation has a self-intersecting boundary. However, this is a limitation of Smith and Schaefer [2015], and we leave a more thorough investigation to future work.

5.6 3D Mesh Deformation

5.6.1 Comparison with ABCD.

We compare with Adaptive Block Coordinate Descent (ABCD) [Naitsat et al. 2020] and FFM in 3D over some challenging test cases. ABCD follows a two-stage procedure in the presence of invalid elements, i.e., map repair and then optimisation. For repair, it uses a flip penalty distinct from the energy in the rest of the domain. In contrast, we minimise the same energy over all elements.

In Figure 15 we compare against ABCD with modified Symmetric Dirichlet. For a twisting bar, ABCD and our framework create inversion-free maps with almost identical isometric quality; the histograms almost overlap. However, our RSD energy creates elements with better Jacobian quality (\(J_{\mathrm{min}}\) in Figure 16). We traced this difference to ABCD’s termination criteria.

For these stress tests, our framework was generally faster than ABCD’s combined map repair and optimisation. With the exception of one example (Figure 16) where our framework was marginally slower (15% for wrench example) our framework converged before ABCD’s map repair completed. For a severely twisted Armadillo, ABCD (PN) failed under fixed-boundary, while our framework minimised both RSD and RAMIPS. As discussed in (Figure 11) however, our framework was statistically over 10\(\times\) faster than ABCD over a large 3D dataset. We also ran FFM for these 3D cases, and while it succeeded for fixed-boundary problems, it failed for the free-boundary Armadillo and wrench.

Fig. 11.

Fig. 12.

Fig. 13.

Fig. 14.

5.7 Pathological Test Cases

The 3D benchmark in (Figure 11) covers a wide range of practical problems, and our framework passes 100%. Nevertheless, it is possible to handcraft pathological tests that are more difficult to optimise. Following on from Figure 15, we optimise a severely twisted bar (Figure 17) by fixing its ends and using FFM, ABCD, and our RAMIPS and RSD energies. This is a corner case where both our RAMIPS energy and FFM generate inversion-free maps, but the final mesh overwinds around itself. ABCD generates a superior map, but with multiple twists and our RSD unfolds the bar to the least feasible number of twists with the same termination criterion. This difference is due to ABCD’s block-wise repair. Our framework instead needs to “blow” the whole mesh in order to un-flip hence, resulting in a global minimum. However, the difference between RAMIPS and RSD, lies in the underlying energy rather than our regularisation scheme.

Fig. 15.

Fig. 16.

Fig. 17.

We further stress test our framework with two randomised ogre heads that are then shrunk by 100\(\times\). In Figure 18, the second initialisation is challenging because intermediate configurations develop around two lines (fixed vertices around ears) and the optimiser has to step through low curvature regions. Neither FFM nor TLC are suitable for such cases, as they perform best under fixed boundary simulations. With our regularisation scheme we were able to minimise RAMIPS/RSD/RSARAP and recover the rest shape, albeit under a much tighter \(10^{-6}\) residual tolerance.

Fig. 18.

5.8 Comparison with Invertible Elasticity

Models that intentionally remove barriers in an attempt to extend energies to the infeasible zone are known as invertible [Irving et al. 2004; Smith et al. 2018; Stomakhin et al. 2012]. Irving et al. [2004] proposed an invertible model at the gradient level (PK1), resulting in an inconsistent treatment of gradients. Energetic consistency was restored by Stomakhin et al. [2012] via \(C^2\) extension, while a state-of-the-art Stable Neo-Hookean (SNH) energy was proposed by Smith et al. [2018] which removed barriers altogether.

Our regularisation also extends distortion energies into the infeasible zone at the singular-value level, so how is it related? Invertible elasticity removes barriers, making it unsuitable for generating inversion-free mappings in the presence of positional constraints. Our regularisation includes the barrier for repaired elements, and only regularises degenerate elements. Once the map is repaired, our approach reduces to those of locally injective methods.

Invertible elasticity is not a suitable substitute for locally injective methods. We deform an elephant in Figure 19 by pushing the boundary vertices to extremes, and compare our RAMIPS energy with SNH [Smith et al. 2018]. While SNH is suitable for invertible simulations, it fails to generate an inversion-free mapping. With a traditional line search, the energy is minimised at a non-injective state, and with a Smith and Schaefer [2015] line search, the simulation stalls.

Fig. 19.

6 Discussion and Conclusion

We have presented a regularisation approach to geometric optimisation of locally non-injective mappings using popular flip-preventing energy. Our regularisation scheme is based on spectral shift in singular values, results in consistent formulation, and preserves convexity when applied on distortion energies (particularly, in 2D).

In addition, we have discussed the analytic eigensystems of the most important energies, such as conformal and isometric distortions. when expressed in terms of regularised stretch invariants. Through numerous examples in 2D and 3D, we have shown that our approach is extremely robust, and its second-order nature allows it to outperform state-of-the-art algorithms. For 2D surface parametrisation problems, our method appears to recover/optimise maps quicker than other methods we tested, while maintaining similar or superior quality. Given that, our framework can be tuned like FFM and ABCD its success rate is no less than these frameworks however, spectral shifting naturally also yields more stable Hessian eigensystems than Quasi-Newton methods which in turn eliminates numerical and divergence issues.

Our framework works well in practice, but has limitations. We do not guarantee bounded distortions. Further, there is not a theoretical guarantee that the global minimum of energy will always be an injective map. As discussed in Section 5.7, it can fail on pathological tests depending on the energy. Regularisation counters the fundamental assumptions of elasticity such as growth condition. The complexity of this issue was explored in Garanzha et al. [2021a] as FFM shares the same limitation. Moreover, most barrier energies where singular value shifting applies are also non-convex. Nevertheless, the shift is a linear transform that does not deter convexity or nonlinearity.

Since our scheme generalises over the regularisation parameter, a future direction is to investigate better estimation of the regularisation parameter for other numerical optimisation problems.

Footnotes

Unless otherwise stated, lower case bold letters \(\boldsymbol {a}\) denote vectors, upper case bold letters \(\boldsymbol {A}\) are second order tensors (matrices), and blackboard bold letters \(\mathbb {A}\) are fourth order tensors (matrix of matrices).

Metric consistency demands that \(\boldsymbol {F}\) and any subsequent derived quantity should also be regularised, not just J.

A function \(\mathcal {D}(\boldsymbol {F})\) is said to be polyconvex if there exists a function \(\bar{\mathcal {D}}(\boldsymbol {F}, \boldsymbol {H}, J) = \mathcal {D}(\boldsymbol {F})\) such that \(\bar{\mathcal {D}}\) is convex with respect to \(\boldsymbol {F}\), \(\boldsymbol {H,}\) and J [Ball 1976]. Many non-convex problems are solvable, but if they are also non-polyconvex they exhibit numerical artefacts such as thinning and unnatural deformations.

Supplementary Material

ZIP File (tog-22-0045-file003.zip)

Download
10.31 MB

A Why the Cofactor Invariants?

Invariants that can be expressed in terms other invariants of the same family are common (e.g. \(II^\star _C = 1/2(I^2_C -II_C)\) with Cauchy-Green) and can be useful when they have clean physical interpretations. In terms of our extended invariants, there are multiple reasons to include the two additional cofactor invariants, \(I_H\) and \(II_H\):

•

A complete set of invariants should complete the characteristic polynomial of \(\boldsymbol {F}\) (note \(I_H\)): \(\sigma _i^3 - I_F\sigma _i^2 + I_H\sigma _i - J = 0\)

•

\(II_H\) corresponds to the Cauchy-Green invariant \(II_C\), which does not appear in the Smith et al. [2019] invariants

•

They appear in well-known energies such as Symmetric Dirichlet, Symmetric ARAP (see Table 1)

•

They appear in our Β parameter (see Algorithm 1)

•

They enable closed-form study of polyconvexity; our supplement proves lack of polyconvexity for Symmetric Dirichlet

•

They significantly simplify the algebra for gradient and Hessian expressions which become forbiddingly long when expressed in terms of singular-values or Smith et al. [2019].

B Gradient and Hessian Tensors

B.1 Gradient and Hessian of Distortion Energies

We can explicitly write generic expressions for the gradient (Equation (3)) and Hessian (Equation (4)) of distortion energies in terms of the extended stretch invariants (Equation (10)) without relying on analytic eigensystems. This allows us to study convexity, and perform convex-concave decompositions in closed-form. To obtain these, we first need the gradients of the invariants

\begin{align*} \boldsymbol {G}_{I_F} & = \frac{\partial I_F}{\partial \boldsymbol {F}} = \boldsymbol {R}, & \boldsymbol {G}_{II_F} & = \frac{\partial II_F}{\partial \boldsymbol {F}} = 2\boldsymbol {F}, & \boldsymbol {G}_{J} & = \frac{\partial J}{\partial \boldsymbol {F}} = \boldsymbol {H} \\ \boldsymbol {G}_{I_H} & = \frac{\partial I_H}{\partial \boldsymbol {F}} = I_F\boldsymbol {R} - \boldsymbol {G}, & \boldsymbol {G}_{II_H} & = \frac{\partial II_H}{\partial \boldsymbol {F}} = 2\boldsymbol {H} \times \boldsymbol {F} & \end{align*}

as well as the Hessians

where \(\mathbb {T}_i\)s are the twist tensor [Smith et al. 2019], \(\boldsymbol {H} \times \boldsymbol {F}\) is the tensor cross product between two \(3\times 3\) tensors, and \(\mathbb {I} \times \boldsymbol {F}\) and \(\boldsymbol {F} \times \mathbb {I} \times \boldsymbol {F}\) are higher-order tensor cross products. Eigen code is given in the supplement. The gradient (Equation (3)) of a distortion energy can be obtained using the chain rule

\begin{align*} \boldsymbol {P} & = \frac{\partial \mathcal {D}}{\partial \boldsymbol {F}} = \boldsymbol {\nabla _F} \hat{\mathcal {D}} (I_F, II_F, I_H, II_H, J) \\ & =\frac{ \partial \hat{\mathcal {D}} }{\partial I_F} \boldsymbol {R} +2\frac{ \partial \hat{\mathcal {D}} }{\partial II_F} \boldsymbol {F} +\frac{ \partial \hat{\mathcal {D}} }{\partial I_H} (I_F \boldsymbol {R} - \boldsymbol {F}) +2\frac{ \partial \hat{\mathcal {D}} }{\partial II_H} \boldsymbol {H} \times \boldsymbol {F} +\frac{ \partial \hat{\mathcal {D}} }{\partial J} \boldsymbol {H}. \end{align*}

This is a \(d \times d\) matrix and can be flattened in to a \(d^2 \times 1\) vector. The Hessian can be expressed as

(15)

\(color{blue}{Blue}\) terms are diagonal entries of the Hessian, and are second derivatives of the energy with respect to the same invariant. \(color{green}{Green}\) terms are off-diagonal entries that are second derivatives of the energy with respect to different invariants. \(color{red}{Red}\) terms are first derivatives of the energy with respect to the invariants, and form the geometric stiffness term [Bonet et al. 2016b]. The “sym” superscript indicates off-diagonal blocks with symmetric counterparts s.t. for every product \(\boldsymbol {A} \odot \boldsymbol {B}\), there exists \(\boldsymbol {B} \odot \boldsymbol {A}\), wherein \(\odot\) is any one of the tensor products along the off-diagonal. Equation (15) can be flattened in to a \(d^2 \times d^2\) matrix.

For polyconvex energies, all terms except the last three \(color{red}{red}\) ones are SPD (the first two, \(\mathbb {I}\) and \(\mathbb {H}_{I_F}\) are SPD), so a convex-concave decomposition can be performed by discarding the last three red terms. Modified-FFM was implemented using this approach.

B.2 Eigenmatrices of the Invariants

Eigenmatrices of stretch invariants have been covered in detail [Kim and Eberle 2022; Kim et al. 2019; Smith et al. 2019], and three types arise: twist (\(\mathbb {T}\)), flip (\(\mathbb {L}\)), and scaling (\(\boldsymbol {D}_i\)) tensors.

B.2.1 General Isotropic Energies.

In 2D, to obtain the analytic eigensystem of an arbitrary isotropic energy, we have

\begin{align} \mathbb {H} = \lambda _1 \boldsymbol {E}_1 \otimes \boldsymbol {E}_1 + \lambda _2 \boldsymbol {E}_2 \otimes \boldsymbol {E}_2 + \lambda _3 \mathbb {L} + \lambda _4 \mathbb {T} \end{align}

(16)

where \(\mathbb {T}\) and \(\mathbb {L}\) are twist and flip tensors and \(\boldsymbol {E}_1\) and \(\boldsymbol {E}_2\) are coupled scaling modes. These can be expressed using decoupled scaling modes \(\boldsymbol {D}_1\) and \(\boldsymbol {D}_2\) as

\begin{align} \boldsymbol {E}_1 = \frac{1}{\sqrt {1 + \tau ^2}} (\tau \boldsymbol {D}_1 + \boldsymbol {D}_2), \quad \boldsymbol {E}_2 = \frac{1}{\sqrt {1 + \tau ^2}} (\boldsymbol {D}_1 - \tau \boldsymbol {D}_2) \end{align}

(17)

where \(\tau = 0\) if the scaling modes decouple. In 3D, the analytic eigensystem can be written as

(20)

Here, \(\boldsymbol {E}_i\) can be computed by solving the \(3\times 3\) eigenvalue problem for the coupling matrix \(\boldsymbol {A}\)

\begin{align} \boldsymbol {E}_i = {\boldsymbol {U}} \mathrm{diag}(\boldsymbol {Q}_i) {\boldsymbol {V}}^T \end{align}

(19)

where \(\boldsymbol {Q}_i\) are the eigenvectors of \(\boldsymbol {A}\) (see Smith et al. [2018, 2019]). We can now present the analytic eigensystem of our two new cofactor invariants, \(I_H\) and \(II_H\).

B.2.2 Analytic Eigensystem of I_H.

For the twist modes we have

\begin{align*} \lambda _{i+3} = \frac{I_F + \sigma _i}{\sigma _j + \sigma _k} \qquad \boldsymbol {E}_{i+3} = \boldsymbol {T}_i \qquad i,j,k = 1,2,3, j\ne i, k \ne i,j \end{align*}

whereas flip modes have constant eigenvalues

\begin{align*} \lambda _{i+6} = -1 \qquad \boldsymbol {E}_{i+6} = \boldsymbol {L}_i \qquad i = 1,2,3 \end{align*}

and the coupled scaling modes

\begin{align*} \lambda _{1} = 2, \quad \lambda _{2,3} = -1 \qquad \boldsymbol {E}_{i} = {\boldsymbol {U}} \mathrm{diag}(\boldsymbol {Q}_i) {\boldsymbol {V}}^T \qquad i = 1,2,3 \end{align*}

where \(\boldsymbol {Q}_{i}\) are the eigenvectors of \(\boldsymbol {A} = [0, 1, 1; 1, 0, 1; 1, 1, 0]\).

B.2.3 Analytic Eigensystem of II_H.

For the twist modes we have

\begin{align*} \lambda _{i+3} =2 \big(\sigma _j^2 + \sigma _i \sigma _k\big) \qquad\kern-3pt \boldsymbol {E}_{i+3} = \boldsymbol {T}_i \qquad\kern-2pt i,j,k = 1,2,3, j\ne i, k \ne i,j \end{align*}

and the flip modes have the same structure

\begin{align*} \lambda _{i+6} =2 \big(\sigma _j^2 - \sigma _i \sigma _k\big) \qquad\kern-2.5pt \boldsymbol {E}_{i+6} = \boldsymbol {L}_i \qquad\kern-2.5pt i,j,k = 1,2,3, j\ne i, k \ne i,j. \end{align*}

The coupled scaling modes are from the eigensystem of:

\begin{align*} \boldsymbol {A} = 2 \begin{bmatrix} \sigma _2^2 + \sigma _3^2 & \quad 2\sigma _1\sigma _2 &\quad 2\sigma _1\sigma _3 \\ &\quad \sigma _1^2 + \sigma _3^2 &\quad 2\sigma _2\sigma _3 \\ sym &\quad &\quad \sigma _1^2 + \sigma _2^2 \end{bmatrix} \end{align*}

where we then apply Equation (19).

References

[1]

Noam Aigerman and Yaron Lipman. 2013. Injective and bounded distortion mappings in 3D. ACM Trans. Graph. 32, 4, Article 106 (July 2013), 14 pages.

Abstract

1 Introduction

2 Related Work

2.1 Problem Statement

2.2 Inversion-free Mappings

2.3 Energy Minimisation

3 Our Regularisation Approach

3.1 Continuum Mechanics Preliminaries

3.2 Regularisation of F

3.3 Regularisation of Stretch Invariants

3.4 Finding the Regulariser, Β

3.4.1 State-of-the-art Approaches.

3.4.2 Our Approach.

3.5 Regularised Energies

3.6 Discussion of Invariant-based Formulation

4 Optimisation Approach

4.0.1 Non-injective States.

5 Numerical Experiments

5.1 Surface Parametrisation

5.1.1 Comparing with FFM.

5.1.2 Comparison with TLC.

5.1.3 Comparison with SA.

5.2 Locally Injective Mappings Benchmark

5.3 Volumetric Polycube Parametrisation

5.4 Free Boundary Surface Parametrisation

5.5 Globally Injective Surface Parametrisation

5.6 3D Mesh Deformation

5.6.1 Comparison with ABCD.

5.7 Pathological Test Cases

5.8 Comparison with Invertible Elasticity

6 Discussion and Conclusion

Footnotes

Supplementary Material

A Why the Cofactor Invariants?

B Gradient and Hessian Tensors

B.1 Gradient and Hessian of Distortion Energies

B.2 Eigenmatrices of the Invariants

B.2.1 General Isotropic Energies.

B.2.2 Analytic Eigensystem of IH.

B.2.3 Analytic Eigensystem of IIH.

References

Cited By

Index Terms

Recommendations

Moreau---Yosida regularization in shape optimization with geometric constraints

Guaranteed Computation of Geometric Mappings Via Convex Optimization

Accelerated quadratic proxy for geometric optimization

Comments

Information

Published In

Publisher

Publication History

Permissions

Check for updates

Author Tags

Qualifiers

Contributors

Other Metrics

Bibliometrics

Article Metrics

Other Metrics

Citations

Cited By

View options

PDF

eReader

HTML Format

Login options

Full Access

Figures

Other

Share

Share this Publication link

Share on social media

Affiliations

B.2.2 Analytic Eigensystem of I_H.

B.2.3 Analytic Eigensystem of II_H.