 # Foldover-free maps in 50 lines of code

Mapping a triangulated surface to 2D space (or a tetrahedral mesh to 3D space) is the most fundamental problem in geometry processing.In computational physics, untangling plays an important role in mesh generation: it takes a mesh as an input, and moves the vertices to get rid of foldovers.In fact, mesh untangling can be considered as a special case of mapping where the geometry of the object is to be defined in the map space and the geometric domain is not explicit, supposing that each element is regular.In this paper, we propose a mapping method inspired by the untangling problem and compare its performance to the state of the art.The main advantage of our method is that the untangling aims at producing locally injective maps, which is the major challenge of mapping.In practice, our method produces locally injective maps in very difficult settings, and with less distortion than the previous work, both in 2D and 3D. We demonstrate it on a large reference database as well as on more difficult stress tests.For a better reproducibility, we publish the code in Python for a basic evaluation, and in C++ for more advanced applications.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

Most geometric objects are represented by a triangulated surface or a tetrahedral mesh. The mapping problem consists in generating a 2D or 3D map of these objects. This is a fundamental problem of computer graphics because it is much easier for many applications to work in this map space than to directly manipulate the object itself. To give few examples, texture mapping stores colors of a surface as images in the map space, remeshing uses global maps in 2D (Bommes et al., 2013) and 3D (Gregson et al., 2011; Nieser et al., 2011). In addition, mapping algorithms can be used to deform volumes (Li et al., 2020), or generate shells from surfaces (Jiang et al., 2020).

#### What is a good map?

Most often maps are represented by the position of the vertices in the map space, and interpolated linearly on each element (triangle or tetrahedron). In a perfect world, the map space would keep the geodesic distances of the object. Unfortunately, this is usually impossible due to Gaussian curvature, and application-specific constraints such as constrained position of vertices or overlaps in the map space. Therefore, the objective of the mapping algorithms is to minimize the distortion between geometric and map spaces, opening the door for numerical optimization approaches as detailed in the survey

(Hormann et al., 2008).

Unfortunately, when high distortion is required to satisfy the constraints, these algorithms may lose the fundamental property of a maps: injectivity. A solution to preserve it (Floater, 1997) relies on Tutte’s theorem (Tutte, 1963), however the surface boundary must be mapped to a convex polygon. Despite this strong limitation, it still remains the reference algorithm to generate injective maps. Lower distortion can be obtained by changing weights of the barycentric coordinates (Eck et al., 1995) (as long as they are not negative), and alternative solutions (Campen et al., 2016; Shen et al., 2019) have been explored to improve robustness to numerical imprecision by modifying the mesh connectivity.

#### Local invertibility.

In many applications, maps are used to access a neighborhood of a point within a coherent local coordinate system. To this end, global injectivity is not required, and we instead look for local injectivity (Schüller et al., 2013; Smith and Schaefer, 2015). Their approach starts from an injective map (Floater, 1997), and maintains the local injectivity when minimizing the distortion and enforcing the constraints. This allows them to optimize at the same time the parameterization and the texture packing (Jiang et al., 2017), with a possibility to scale to larger meshes (Rabinovich et al., 2017).

#### Recover local injectivity.

Local injectivity can also be recovered for a map with few foldovers present. For example, in 2D (Lipman, 2012) and 3D (Aigerman and Lipman, 2013), the map is projected on a class of bounded distortion maps. The numerical method is however unlikely to succeed for stiff problems. Recovering local injectivity is also known as mesh untangling. Originally related to Arbitrary Lagrangian-Eulerian moving mesh approach, the mesh untangling problem considers a simplicial complex with misoriented elements and attempts to flip them back by optimizing the position of the vertices. There is an abundant literature on mesh untangling (Du et al., 2020; Knupp, 2001; Freitag and Plassmann, 2000; Escobar et al., 2003; Toulorge et al., 2013), however the common opinion is that untangling is a very hard problem and algorithms are not robust enough. As a manifestation of frustration over this problem (Danczyk and Suresh, 2013) investigates a finite element method working directly on tangled (sic!) meshes.

#### Elastic deformations.

To recover local injectivity, we propose a method stemming from the computational physics. It is very important to note that there is rich literature on mesh deformation in the community working on grid generation for scientific computation. Numerical simulation of hydrodynamic instability of layered structures requires sound mathematical foundations behind moving deforming mesh algorithms. In the ’60s Winslow and Crowley, independently one from another, introduced mesh generation methods based on inverse harmonic maps (Crowley, 1962; Winslow, 1966).

Since then, a lot of effort was spent on mesh generation based on elastic deformations (Jacquotte, 1988), but mostly for regular grids. In 1988, at the time of domination of finite difference mapped grid generation methods, S. Ivanenko introduced the pioneering concept of barrier variational grid generations methods guaranteeing construction of non-degenerate grids (Ivanenko, 1988; Charakhch’yan and Ivanenko, 1997). To generate deformations with bounded global distortion (bounded quasi-isometry constant), Garanzha proposed to minimize an elastic energy for a hyperelastic material with stiffening suppressing singular deformations (Garanzha, 2000). Invertibility theorem for deformation of this material was established in the 3D case as well (Garanzha et al., 2014).

A solid mathematical ground for these methods was laid by J. Ball who introduced in 1976 his theory of finite elasticity based on the concept of polyconvex distortion energies (Ball, 1976). He not only proved Weierstrass-style existence theorem for this class of variational problems, but also formulated a theorem on invertibility of elastic deformations for quite general 3D domains (Ball, 1981). It is important to note that Ball invertibility theorem is proved for Sobolev mappings and can be applied directly for finite element spaces, i.e. to deformation of meshes, as was pointed out in (Rumpf, 1996).

#### Our contributions

Inspired by these results on untangling and elastic deformations, we propose a method outperforming recent state of the art on locally injective parameterization (Du et al., 2020) in terms of robustness, quality and supported features. To sum up, our contributions are:

• we propose a method that optimizes the map distortion directly, with a parameter for a tradeoff between the angle and area preservation,

• we provide a rigourous analysis of the foundations of our numerical resolution scheme,

• our method supports mapping with free boundaries,

• we observe better robustness for high distortions as well as with respect to a poor quality initialization.

• To ease the reproducibility, we publish the code in Python (refer to Listing 1) for a basic evaluation, and a C++ code in the supplemental material for more advanced applications.

The rest of the paper is organized as follows: we start with presenting our method in § 2, then we evaluate its performance (§ 3.1 and § 3.2) as well as its limitations (§ 3.3). Then we present theoretical guarantees for our resolution scheme: in § 4.1 we prove that our approximation of Hessian matrix is definite positive, and finally in § 4.2 we prove that our choice of the regularization parameter sequence guarantees that a minimization algorithm111The minimization algorithm is subject to conditions of Th. (1); from our numerical experiments we observe that our minimization algorithm almost always respects the conditions. can find a mesh free of inverted elements in a finite number of steps.

## 2. Penalty method for mesh untangling

In this rather short section we present our method of computing a foldover-free map , i.e. we map the domain to a parametric domain. This presentation is unified both for 2D and 3D settings, and by

we denote the number of dimensions; in our notations we use arrows for all vectors of dimension

.

The section is organized as follows: in § 2.1 we give a primer on the variational formulation of mapping problem in continuous settings, then we state our problem in § 2.2 as a regularization of this variational formulation, and finally we present our numerical resolution scheme in § 2.3.

### 2.1. Variational formulation for grid generation

Let us denote by a map to a parametric domain: for the flat 2D case we can write as , and for a 3D map .

Consider the following variational problem:

 (1) argmin→u∫Ω(f(J)+λg(J))dx,

where is the Jacobian matrix of the mapping , and

 f(J):=⎧⎪⎨⎪⎩trJ⊤J(detJ)2d,detJ>0+∞,detJ≤0
 g(J):={detJ+1detJ,detJ>0+∞,detJ≤0

Problem (1) may be subject to some constraints that we do not write explicitly. To give an example, one may pin some points in the map. In this formulation, functions and have concurrent goals, one preserves angles and the other preserves the area, and thus serves as a trade-off parameter.

As a side note, with and , Prob. (1) presents a variational formulation of an inverse harmonic map problem. Namely, if we write down the Euler-Lagrange equations for Prob. (1) and interchange the dependent and independent variables222Attention, this step assumes that the solution of Prob. (1) is a diffeomorphism., we obtain the Laplace equation (not to be confused with omnipresent !). For this case, Prob. (1) is often reffered as to Winslow’s functional, however Winslow himself has never formulated the variational problem, working with inverse Laplace equations. To the best of our knowledge, the first publication of the variational problem is made by Brackbill and Saltzman (Brackbill and Saltzman, 1982).

Prob. (1) provides a very powerful tool based on the theory of finite hyperelasticity introduced by J. Ball (Ball, 1976). From a computational point of view, starting with a map without foldovers, this tool allows us to optimize the quality of the map. Note that the energy (1) is a polyconvex function (refer to App. B–Rem. 3 for a proof) satisfying the ellipticity conditions, and therefore is very well suited for a numerical optimization provided that we have an initial guess in the admissible domain .

The problem, however, is that while being theoretically sound, this problem statement does not offer any practical way to get rid of foldovers in a map, because for a map with foldovers the energy is infinite and provides no indications on how to improve the situation, hence we propose to alter a little the problem statement.

### 2.2. Penalty method

We can avoid nonpositive denominators in and using a regularization function for a positive value of (Fig. 2):

 (2) χ(D,ε):=D+√ε2+D22, Figure 2. Regularization function for the denominator in Eq. (3). When ε tends to zero, χ(ε,D) tends to D for positive values of D, and to 0+ for negative values of D.

Then we define a regularized version of functions and :

 (3) fε(J):=trJ⊤J(χ(detJ,ε))2d,gε(J):=det2J+1χ(detJ,ε),

so that Prob (1) is reformulared as

 (4) limε→0+argmin→u∫Ω(fε(J)+λgε(J))dx

Under certain assumptions333The fact that the solution of Prob. (4) is a diffeomorphism is sufficient (but not necessary) for the equivalence. solutions of Prob. (4) are solutions of Prob. (1), however, Prob. (4) does offer a way of getting rid of foldovers if a foldover-free initialization is not available.

In practice, the map is piecewise affine with the Jacobian matrix being piecewise constant, and can be represented by the coordinates of the vertices in the parametric domain . Let us denote the vector of all variables as , then our optimization problem can be discretized as follows:

 (5) limε→0+argminUF(U,ε), where F(U,ε):=#T∑t=1(fε(Jt)+λgε(Jt))vol(Tt),

is the number of vertices, is the number of simplices, is the Jacobian matrix for the simplex and is the volume of the simplex in the original domain.

### 2.3. Resolution scheme

To solve Prob. (5), we use an iterative descent method. We start from an initial guess , and we build a sequence of approximations . For each iteration we need to carefully choose the regularization parameter . Starting from , we define the sequence as follows:

 (6) εk+1:=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩(1−σk√(Dk+1−)2+(εk)2|Dk+1−|+√(Dk+1−)2+(εk)2)εk, if Dk+1−<0,(1−σk)εk, if Dk+1−≥0,

where is the minimum value of the Jacobian determinant over all cells of the mesh at the iteration , and . Note that Eq. (6) is justified by Th. 1 (§ 4.2

) on finite untangling sequence; while this formula provides some guarantees, one may use a simpler heuristic to accelerate the computations, more on this in §

3.2.

The simplest way to find is to call a quasi-Newtonian solver such as L-BFGS-B (Zhu et al., 1997). The only thing we need to implement is the computation of the function and its gradient .

Another option is to compute analytically the Hessian matrix instead of estimating it. The problem, however, is that the Hessian matrix

is not positive definite. In this paper we propose its approximation that ensures the positive definiteness. The modified Hessian matrix of the function with respect to at the point is built out of blocks

 H+ij≈∂2F∂→ui∂→u⊤j(Uk,εk).

Here, the matrix is placed on the intersection of -th block row and -th block column; the symbol means that we remove all the terms depending on the second derivative of and second derivatives of to keep positive definite. Refer to Appendix A for the formulæ, and to § 4.1 for the proof of the positive definiteness of .

A detailed description of the resolution scheme is given in Alg. 1. Refer to List. 1 and Fig. 3 for a complete working example of Python implementation and the corresponding input/output generated by the code. Note that our method is not limited to simplicial meshes only: in this particular example we evaluate the Jacobian matrix for every triangle forming quad corners, what corresponds to the trapezoidal quadrature rule. Figure 3. The input mesh with foldovers and the untangling produced by the Listings 1 and 2. Left: a quad 2D mesh to untangle. The boundary (in red) is locked, and the black mesh is free to move. Right: fold-free result. Figure 4. Injective mapping sanity check: make two points inside a square switch places. Left: the input problem, all locked points are shown in red. Right: foldover-free result obtained with our method. Figure 5. Constrained boundary injective mapping challenge proposed by (Du et al., 2020): the “Lucy” mesh is mapped to a P-shaped domain by constraining the vertices shown in red. Left column: three different initializations for the same problem. Middle column: our method produces the same (up to a numerical precision) result on all three initializations. Right column: total lifted content method (Du et al., 2020) fails to solve for the randomly initialized interior vertices, and produces very different results on other two initializations. Three miniature images of “Lucy” show in red the portion of the surface visible in the corresponding close-ups.

## 3. Results and discussion

In this section we provide an experimental evaluation of the method. In the field of computer graphics, any claim about map injectivity always faces a simple sanity check (Fig. 4): take a square and make two points inside switch places. Our method successfully avoids the desk-reject, so we start this section (§ 3.1) by testing our method on the benchmark (Du et al., 2020), then we continue with further tests we have found relevant (§ 3.2), and finally we discuss the limitations of the approach in § 3.3.

### 3.1. Benchmark database

Along with their paper, Du et al. have published a valuable benchmark database. It contains a huge number of 2D and 3D constrained boundary injective mapping challenges. For 2D challenges, the benchmark contains 3D triangulated surfaces to flatten, and not flat 2D meshes as we have described in § 2.1. Nevertheless, our method can handle it directly because the mapping is still on each triangle. To the best of our knowledge, Total Lifted Content (TLC) (Du et al., 2020) and our method are the only ones passing the benchmark without any fail.

A representative example from the database is given in top row of Fig. 6. The challenge is to map the “Lucy” mesh statuette from the Stanford Computer Graphics Laboratory to a P-shaped domain. This mesh has a topology of a disk, and its boundary vertices are uniformly spaced on the P-shape boundary. As an initialization to the problem, Du et al. have computed the corresponding (flat) minimal surface that obviously contains a foldover (Fig. 6–top left). Then the problem boils down to a mesh untangling with locked boundary.

#### Mapping quality measure

How to measure quality of a map? Well, it depends on the goal. An identity is an unreachable ideal; traditional competing goals are (as much as possible) angle preserving and area preserving maps. Thus, we can measure the extreme values of the failure of a map to be conformal or authalic. Our maps being piecewise affine, the Jacobian matrix

is constant per element. Let us define the largest singular value of

as , and the smallest singular value as ; then the quality of a mapping can be reduced to extreme values of the stretch () and the scaling ().

For the “Lucy-to-P” challenge (Fig. 6–top row) our map differs from the TLC result by 12 orders of magnitude in terms of minimum scaling, and by two orders of magnitude in terms of maximum stretch. To visualize this difference in scaling, we have provided the close-ups: Fig. 6

–top middle shows a map of the Lucy’s torch, whereas the same level of zoom on the result by Du et al. (Fig.

6–top right) contains not only the torch, but also both wings, the head and the right arm!

Note also that the input “Lucy” mesh is slightly anisotropic; our method allows us to prescribe the element target shape, so the dress pleats are clearly visible in our mapping. Figure 7. Quality plot of the resulting locally injective maps for every challenge from the database provided by (Du et al., 2020). Our results are shown in blue, whereas the results by Du et al. are shown in red. Each dot corresponds to a quality of the corresponding map reduced to two numbers: the maximum stretch and the mininum scale. Top row: mapping quality on the 2D dataset. Bottom row: mapping quality on the 3D dataset. Left column shows the absolute maximum stretch and absolute minimum scale, whereas the right column shows the maximum stretch and minimum scale for the top 95% of measurements.

#### Benchmark database

Our method successfully passes all challenges from the benchmark (Du et al., 2020). The benchmark consists of 10743 meshes to untangle in 2D and 904 meshes in 3D under locked boundary constraints. In Fig. 7 we provide quality plots of the resulting locally injective maps. These are scatter plots: each dot corresponds to a quality of the corresponding map reduced to two numbers: the maximum stretch () and the minimum scaling (). Left column of Fig. 7 shows the worst quality measurements for every 2D problem (top) as well as for every 3D challenge (bottom) of the dataset. Our results are shown in blue, whereas TLC results are shown in red. To illustrate the distribution of the elements’ quality, for each injective map we have removed 5% of worst measurements: the right column of Fig. 7 shows the maximum stretch and the minimum scaling for the top 95% of measurements.

Note the dot arrangements forming lines in the plot: these dots correspond to the few sequences of deformation present in the database. Figure 8. Performance of our method tested on the benchmark (Du et al., 2020). Each dot corresponds to a challenge from the database (10743 in 2D and 904 in 3D). Blue dots show the running times obtained using Eq. (7), green dots correspond to Eq. (6).

#### Timings

Fig. 8 provides a scatter plot of our running time vs mesh size for all the challenges from the database (Du et al., 2020): for each run, the time varies from a fraction of a second to several minutes for the largest meshes. These times were obtained with a 12 cores i7-6800K CPU @ 3.40 GHz. As in Fig. 8, the vertical lines in the 3D dataset plot correspond to the sequences of deformation in the benchmark.

There are two scatter plots superposed, both represent the same resolution scheme with an exception corresponding to the way we compute (Alg. 1–line 3). The green scatter plot corresponds to a conservative update rule (Eq. (6)) offering guarantees on untangling in a finite number of steps (refer to Th. 1), whereas the blue scatter plot is obtained using the following update rule:

 (7) εk:=√10−12+4⋅10−2⋅[min(0,mint∈1…#TdetJkt)]2,

This formula was chosen empirically, and it performs very well in the vast majority of situations. In particular, it allows for all the database (Du et al., 2020) to pass the injectivity test.

### 3.2. Further testing

#### Sensitivity to initialization

Our next test is the sensitivity to the initialization. We have generated two other initializations for the “Lucy-to-P” challenge: the one with all interior vertices collapsed onto a single point (Fig. 6–middle row), and with the interior vertices being randomly placed withing a bounding square (Fig. 6–bottom row).

Our method produces virtually the same result on all three initializations, whereas TLC generates very different results for the first two, and fails for the third one. It is interesting to note that TLC is heavily depending on the initialization, it alters very little the input geometry.

#### Large deformation stress test

For our next test we have generated an isotropic tetrahedral mesh of a cube with a cavity, and we rotated the inner boundary to test the robustness of our method to large deformations. Figure 6 shows the results. Our L-BFGS-based optimization scheme succeeds up to the rotation of 135°, and we had to switch to the Newton method to reach the 180° rotation. TLC method had succeded on 45° and 90°, and failed for the 135° and 180°. Note that as in the previous test, even when the untangling succeeds, TLC alters very little the input map, thus producing heavily stretched tetrahedra, whereas our method evenly dissipates the stress over all the domain. Figure 9. Free boundary injective mapping. The vertices shown in red are constrained, all other vertices are free to move. (a): a compression test, (b): a bend test. Refer to Fig 10–a for the rest shape.

#### Free boundary injective mapping

To the best of our knowledge, our method is the only one able to produce inversion-free maps with free boundary. Since TLC tries to minimize the overall volume, relaxing the boundary constraints results in degenerate maps.

Fig. 9 shows two maps obtained with our method: a 2D shape being compressed and the same shape being bent. The boundary is free to move, we lock the vertices shown in red. Refer to Fig. 10-a for the rest shape. The shape behaves exactly as a human would expect it: upon compression the shape chooses one of the two possible results (Fig. 9–a), and successfully passes the bend test (Fig. 9–b), note the geometrical details that are naturally rotated. Figure 10. Free boundary injective mapping: influence of the parameter λ in Prob. (5). The vertices shown in red are constrained, all other vertices are free to move. (a): the rest shape; (b): a stretch preserving the shape of the elements (λ=0); (c): an area-preserving stretch (λ=104); (c): a trade-off between the shape and the area preservation (λ=1).

Our final tests illustrate the influence of the parameter in Prob. (5) on the resulting map. We have computed three free boundary maps of the rest shape (Fig. 10–a) being stretched. First we chose , that is, only the shape quality term is taken into account in Prob. (5). When we optimize for the angles, the area of the triangles is forced to change, refer to Fig. 10–b for the resulting map. Naturally, an area preserving map () must deform the elements to satisfy the area constraint (Fig. 10–c). Finally, in Fig. 10–d we show an example with a tradeoff between the area and angles preservation.

### 3.3. Limitations

While globally performing very well in practice, our method still presents some limitations. We have two main sources of limitations: theoretical limitations as well as very practical ones related to numerical stability of our resolution scheme.

#### Overlaps

First of all, an inversion-free map does not imply global injectivity. Fig. 11–a provides an example of an inversion-free map with two cases of non-injectivity when optimizing for a map with free boundaries: the map can present global overlaps as well as the boundary can “wind up” around boundary vertices, i.e. the total angle of triangles incident to a vertex can be superior to . Moreover, while being less frequent, similar situations may occur on interior vertices. Typically this situation happens near constraints causing a local compression in the shape.

Let us illustrate this behavior on a very simplistic mesh consisting of a single fan of 12 triangles. All vertices are free to move, the target shape is set to be the unit equilateral triangle for all elements. For this problem Fig. 12–a shows a local minimum, and the Fig. 12–c shows the global minimum respecting perfectly the prescribed total angle of around the center vertex. Both are inversion-free maps, but only the map in Fig. 12–a is a globally injective one. Depending on the initialization and the resolution scheme chosen, we can converge to either minimum. Note, however, that the center vertex has the winding number 1 in one map and 2 in the other, and thus we can not deform continuously one to the other without inverting some elements. Note also that the configurations like in the Fig. 12–b present inverted elements and thus can not be generated by our method.

It is possible to avoid all overlaps altogether by embedding our shape to optimize into an outer triangulation, and performing a “bi-material” optimization. In this case, both global overlaps and fold-2-coverings are prohibited by the the fact that the outer material must not have inverted elements (refer to Fig. 11–b). The thick prismatic layer in Fig. 1 was generated by a similar procedure: we have generated a very thin layer of triangular prisms around the dragon, and tetrahedralized the exterior bounded by a cube. After calling the untangling procedure, we have obtained an offset surface with exactly the same mesh connectivity as the original dragon mesh.

While this embedding kind of approach works well for certain applications, for other it may be hard to apply. We are currently exploring simpler practical ways to fix the winding problem. Figure 11. Free boundary mapping limitations. (a): this mesh presents two kinds of problems, namely, a global overlap and the mesh wrapped around a boundary vertex. (b): Both problems can be avoided by embedding the mesh into an outer triangulation. Figure 12. Free boundary mapping limitations: three maps of a very simplistic mesh made of 12 triangles. (a) and (c) both are inversion-free maps and thus allowed by our method, whereas the map (b) has inverted elements, and thus is prohibited by our method.

#### Numerical challenges

Even when the problem is well-posed, a robust implementation may present significant difficulties. As we have said above, the quasi-Newtonian optimization scheme performs well for “simple” problems (it passes all the benchmark database!), but may fail for large deformations, where Newton iterations are necessary. While our modified Hessian matrix is symmetric positive definite, note that for stiff problems the Jacobi preconditioned conjugate gradients can fail and one might need the incomplete Choleski decomposition and beyond.

In practice we have found the method being very robust in 2D settings: we have not encountered a practical test case we were not able to treat with our method. In 3D, however, it can fail due to the numerical challenges in very anisotropic and highly twisted meshes.

## 4. Analysis

This section presents in two parts a rigorous analysis of the penalty method. First in § 4.1 we prove that the modified Hessian matrix is indeed positive definite, and then in § 4.2 we show the origins of Eq. (6) for the regularization parameter sequence . Namely, we prove that if the problem has a solution, then an idealized minimization method can reach the admissible set in a finite number of steps. An immediate consequence of this theorem is that if the problem has a solution, then for some the solution belongs to the admissible set.

### 4.1. Modified Hessian matrix

Recall that in our resolution scheme we use the modified Hessian matrix of the function built out of blocks placed on the intersection of -th block row and -th block column. It is a common practice to add some regularization terms to the Hessian matrix to make it positive definite, but we propose to modify the finite element (FE) matrix assembly procedure by eliminating some terms potentially leading to an indefinite FE matrix.

To this end, we restrict our attention to a single simplex and we study a function of the Jacobian matrix defined as follows:

 (8) ϕ(J):=fε(J)+λgε(J)=trJ⊤J(χ(detJ,ε))2d+λdet2J+1χ(detJ,ε)

Let us denote by the (columnwise) flattening of the Jacobian matrix , i.e. the vector composed of the elements of . We decompose the Hessian matrix of with respect to the Jacobian matrix entries into two parts: , where is a positive definite matrix, and can be an indefinite matrix that we neglect. The matrix contains all terms depending on and second derivatives of with respect to elements of the Jacobian matrix . Our map is affine on the simplex of interest, therefore its Jacobian matrix is a linear function of the vertices of the simplex. Thus, the idea is to compute a positive definite matrix

, and use the chain rule to get the Hessian matrix with respect to our variables

and assemble the matrix .

So, we choose some arbitrary point and we want to show the way to decompose into a sum of and with . To do so, first we write down the first order Taylor expansion of the function around some point :

 q(D):=χ(D0,ε)+∂χ∂D(D0,ε)(D−D0).

Next we define a function as follows:

 Φ(a,D):=|a|2(q(D))2d+λD2+1q(D).

Note that differs a bit from : it has one more argument and is replaced by its linearization in the denominator. While this maneuver might seem obscure, light will be shed very shortly. has a major virtue of being convex! The convexity is easy to prove, refer to Appendix B for a formal proof.

Having built a convex function , it is straightforward to verify that the following decomposition holds:

 (9) ∂2ϕ∂a∂a⊤(J0)=M+(J0)+M±(J0),

where

 M+ :=(I∂D∂a)⎛⎜⎝∂2Φ∂a∂a⊤∂2Φ∂a∂D∂2Φ∂D∂a⊤∂2Φ∂D2⎞⎟⎠(I∂D∂a)⊤, and M± :=∂Φ∂D∂2D∂a∂a⊤−χ′′χ(2dfε+λgε)∂D∂a∂D∂a⊤.

The easiest way to check that the equality (9) holds is to note that at the point we have , , and therefore we have

 ∂ϕ(a)∂a=∂Φ(a,D(a))∂a+∂Φ(a,D(a))∂D∂D(a)∂a.

To calculate the Hessian , it suffices to differentiate this expression one more time and add the terms in that were zeroed out by the linearization.

To sum up, in our computations, for each simplex we approximate the Hessian matrix by the matrix and we neglect the term . Thanks to the convexity of , it is trivial to verify that for any choice of the matrix is definite positive. Then we use the chain rule over to get the Hessian matrix with respect to our variables , and we assemble a approximation of the Hessian matrix for the energy function . Matrix is positive definite provided that at least mesh vertices are fixed. If less than points are fixed, rigid body transformations are allowed. The energy is invariant w.r.t rigid body transformations, so when constraints allow for such transformations, matrix becomes positive semidefinite Note that the leading blocks are always positive definite. Refer to Appendix A for further details on the finite element assembly procedure.

### 4.2. Choice of εk

In this section we provide a strategy for the choice of the regularization parameter at each iteration. Namely, we prove that an idealized minimization algorithm reaches the admissible set in a finite number of iterations.

###### Theorem 1 ().

Let us suppose that the admissible set is not empty, namely there exists a mesh satisfying . We also suppose that we have a minimization algorithm satisfying one of the following efficiency conditions for some :

• Either the essential descent condition holds

 (10) F(Uk+1,εk)≤(1−σ)F(Uk,εk),
• or the vector satisfies the quasi-minimality condition:

 (11) minUF(U,εk)>(1−σ)F(Uk,εk).

Then the admissible set is reachable by solving a finite number of minimization problems in with fixed for each problem. In other words, under a proper choice of the regularization parameter sequence , we obtain .

###### Proof.

The main idea is to expose an explicit way to build a decreasing sequence such that the sequence is bounded from above. Then we can prove that the admissible set is reachable in a finite number of steps by a simple reductio ad absurdum argument.

First of all, the function can be rewritten as follows

 (12) F(U,ε)=∑iαiψi(U)χ(Di,ε),

where denotes the Jacobian determinant for -th simplex of the mesh (), and are positive, separated from zero weights assigned to each simplex. The functions

 ψi(U,ε):=χ(Di,ε)1−2dtrJ⊤iJi+λ(D2i+1)

defined according to Eq.(3) are positive and bounded from below as

 αiψi(U,ε)≥λminiαi.

Note also that are increasing functions of .

Our goal is to build a decreasing sequence such that the sequence is bounded from above. We split the construction into two parts: first we suppose that at some iteration the essential condition (10) is satisfied, and then we explore the case (11).

Suppose that the condition (10) holds at iteration . In order to guarantee that the function does not increase, it suffices to establish the following inequality:

 (13) (1−σ)F(Uk+1,εk+1)≤F(Uk+1,εk).

By noting that , Ineq. (13) is implied if the following condition holds:

 (14) ∀i:(1−σ)χ(Dk+1i,εk)≤χ(Dk+1i,εk+1)

where denotes the Jacobian determinant of simplex at iteration .

Let us show a constructive way to build such that Ineq. (14) is satisfied. To do so, first note that is convex with respect to the second argument, hence

 χ(Dk+1i,εk+1)−χ(Dk+1i,εk)≥(εk+1−εk)∂∂εχ(Dk+1i,εk).

Thus, since and are both positive, Ineq. (14) holds provided that

 (15) ∀i:εk+1≥εk−σχ(Dk+1i,εk)∂∂εχ(Dk+1i,εk)

In its turn, since is an increasing function of , condition (15) can be simplified to

 εk+1≥εk−σχ(Dk+1−,εk)∂∂εχ(Dk+1−,εk),

where is the minimum value of the Jacobian determinant over all simplices.

Hence, if is an approximate solution of the minimization problem with fixed parameter , we can use the following update rule for :

 (16) εk+1=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩εk−σχ(Dk+1−,εk)∂∂εχ(Dk+1−,εk), if Dk+1−<0,εk−σχ(0,εk)∂∂εχ(0,εk), if Dk+1−≥0,

This update rule guarantees that Ineq. (13) is satisfied; coupled with the assumption (10) of the theorem, this implies the required non-growth property of the function values sequence:

 (17) F(Uk+1,εk+1)≤F(Uk,εk).

Consider now the case where condition (11) holds at iteration . Note that condition (11) essentially means that our current solution is very close to the global minimum of , and thus Ineq. (10) cannot be satisfied. Nevertheless, we can use the same update rule (16) for computation of . Indeed, with this choice we have

 F(Uk+1,εk+1) ≤1(1−σ)F(Uk+1,εk)≤1(1−σ)F(Uk,εk)< <1(1−σ)2minUF(U,εk)<1(1−σ)2minUF(U,0).

Here the last inequality provides a global bound on the function values sequence, and it is based on the observation .

To sum up, we have shown a way to build a sequence such that the sequence is bounded from above. Now let us prove that the update rule (16) allows to reach the admissible set in a finite number of steps. To do so, we use a reductio ad absurdum argument.

Suppose that the admissible set is never reached for an infinite decreasing sequence built using the update rule (16), i.e. .

Then, using the facts that and form a decreasing sequence of non-negative values, we can rewrite the update rule (16) for some :

 ε0−εK≥2σK−1∑k=0χ(Dk+1−,εk)≥2σKmin0≤k

with an immediate consequence that for an arbitrarily large we have the following inequality:

 max0≤k

Since all terms in (12) are bounded from below, the resulting estimate contradicts the boundedness of , thus concluding our proof.

###### Remark 1 ().

An important corollary of Th. 1 is that, provided that the admissible set is not empty, there exists an iteration such that the global minimum of the function belongs to the admissible set. The proof is rather obvious: suppose we have an idealized minimizer such that . This minimizer always satisfies the conditions of Th. 1, therefore it can untangle the mesh in a finite number of steps.

###### Remark 2 ().

For our choice of the regularization function , the update rule (16) can be instatiated as follows:

 (18) εk+1=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩(1−σ√(Dk+1−)2+(εk)2|Dk+1−|+√(Dk+1−)2+(εk)2)εk, if Dk+1−<0,(1−σ)εk, if Dk+1−≥0,

In practice the global estimate is not known in advance, and the optimization routine may be far from the ideal. For each minimization step we compute the local descent coefficient :

 σk:=1−F(Uk+1,εk)F(Uk,εk).

When one can use the update rule (18) using the local value guaranteeing that Ineq. (17) holds. In the case one should check that condition (11) holds for prescribed . If positive, we can assign and use update rule (18). If one cannot assure (11), it means that minimization procedure for failed and theorem cannot be applied (it does not mean that Alg. 1 will not reach the admissible set!). In numerical experiments we use value .

## 5. Conclusion

Producing maps without reverted elements is a major challenge in geometry processing. Inspired by untangling solutions in computational physics, our solution outperforms the state of the art both in terms of robustness and distortion. It is easy to use since we provide a simple implementation that is free of commercial product dependency (compiler, library, etc.). Moreover, the energy is estimated independently on each triangle / tetrahedra, making it a good candidate to be adapted to more difficult settings including free boundary and global parameterization.

## References

• (1)
• Aigerman and Lipman (2013) Noam Aigerman and Yaron Lipman. 2013. Injective and Bounded Distortion Mappings in 3D. ACM Trans. Graph. 32, 4, Article 106 (July 2013), 14 pages.
• Ball (1976) John M Ball. 1976. Convexity conditions and existence theorems in nonlinear elasticity. Archive for rational mechanics and Analysis 63, 4 (1976), 337–403.
• Ball (1981) J. M. Ball. 1981. Global invertibility of Sobolev functions and the interpenetration of matter. Proceedings of the Royal Society of Edinburgh: Section A Mathematics 88, 3-4 (1981), 315–328.
• Bommes et al. (2013) David Bommes, Marcel Campen, Hans-Christian Ebke, Pierre Alliez, and Leif Kobbelt. 2013. Integer-Grid Maps for Reliable Quad Meshing. ACM Trans. Graph. 32, 4, Article 98 (July 2013), 12 pages.
• Brackbill and Saltzman (1982) J.U Brackbill and J.S Saltzman. 1982. Adaptive zoning for singular problems in two dimensions. J. Comput. Phys. 46, 3 (1982), 342 – 368.
• Campen et al. (2016) Marcel Campen, Cláudio T. Silva, and Denis Zorin. 2016. Bijective Maps from Simplicial Foliations. ACM Trans. Graph. 35, 4, Article 74 (July 2016), 15 pages.
• Charakhch’yan and Ivanenko (1997) AA Charakhch’yan and SA Ivanenko. 1997. A variational form of the Winslow grid generator. J. Comput. Phys. 136, 2 (1997), 385–398.
• Crowley (1962) WP Crowley. 1962. An equipotential zoner on a quadrilateral mesh. Memo, Lawrence Livermore National Lab 5 (1962).
• Danczyk and Suresh (2013) Josh Danczyk and Krishnan Suresh. 2013. Finite element analysis over tangled simplicial meshes: Theory and implementation. Finite Elements in Analysis and Design 70-71 (2013), 57 – 67.
• Du et al. (2020) Xingyi Du, Noam Aigerman, Qingnan Zhou, Shahar Z. Kovalsky, Yajie Yan, Danny M. Kaufman, and Tao Ju. 2020. Lifting Simplices to Find Injectivity. ACM Trans. Graph. 39, 4, Article 120 (July 2020), 17 pages.
• Eck et al. (1995) Matthias Eck, Tony DeRose, Tom Duchamp, Hugues Hoppe, Michael Lounsbery, and Werner Stuetzle. 1995. Multiresolution Analysis of Arbitrary Meshes. In Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Techniques (SIGGRAPH ’95). Association for Computing Machinery, New York, NY, USA, 173–182.
• Escobar et al. (2003) José Marıa Escobar, Eduardo Rodrıguez, Rafael Montenegro, Gustavo Montero, and José Marıa González-Yuste. 2003. Simultaneous untangling and smoothing of tetrahedral meshes. Computer Methods in Applied Mechanics and Engineering 192, 25 (2003), 2775–2787.
• Floater (1997) Michael S. Floater. 1997. Parametrization and Smooth Approximation of Surface Triangulations. Comput. Aided Geom. Des. 14, 3 (April 1997), 231–250.
• Freitag and Plassmann (2000) Lori A Freitag and Paul Plassmann. 2000. Local optimization-based simplicial mesh untangling and improvement. Internat. J. Numer. Methods Engrg. 49, 1-2 (2000), 109–125.
• Garanzha (2000) VA Garanzha. 2000. The barrier method for constructing quasi-isometric grids. Computational Mathematics and Mathematical Physics 40 (2000), 1617–1637.
• Garanzha et al. (2014) V.A. Garanzha, L.N. Kudryavtseva, and S.V. Utyuzhnikov. 2014. Variational method for untangling and optimization of spatial meshes. J. Comput. Appl. Math. 269 (2014), 24 – 41.
• Gregson et al. (2011) J. Gregson, A. Sheffer, and E. Zhang. 2011. All-Hex Mesh Generation via Volumetric PolyCube Deformation. Computer Graphics Forum (Special Issue of Symposium on Geometry Processing 2011) 30, 5 (2011).
• Hormann et al. (2008) Kai Hormann, Konrad Polthier, and Alia Sheffer. 2008. Mesh Parameterization: Theory and Practice. In ACM SIGGRAPH ASIA 2008 Courses (Singapore) (SIGGRAPH Asia ’08). Association for Computing Machinery, New York, NY, USA, Article 12, 87 pages.
• Ivanenko (1988) SA Ivanenko. 1988. Construction of nondegenerate grids. Zh. Vychisl. Mat. Mat. Fiz 28, 10 (1988), 1498.
• Jacquotte (1988) Olivier-P Jacquotte. 1988. A mechanical model for a new grid generation method in computational fluid dynamics. Computer methods in applied mechanics and engineering 66, 3 (1988), 323–338.
• Jiang et al. (2017) Zhongshi Jiang, Scott Schaefer, and Daniele Panozzo. 2017. Simplicial Complex Augmentation Framework for Bijective Maps. ACM Trans. Graph. 36, 6, Article 186 (Nov. 2017), 9 pages.
• Jiang et al. (2020) Zhongshi Jiang, Teseo Schneider, Denis Zorin, and Daniele Panozzo. 2020. Bijective Projection in a Shell. ACM Trans. Graph. 39, 6, Article 247 (Nov. 2020), 18 pages.
• Knupp (2001) Patrick M Knupp. 2001. Hexahedral and tetrahedral mesh untangling. Engineering with Computers 17, 3 (2001), 261–268.
• Li et al. (2020) Minchen Li, Zachary Ferguson, Teseo Schneider, Timothy Langlois, Denis Zorin, Daniele Panozzo, Chenfanfu Jiang, and Danny M. Kaufman. 2020. Incremental Potential Contact: Intersection-and Inversion-Free, Large-Deformation Dynamics. ACM Trans. Graph. 39, 4, Article 49 (July 2020), 20 pages.
• Lipman (2012) Yaron Lipman. 2012. Bounded Distortion Mapping Spaces for Triangular Meshes. ACM Trans. Graph. 31, 4, Article 108 (July 2012), 13 pages.
• Nieser et al. (2011) Matthias Nieser, Ulrich Reitebuch, and Konrad Polthier. 2011. CubeCover - Parameterization of 3D Volumes. Computer Graphics Forum (2011).
• Rabinovich et al. (2017) Michael Rabinovich, Roi Poranne, Daniele Panozzo, and Olga Sorkine-Hornung. 2017. Scalable Locally Injective Mappings. ACM Trans. Graph. 36, 2, Article 16 (April 2017), 16 pages.
• Rumpf (1996) Martin Rumpf. 1996. A variational approach to optimal meshes. Numer. Math. 72, 4 (1996), 523–540.
• Schüller et al. (2013) Christian Schüller, Ladislav Kavan, Daniele Panozzo, and Olga Sorkine-Hornung. 2013. Locally Injective Mappings. Computer Graphics Forum (proceedings of Symposium on Geometry Processing) 32, 5 (2013).
• Shen et al. (2019) Hanxiao Shen, Zhongshi Jiang, Denis Zorin, and Daniele Panozzo. 2019. Progressive Embedding. ACM Trans. Graph. 38, 4, Article 32 (July 2019), 13 pages.
• Smith and Schaefer (2015) Jason Smith and Scott Schaefer. 2015. Bijective Parameterization with Free Boundaries. ACM Trans. Graph. 34, 4, Article 70 (July 2015), 9 pages.
• Toulorge et al. (2013) Thomas Toulorge, Christophe Geuzaine, Jean-François Remacle, and Jonathan Lambrechts. 2013. Robust untangling of curvilinear meshes. J. Comput. Phys. 254 (2013), 8–26.
• Tutte (1963) W. T. Tutte. 1963. How to Draw a Graph. Proceedings of the London Mathematical Society s3-13, 1 (01 1963), 743–767.
• Winslow (1966) Alan M Winslow. 1966. Numerical solution of the quasilinear Poisson equation in a nonuniform triangle mesh. Journal of computational physics 1, 2 (1966), 149–172.
• Zhu et al. (1997) C. Zhu, R. H. Byrd, and J. Nocedal. 1997. L-BFGS-B: Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization. ACM Trans. Math. Software 23, 4 (1997), 550–560.

## Appendix A Comprehensive design formulæ Figure 13. On each simplex the map →u(→x) is affine and is entirely defined by the position of the vertices of the domain simplex {→xi} and its image {→ui}.

Given a map , let us denote by the tangent basis, i.e. vectors forming the columns of the Jacobian matrix . For example, in 2D we have and . Let us denote by the dual basis, i.e. vectors chosen in the way that for all indices . In particular, for the 2D settings the dual basis can be written as and . In the 3D case , where is cyclic permutation from . It is a handy choice of variables, in particular, and . For further simplification of notations we will use for , for and .

In order to derive expressions for the gradient and the Hessian matrix of , we write down explicitly the Jacobian matrix for the affine map of a simplex with vertices :

where

 S:=(→x1−→x0 →x2−→x0…→xd−→x0),  detS>0

is the shape matrix, are vertices of “ideal” or “target” shape for the image of the simplex , and is a matrix defined as

 Z:={zij}:=(−1…−1I)S−1

Since the Jacobian matrix is a linear function of , we have

 ∂→ai∂→u⊤j=zjiI,i=1,…,d, j=0,…,d.

The additive contribution to gradient of from the simplex can be written using correspondence of local indices and global indices in the list of vertices:

 (∇F)gj+= detSd!d∑i=1∂→a⊤i∂→uj∂ϕ∂→ai= = detSd!d∑i=1zji∂ϕ∂→ai, j=0,…,d,

where function is defined in § 4.1. Let us provide an explicit expression for :

 ∂ϕ∂→ai=2χ2d→ai−1χ(2dfεχ′−2λdetJ+λgεχ′)→bi

### a.2. Hessian

The blocks of the nonnegative definite part of Hessian matrix of can be updated using the following general formula

 H+gjgi+=detSd!∑m,l∂→a⊤m∂→ujM+ml∂→al∂→u⊤i,

where denotes a block of positive definite matrix defined in Eq. (9). Let us provide an explicit expression for the matrix:

 M+=(Ib)⎛⎜⎝∂2Φ∂a∂a⊤∂2Φ∂a∂D∂2Φ∂D∂a⊤∂2Φ∂D2⎞⎟⎠(Ib)⊤, where
 ∂2Φ∂a∂a⊤ =2χ2dI ∂2Φ∂D2 =2d(1+2d)|a|2χ′2χ2+2d+λ(2χ−4Dχ′χ2+2(1+D2)χ′2χ3) ∂2Φ