Foldover-free maps in 50 lines of code
Vladimir Garanzha, Igor Kaporin, Liudmila Kudryavtseva, François Protais, Nicolas Ray, Dmitry Sokolov
FFoldover-free maps in 50 lines of code
VLADIMIR GARANZHA ∗ , IGOR KAPORIN ∗ , and LIUDMILA KUDRYAVTSEVA ∗ , Dorodnicyn Computing Cen-ter FRC CSC RAS, Moscow, Russia, Moscow Institute of Physics and Technology, Moscow, Russia
FRANÇOIS PROTAIS † , NICOLAS RAY, and DMITRY SOKOLOV, Université de Lorraine, CNRS, Inria, LORIA,F-54000 Nancy, France
Fig. 1. Our method of constructing injective maps opens a door for a large variety of applications. This figure shows an example of a thick prismatic meshlayer (shown in green) built around a triangulated surface, a very challenging problem for highly curved objects. Thanks to our method, we are able to computesuch a layer free of folds and self-intersections.
Mapping a triangulated surface to 2D space (or a tetrahedral mesh to 3Dspace) is the most fundamental problem in geometry processing. In com-putational 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 mappingwhere the geometry of the object is to be defined in the map space and thegeometric domain is not explicit, supposing that each element is regular.In this paper, we propose a mapping method inspired by the untanglingproblem and compare its performance to the state of the art. The mainadvantage of our method is that the untangling aims at producing locallyinjective maps, which is the major challenge of mapping. In practice, ourmethod produces locally injective maps in very difficult settings, and withless distortion than the previous work, both in 2D and 3D. We demonstrateit on a large reference database as well as on more difficult stress tests. For abetter reproducibility, we publish the code in Python for a basic evaluation,and in C++ for more advanced applications.CCS Concepts: • Computing methodologies → Mesh models . Additional Key Words and Phrases:
Parameterization, injective mapping,mesh untangling
Most geometric objects are represented by a triangulated surfaceor a tetrahedral mesh. The mapping problem consists in generatinga 2D or 3D map of these objects. This is a fundamental problem ofcomputer graphics because it is much easier for many applicationsto work in this map space than to directly manipulate the objectitself. To give few examples, texture mapping stores colors of asurface as images in the map space, remeshing uses global mapsin 2D [Bommes et al. 2013] and 3D [Gregson et al. 2011; Nieser ∗ This work is supported by the Ministry of Science and Higher Education of the RussianFederation, project No 075-15-2020-799 † Corresponding author: [email protected] are listed in alphabetical order. et al. 2011]. In addition, mapping algorithms can be used to deformvolumes [Li et al. 2020], or generate shells from surfaces [Jiang et al.2020].
What is a good map?
Most often maps are represented by theposition of the vertices in the map space, and interpolated linearlyon each element (triangle or tetrahedron). In a perfect world, themap space would keep the geodesic distances of the object. Un-fortunately, this is usually impossible due to Gaussian curvature,and application-specific constraints such as constrained position ofvertices or overlaps in the map space. Therefore, the objective of themapping algorithms is to minimize the distortion between geomet-ric and map spaces, opening the door for numerical optimizationapproaches as detailed in the survey [Hormann et al. 2008].
What about the invertibility?
Unfortunately, when high distortionis required to satisfy the constraints, these algorithms may lose thefundamental property of a maps: injectivity. A solution to preserveit [Floater 1997] relies on Tutte’s theorem [Tutte 1963], however thesurface boundary must be mapped to a convex polygon. Despitethis strong limitation, it still remains the reference algorithm to gen-erate injective maps. Lower distortion can be obtained by changingweights of the barycentric coordinates [Eck et al. 1995] (as longas they are not negative), and alternative solutions [Campen et al.2016; Shen et al. 2019] have been explored to improve robustness tonumerical imprecision by modifying the mesh connectivity.
Local invertibility.
In many applications, maps are used to accessa neighborhood of a point within a coherent local coordinate system.To this end, global injectivity is not required, and we instead lookfor local injectivity [Schüller et al. 2013; Smith and Schaefer 2015].Their approach starts from an injective map [Floater 1997], andmaintains the local injectivity when minimizing the distortion andenforcing the constraints. This allows them to optimize at the same
Technical Report, 2021. a r X i v : . [ c s . C G ] F e b • Vladimir Garanzha, Igor Kaporin, Liudmila Kudryavtseva, François Protais, Nicolas Ray, and Dmitry Sokolov 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 recoveredfor a map with few foldovers present. For example, in 2D [Lipman2012] and 3D [Aigerman and Lipman 2013], the map is projectedon a class of bounded distortion maps. The numerical method ishowever unlikely to succeed for stiff problems. Recovering localinjectivity is also known as mesh untangling. Originally related toArbitrary Lagrangian-Eulerian moving mesh approach, the meshuntangling problem considers a simplicial complex with misorientedelements and attempts to flip them back by optimizing the positionof the vertices. There is an abundant literature on mesh untangling[Du et al. 2020; Escobar et al. 2003; Freitag and Plassmann 2000;Knupp 2001; Toulorge et al. 2013], however the common opinionis that untangling is a very hard problem and algorithms are notrobust enough. As a manifestation of frustration over this problem[Danczyk and Suresh 2013] investigates a finite element methodworking directly on tangled ( sic! ) meshes.
Elastic deformations.
To recover local injectivity, we propose amethod stemming from the computational physics. It is very im-portant to note that there is rich literature on mesh deformation inthe community working on grid generation for scientific computa-tion. Numerical simulation of hydrodynamic instability of layeredstructures requires sound mathematical foundations behind movingdeforming mesh algorithms. In the ’60s Winslow and Crowley, inde-pendently one from another, introduced mesh generation methodsbased on inverse harmonic maps [Crowley 1962; Winslow 1966].Since then, a lot of effort was spent on mesh generation based onelastic deformations [Jacquotte 1988], but mostly for regular grids. In1988, at the time of domination of finite difference mapped grid gen-eration methods, S. Ivanenko introduced the pioneering concept ofbarrier variational grid generations methods guaranteeing construc-tion of non-degenerate grids [Charakhch’yan and Ivanenko 1997;Ivanenko 1988]. To generate deformations with bounded global dis-tortion (bounded quasi-isometry constant), Garanzha proposed tominimize an elastic energy for a hyperelastic material with stiffeningsuppressing singular deformations [Garanzha 2000]. Invertibilitytheorem for deformation of this material was established in the 3Dcase 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 onthe concept of polyconvex distortion energies [Ball 1976]. He notonly proved Weierstrass-style existence theorem for this class ofvariational problems, but also formulated a theorem on invertibilityof elastic deformations for quite general 3D domains [Ball 1981].It is important to note that Ball invertibility theorem is proved forSobolev mappings and can be applied directly for finite elementspaces, i.e. to deformation of meshes, as was pointed out in [Rumpf1996].
Our contributions.
Inspired by these results on untangling andelastic deformations, we propose a method outperforming recentstate 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 di-rectly, with a parameter for a tradeoff between the angle andarea preservation, • we provide a rigourous analysis of the foundations of ournumerical resolution scheme, • our method supports mapping with free boundaries, • we observe better robustness for high distortions as well aswith 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 inthe supplemental material for more advanced applications.The rest of the paper is organized as follows: we start with pre-senting our method in § 2, then we evaluate its performance (§ 3.1and § 3.2) as well as its limitations (§ 3.3). Then we present theoreti-cal guarantees for our resolution scheme: in § 4.1 we prove that ourapproximation of Hessian matrix is definite positive, and finally in§ 4.2 we prove that our choice of the regularization parameter se-quence guarantees that a minimization algorithm can find a meshfree of inverted elements in a finite number of steps. In this rather short section we present our method of computing afoldover-free map (cid:174) 𝑢 : Ω ⊂ R 𝑑 → R 𝑑 , i.e. we map the domain Ω toa parametric domain. This presentation is unified both for 2D and3D settings, and by 𝑑 we denote the number of dimensions; in ournotations we use arrows for all vectors of dimension 𝑑 .The section is organized as follows: in § 2.1 we give a primeron the variational formulation of mapping problem in continuoussettings, then we state our problem in § 2.2 as a regularization ofthis variational formulation, and finally we present our numericalresolution scheme in § 2.3. Let us denote by (cid:174) 𝑢 ( (cid:174) 𝑥 ) a map to a parametric domain: for the flat 2Dcase we can write as (cid:174) 𝑢 ( 𝑥, 𝑦 ) = ( 𝑢 ( 𝑥, 𝑦 ) , 𝑣 ( 𝑥, 𝑦 )) , and for a 3D map (cid:174) 𝑢 ( 𝑥, 𝑦, 𝑧 ) = ( 𝑢 ( 𝑥, 𝑦, 𝑧 ) , 𝑣 ( 𝑥, 𝑦, 𝑧 ) , 𝑤 ( 𝑥, 𝑦, 𝑧 )) .Consider the following variational problem:arg min (cid:174) 𝑢 ∫ Ω ( 𝑓 ( 𝐽 ) + 𝜆𝑔 ( 𝐽 )) 𝑑𝑥, (1)where 𝐽 is the Jacobian matrix of the mapping (cid:174) 𝑢 ( (cid:174) 𝑥 ) , and 𝑓 ( 𝐽 ) : = (cid:40) tr 𝐽 ⊤ 𝐽 ( det 𝐽 ) 𝑑 , det 𝐽 > +∞ , det 𝐽 ≤ 𝑔 ( 𝐽 ) : = (cid:26) det 𝐽 + 𝐽 , det 𝐽 > +∞ , det 𝐽 ≤ 𝑓 ( 𝐽 ) and 𝑔 ( 𝐽 ) have concurrent goals,one preserves angles and the other preserves the area, and thus 𝜆 serves as a trade-off parameter. The minimization algorithm is subject to conditions of Th. (1); from our numericalexperiments we observe that our minimization algorithm almost always respects theconditions.Technical Report, 2021. oldover-free maps in 50 lines of code • 3
Fig. 2. Regularization function for the denominator in Eq. (3) . When 𝜀 tendsto zero, 𝜒 ( 𝜀, 𝐷 ) tends to 𝐷 for positive values of 𝐷 , and to + for negativevalues of 𝐷 . As a side note, with 𝜆 = 𝑑 =
2, Prob. (1) presents a varia-tional formulation of an inverse harmonic map problem. Namely, ifwe write down the Euler-Lagrange equations for Prob. (1) and in-terchange the dependent and independent variables , we obtain theLaplace equation Δ (cid:174) 𝑥 ( (cid:174) 𝑢 ) = (cid:174) Δ (cid:174) 𝑢 ( (cid:174) 𝑥 ) = (cid:174) Ω 𝐽 ( (cid:174) 𝑢 ) > We can avoid nonpositive denominators in 𝑓 and 𝑔 using a regular-ization function 𝜒 for a positive value of 𝜀 (Fig. 2): 𝜒 ( 𝐷, 𝜀 ) : = 𝐷 + √ 𝜀 + 𝐷 , (2)Then we define a regularized version 𝑓 𝜀 , 𝑔 𝜀 of functions 𝑓 and 𝑔 : 𝑓 𝜀 ( 𝐽 ) : = tr 𝐽 ⊤ 𝐽 ( 𝜒 ( det 𝐽, 𝜀 )) 𝑑 , 𝑔 𝜀 ( 𝐽 ) : = det 𝐽 + 𝜒 ( det 𝐽, 𝜀 ) , (3)so that Prob (1) is reformulared aslim 𝜀 → + arg min (cid:174) 𝑢 ∫ Ω ( 𝑓 𝜀 ( 𝐽 ) + 𝜆𝑔 𝜀 ( 𝐽 )) 𝑑𝑥 (4) Attention, this step assumes that the solution of Prob. (1) is a diffeomorphism.
Under certain assumptions solutions of Prob. (4) are solutionsof Prob. (1), however, Prob. (4) does offer a way of getting rid offoldovers if a foldover-free initialization is not available.In practice, the map (cid:174) 𝑢 is piecewise affine with the Jacobian matrix 𝐽 being piecewise constant, and can be represented by the coor-dinates of the vertices in the parametric domain {(cid:174) 𝑢 𝑖 } 𝑉𝑖 = . Let usdenote the vector of all variables as 𝑈 : = (cid:16) (cid:174) 𝑢 ⊤ . . . (cid:174) 𝑢 ⊤ 𝑉 (cid:17) ⊤ , then ouroptimization problem can be discretized as follows:lim 𝜀 → + arg min 𝑈 𝐹 ( 𝑈 , 𝜀 ) , (5)where 𝐹 ( 𝑈 , 𝜀 ) : = 𝑇 ∑︁ 𝑡 = ( 𝑓 𝜀 ( 𝐽 𝑡 ) + 𝜆𝑔 𝜀 ( 𝐽 𝑡 )) vol ( 𝑇 𝑡 ) , 𝑉 is the number of vertices, 𝑇 is the number of simplices, 𝐽 𝑡 is theJacobian matrix for the simplex 𝑡 and vol ( 𝑇 𝑡 ) is the volume of thesimplex 𝑇 𝑡 in the original domain. To solve Prob. (5), we use an iterative descent method. We startfrom an initial guess 𝑈 , and we build a sequence of approximations 𝑈 𝑘 + : = 𝑈 𝑘 + Δ 𝑈 𝑘 . For each iteration we need to carefully choosethe regularization parameter 𝜀 𝑘 . Starting from 𝜀 : =
1, we definethe sequence as follows: 𝜀 𝑘 + : = (cid:18) − 𝜎 𝑘 √ ( 𝐷 𝑘 + − ) +( 𝜀 𝑘 ) | 𝐷 𝑘 + − |+ √ ( 𝐷 𝑘 + − ) +( 𝜀 𝑘 ) (cid:19) 𝜀 𝑘 , if 𝐷 𝑘 + − < , ( − 𝜎 𝑘 ) 𝜀 𝑘 , if 𝐷 𝑘 + − ≥ , (6)where 𝐷 𝑘 + − : = min 𝑡 ∈ ... 𝑇 det 𝐽 𝑘 + 𝑡 is the minimum value of the Jacobiandeterminant over all cells of the mesh at the iteration 𝑘 +
1, and 𝜎 𝑘 : = max (cid:16) , − 𝐹 ( 𝑈 𝑘 + ,𝜀 𝑘 ) 𝐹 ( 𝑈 𝑘 ,𝜀 𝑘 ) (cid:17) . Note that Eq. (6) is justified by Th. 1(§ 4.2) on finite untangling sequence; while this formula providessome guarantees, one may use a simpler heuristic to accelerate thecomputations, more on this in § 3.2.The simplest way to find Δ 𝑈 𝑘 is to call a quasi-Newtonian solversuch as L-BFGS-B [Zhu et al. 1997]. The only thing we need toimplement is the computation of the function 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 ) and itsgradient ∇ 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 ) .Another option is to compute analytically the Hessian matrixinstead of estimating it. The problem, however, is that the Hessianmatrix 𝜕 𝐹𝜕 𝑈 𝜕 𝑈 ⊤ is not positive definite. In this paper we propose itsapproximation that ensures the positive definiteness. The modifiedHessian matrix 𝐻 + ( 𝑈 𝑘 , 𝜀 𝑘 ) of the function 𝐹 with respect to 𝑈 atthe point 𝑈 𝑘 is built out of 𝑑 × 𝑑 blocks 𝐻 + 𝑖 𝑗 ≈ 𝜕 𝐹𝜕 (cid:174) 𝑢 𝑖 𝜕 (cid:174) 𝑢 ⊤ 𝑗 ( 𝑈 𝑘 , 𝜀 𝑘 ) . Here, the matrix 𝐻 + 𝑖 𝑗 is placed on the intersection of 𝑖 -th blockrow and 𝑗 -th block column; the ≈ symbol means that we removeall the terms depending on the second derivative of 𝜒 and second The fact that the solution of Prob. (4) is a diffeomorphism is sufficient (but not neces-sary) for the equivalence. Technical Report, 2021. • Vladimir Garanzha, Igor Kaporin, Liudmila Kudryavtseva, François Protais, Nicolas Ray, and Dmitry Sokolov
ALGORITHM 1:
Computation of a locally injective map
Input: 𝑈 ; // initial guess (vector of size 𝑉 × 𝑑 ) Input: useQuasiNewton ; // boolean to choose the optimization scheme
Output: 𝑈 ; // final locally injective map (vector of size 𝑉 × 𝑑 ) 𝑘 ← repeat compute 𝜀 𝑘 ; // regularization parameter, Eq. (6) [variant: Eq. (7) ] if useQuasiNewton then 𝑈 𝑘 + ← L-BFGS-B ( 𝑈 𝑘 , 𝜀 𝑘 ) ; // inner L-BFGS-B loop else compute a modified Hessian matrix 𝐻 + ( 𝑈 𝑘 , 𝜀 𝑘 ) ; Δ 𝑈 𝑘 ← ( 𝐻 + ) − ∇ 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 ) ; // conjugate gradients 𝑈 𝑘 + ← arg min 𝜏 𝐹 ( 𝑈 𝑘 + 𝜏 Δ 𝑈 𝑘 , 𝜀 𝑘 ) ; // line search end 𝑘 ← 𝑘 + until min 𝑡 ∈ ... 𝑇 det 𝐽 𝑘𝑡 > and 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 ) > ( − − ) 𝐹 ( 𝑈 𝑘 − , 𝜀 𝑘 − ) ; 𝑈 ← 𝑈 𝑘 ; Fig. 3. The input mesh with foldovers and the untangling produced by theListings 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.Fig. 4. Injective mapping sanity check: make two points inside a squareswitch places. Left: the input problem, all locked points are shown in red.Right: foldover-free result obtained with our method. derivatives of det 𝐽 to keep 𝐻 + positive definite. Refer to Appendix Afor the formulæ, and to § 4.1 for the proof of the positive definitenessof 𝐻 + .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 Pythonimplementation and the corresponding input/output generated bythe code. Note that our method is not limited to simplicial meshesonly: in this particular example we evaluate the Jacobian matrixfor every triangle forming quad corners, what corresponds to thetrapezoidal quadrature rule. In this section we provide an experimental evaluation of the method.In the field of computer graphics, any claim about map injectivityalways faces a simple sanity check (Fig. 4): take a square and maketwo points inside switch places. Our method successfully avoids thedesk-reject, so we start this section (§ 3.1) by testing our method onthe benchmark [Du et al. 2020], then we continue with further testswe have found relevant (§ 3.2), and finally we discuss the limitationsof the approach in § 3.3.
Along with their paper, Du et al. have published a valuable bench-mark database. It contains a huge number of 2D and 3D constrainedboundary injective mapping challenges. For 2D challenges, thebenchmark contains 3D triangulated surfaces to flatten, and not flat2D meshes as we have described in § 2.1. Nevertheless, our methodcan handle it directly because the mapping is still R → R oneach triangle. To the best of our knowledge, Total Lifted Content(TLC) [Du et al. 2020] and our method are the only ones passing thebenchmark without any fail.A representative example from the database is given in top rowof Fig. 5. The challenge is to map the “Lucy” mesh statuette fromthe Stanford Computer Graphics Laboratory to a P-shaped domain.This mesh has a topology of a disk, and its boundary vertices areuniformly spaced on the P-shape boundary. As an initializationto the problem, Du et al. have computed the corresponding (flat)minimal surface that obviously contains a foldover (Fig. 5–top left).Then the problem boils down to a mesh untangling with lockedboundary. Mapping quality measure.
How to measure quality of a map? Well,it depends on the goal. An identity is an unreachable ideal; tradi-tional competing goals are (as much as possible) angle preservingand area preserving maps. Thus, we can measure the extreme valuesof the failure of a map to be conformal or authalic. Our maps beingpiecewise affine, the Jacobian matrix 𝐽 is constant per element. Letus define the largest singular value of 𝐽 as 𝜎 ( 𝐽 ) , and the smallestsingular value as 𝜎 𝑑 ( 𝐽 ) ; then the quality of a mapping can be re-duced to extreme values of the stretch (max 𝜎 ( 𝐽 ) 𝜎 𝑑 ( 𝐽 ) ) and the scaling(min det 𝐽 ).For the “Lucy-to-P” challenge (Fig. 5–top row) our map differsfrom the TLC result by 12 orders of magnitude in terms of minimumscaling, and by two orders of magnitude in terms of maximumstretch. To visualize this difference in scaling, we have providedthe close-ups: Fig. 5–top middle shows a map of the Lucy’s torch,whereas the same level of zoom on the result by Du et al. (Fig. 5–topright) contains not only the torch, but also both wings, the head andthe right arm!Note also that the input “Lucy” mesh is slightly anisotropic; ourmethod allows us to prescribe the element target shape, so the dresspleats are clearly visible in our mapping. Benchmark database.
Our method successfully passes all chal-lenges from the benchmark [Du et al. 2020]. The benchmark con-sists of 10743 meshes to untangle in 2D and 904 meshes in 3D underlocked boundary constraints. In Fig. 7 we provide quality plots
Technical Report, 2021. oldover-free maps in 50 lines of code • 5
Initialization Our result [Du et al. 2020] min det 𝐽 ≈ · − max 𝜎 ( 𝐽 ) 𝜎 ( 𝐽 ) ≈ 𝐽 ≈ · − max 𝜎 ( 𝐽 ) 𝜎 ( 𝐽 ) ≈ 𝐽 ≈ · − max 𝜎 ( 𝐽 ) 𝜎 ( 𝐽 ) ≈
70 min det 𝐽 ≈ · − max 𝜎 ( 𝐽 ) 𝜎 ( 𝐽 ) ≈ 𝐽 ≈ · − max 𝜎 ( 𝐽 ) 𝜎 ( 𝐽 ) ≈ Fig. 5. Constrained boundary injective mapping challenge proposed by [Du et al. 2020]: the “Lucy” mesh is mapped to a P-shaped domain by constraining thevertices shown in red. Left column: three different initializations for the same problem. Middle column: our method produces the same (up to a numericalprecision) result on all three initializations. Right column: total lifted content method [Du et al. 2020] fails to solve for the randomly initialized interiorvertices, and produces very different results on other two initializations. Three miniature images of “Lucy” show in red the portion of the surface visible in thecorresponding close-ups. max 𝜎 ( 𝐽 ) 𝜎 ( 𝐽 ) ≈ 𝐽 ≈ . (a) max 𝜎 ( 𝐽 ) 𝜎 ( 𝐽 ) ≈ 𝐽 ≈ . (b) max 𝜎 ( 𝐽 ) 𝜎 ( 𝐽 ) ≈ 𝐽 ≈ . (c) max 𝜎 ( 𝐽 ) 𝜎 ( 𝐽 ) ≈ 𝐽 ≈ . (d) max 𝜎 ( 𝐽 ) 𝜎 ( 𝐽 ) ≈ 𝐽 ≈ . (e) max 𝜎 ( 𝐽 ) 𝜎 ( 𝐽 ) ≈ 𝐽 ≈ . (f) Fig. 6. Constrained boundary injective mapping stress test. We have generated an isotropic tetrahedral mesh of a cube subtracted from a larger cube, andwe rotate the inner cube’s boundary to test the robustness. Top row and bottom row correspond to two different slices of the same mesh. Columns (a)–(d) :injective maps produced by our method, columns (e) and (f) : injective maps produced by [Du et al. 2020]. The method by [Du et al. 2020] fails to generateinjective maps for 135° and 180° inner cube rotations. The colormap illustrates the relative volume scaling: green for det 𝐽 ≈ , red for inflation, blue forcompression. Technical Report, 2021. • Vladimir Garanzha, Igor Kaporin, Liudmila Kudryavtseva, François Protais, Nicolas Ray, and Dmitry Sokolov
2D dataset (10743 challenges)3D dataset (904 challenges)
Fig. 7. Quality plot of the resulting locally injective maps for every challengefrom 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 toa quality of the corresponding map reduced to two numbers: the maximumstretch and the mininum scale.
Top row: mapping quality on the 2D dataset.
Bottom row: mapping quality on the 3D dataset. Left column shows theabsolute maximum stretch and absolute minimum scale, whereas the rightcolumn shows the maximum stretch and minimum scale for the top 95% ofmeasurements. of the resulting locally injective maps. These are log − log scatterplots: each dot corresponds to a quality of the corresponding mapreduced to two numbers: the maximum stretch (max 𝜎 ( 𝐽 ) 𝜎 𝑑 ( 𝐽 ) ) andthe minimum scaling (min det 𝐽 ). Left column of Fig. 7 shows theworst quality measurements for every 2D problem (top) as wellas for every 3D challenge (bottom) of the dataset. Our results areshown in blue, whereas TLC results are shown in red. To illustratethe distribution of the elements’ quality, for each injective map wehave removed 5% of worst measurements: the right column of Fig. 7shows the maximum stretch and the minimum scaling for the top95% of measurements.Note the dot arrangements forming lines in the plot: these dotscorrespond to the few sequences of deformation present in thedatabase. Timings.
Fig. 8 provides a log − log scatter plot of our runningtime 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 toseveral minutes for the largest meshes. These times were obtainedwith a 12 cores i7-6800K CPU @ 3.40 GHz. As in Fig. 8, the ver-tical lines in the 3D dataset plot correspond to the sequences ofdeformation in the benchmark. Fig. 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 and904 in 3D). Blue dots show the running times obtained using Eq. (7) , greendots correspond to Eq. (6) . There are two scatter plots superposed, both represent the sameresolution scheme with an exception corresponding to the way wecompute 𝜀 𝑘 (Alg. 1–line 3). The green scatter plot corresponds to aconservative update rule (Eq. (6)) offering guarantees on untanglingin a finite number of steps (refer to Th. 1), whereas the blue scatterplot is obtained using the following update rule: 𝜀 𝑘 : = √︄ − + · − · (cid:20) min ( , min 𝑡 ∈ ... 𝑇 det 𝐽 𝑘𝑡 ) (cid:21) , (7)This formula was chosen empirically, and it performs very well inthe vast majority of situations. In particular, it allows for all thedatabase [Du et al. 2020] to pass the injectivity test. Sensitivity to initialization.
Our next test is the sensitivity to theinitialization. We have generated two other initializations for the“Lucy-to-P” challenge: the one with all interior vertices collapsedonto a single point (Fig. 5–middle row), and with the interior verticesbeing randomly placed withing a bounding square (Fig. 5–bottomrow).Our method produces virtually the same result on all three initial-izations, whereas TLC generates very different results for the firsttwo, and fails for the third one. It is interesting to note that TLC isheavily depending on the initialization, it alters very little the inputgeometry.
Large deformation stress test.
For our next test we have generatedan isotropic tetrahedral mesh of a cube with a cavity, and we ro-tated the inner boundary to test the robustness of our method tolarge deformations. Figure 6 shows the results. Our L-BFGS-basedoptimization scheme succeeds up to the rotation of 135°, and wehad to switch to the Newton method to reach the 180° rotation. TLCmethod had succeded on 45° and 90°, and failed for the 135° and180°. Note that as in the previous test, even when the untanglingsucceeds, TLC alters very little the input map, thus producing heav-ily stretched tetrahedra, whereas our method evenly dissipates thestress over all the domain.
Free boundary injective mapping.
To the best of our knowledge,our method is the only one able to produce inversion-free maps
Technical Report, 2021. oldover-free maps in 50 lines of code • 7 (a) (b)
Fig. 9. Free boundary injective mapping. The vertices shown in red areconstrained, all other vertices are free to move. (a): a compression test, (b): a bend test. Refer to Fig 10–a for the rest shape. (a) (b)(c) (d)
Fig. 10. Free boundary injective mapping: influence of the parameter 𝜆 inProb. (5) . The vertices shown in red are constrained, all other vertices arefree to move. (a): the rest shape; (b): a stretch preserving the shape of theelements ( 𝜆 = ); (c): an area-preserving stretch ( 𝜆 = ); (c): a trade-offbetween the shape and the area preservation ( 𝜆 = ). 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 shapebeing compressed and the same shape being bent. The boundary isfree to move, we lock the vertices shown in red. Refer to Fig. 10-a forthe rest shape. The shape behaves exactly as a human would expectit: upon compression the shape chooses one of the two possibleresults (Fig. 9–a), and successfully passes the bend test (Fig. 9–b),note the geometrical details that are naturally rotated. Shape-area tradeoff 𝜆 . Our final tests illustrate the influence ofthe parameter 𝜆 in Prob. (5) on the resulting map. We have com-puted three free boundary maps of the rest shape (Fig. 10–a) beingstretched. First we chose 𝜆 =
0, that is, only the shape quality termis 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 forthe resulting map. Naturally, an area preserving map ( 𝜆 = ) mustdeform the elements to satisfy the area constraint (Fig. 10–c). Finally,in Fig. 10–d we show an example with a tradeoff between the areaand angles preservation. While globally performing very well in practice, our method stillpresents some limitations. We have two main sources of limita-tions: theoretical limitations as well as very practical ones relatedto numerical stability of our resolution scheme.
Overlaps.
First of all, an inversion-free map does not imply globalinjectivity. Fig. 11–a provides an example of an inversion-free mapwith two cases of non-injectivity when optimizing for a map withfree boundaries: the map can present global overlaps as well asthe boundary can “wind up” around boundary vertices, i.e. thetotal angle of triangles incident to a vertex can be superior to 2 𝜋 .Moreover, while being less frequent, similar situations may occur oninterior vertices. Typically this situation happens near constraintscausing a local compression in the shape.Let us illustrate this behavior on a very simplistic mesh consistingof a single fan of 12 triangles. All vertices are free to move, the targetshape is set to be the unit equilateral triangle for all elements. For thisproblem Fig. 12–a shows a local minimum, and the Fig. 12–c showsthe global minimum respecting perfectly the prescribed total angleof 4 𝜋 around the center vertex. Both are inversion-free maps, butonly the map in Fig. 12–a is a globally injective one. Depending onthe initialization and the resolution scheme chosen, we can convergeto either minimum. Note, however, that the center vertex has thewinding number 1 in one map and 2 in the other, and thus we cannot deform continuously one to the other without inverting someelements. Note also that the configurations like in the Fig. 12–bpresent inverted elements and thus can not be generated by ourmethod.It is possible to avoid all overlaps altogether by embedding ourshape 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 materialmust not have inverted elements (refer to Fig. 11–b). The thickprismatic layer in Fig. 1 was generated by a similar procedure: wehave generated a very thin layer of triangular prisms around thedragon, and tetrahedralized the exterior bounded by a cube. Aftercalling the untangling procedure, we have obtained an offset surfacewith exactly the same mesh connectivity as the original dragonmesh.While this embedding kind of approach works well for certainapplications, for other it may be hard to apply. We are currentlyexploring simpler practical ways to fix the winding problem. Numerical challenges.
Even when the problem is well-posed, a ro-bust implementation may present significant difficulties. As we havesaid above, the quasi-Newtonian optimization scheme performs wellfor “simple” problems (it passes all the benchmark database!), butmay fail for large deformations, where Newton iterations are nec-essary. While our modified Hessian matrix is symmetric positivedefinite, note that for stiff problems the Jacobi preconditioned conju-gate gradients can fail and one might need the incomplete Choleskidecomposition and beyond.In practice we have found the method being very robust in 2Dsettings: we have not encountered a practical test case we werenot able to treat with our method. In 3D, however, it can fail due
Technical Report, 2021. • Vladimir Garanzha, Igor Kaporin, Liudmila Kudryavtseva, François Protais, Nicolas Ray, and Dmitry Sokolov (a) (b)
Fig. 11. Free boundary mapping limitations. (a): this mesh presents twokinds of problems, namely, a global overlap and the mesh wrapped arounda boundary vertex. (b):
Both problems can be avoided by embedding themesh into an outer triangulation. (a) (b) (c)
Fig. 12. Free boundary mapping limitations: three maps of a very simplisticmesh made of 12 triangles. (a) and (c) both are inversion-free maps andthus allowed by our method, whereas the map (b) has inverted elements,and thus is prohibited by our method. to the numerical challenges in very anisotropic and highly twistedmeshes.
This section presents in two parts a rigorous analysis of the penaltymethod. First in § 4.1 we prove that the modified Hessian matrix 𝐻 + ( 𝑈 , 𝜀 ) is indeed positive definite, and then in § 4.2 we show theorigins 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 min det 𝐽 > 𝐾 < ∞ thesolution arg min 𝑈 𝐹 ( 𝑈 , 𝜀 𝐾 ) belongs to the admissible set. 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 𝑗 -thblock column. It is a common practice to add some regularizationterms to the Hessian matrix to make it positive definite, but wepropose to modify the finite element (FE) matrix assembly procedureby eliminating some terms potentially leading to an indefinite FEmatrix. To this end, we restrict our attention to a single simplex and westudy a function 𝜙 ( 𝐽 ) of the Jacobian matrix defined as follows: 𝜙 ( 𝐽 ) : = 𝑓 𝜀 ( 𝐽 ) + 𝜆𝑔 𝜀 ( 𝐽 ) = tr 𝐽 ⊤ 𝐽 ( 𝜒 ( det 𝐽, 𝜀 )) 𝑑 + 𝜆 det 𝐽 + 𝜒 ( det 𝐽, 𝜀 ) (8)Let us denote by 𝑎 ∈ R 𝑑 the (columnwise) flattening of theJacobian matrix 𝐽 , i.e. the vector composed of the elements of 𝐽 .We decompose the 𝑑 × 𝑑 Hessian matrix of 𝜙 with respect to theJacobian matrix entries into two parts: 𝜕 𝜙𝜕 𝑎 𝜕 𝑎 ⊤ = 𝑀 + + 𝑀 ± , where 𝑀 + is a positive definite matrix, and 𝑀 ± can be an indefinite matrixthat we neglect. The matrix 𝑀 ± contains all terms depending on 𝜒 ′′ and second derivatives of det 𝐽 with respect to elements of theJacobian matrix 𝐽 . Our map is affine on the simplex of interest,therefore its Jacobian matrix 𝐽 is a linear function of the vertices ofthe simplex. Thus, the idea is to compute a positive definite matrix 𝑀 + ( 𝐽 ) , and use the chain rule to get the Hessian matrix with respectto our variables 𝑈 and assemble the matrix 𝐻 + .So, we choose some arbitrary point 𝐽 and we want to show theway to decompose 𝜕 𝜙𝜕 𝑎 𝜕 𝑎 ⊤ ( 𝐽 ) into a sum of 𝑀 + ( 𝐽 ) and 𝑀 ± ( 𝐽 ) with 𝑀 + ( 𝐽 ) >
0. To do so, first we write down the first orderTaylor expansion 𝑞 ( 𝐷 ) of the function 𝜒 ( 𝐷, 𝜀 ) around some point 𝐷 = det 𝐽 : 𝑞 ( 𝐷 ) : = 𝜒 ( 𝐷 , 𝜀 ) + 𝜕 𝜒𝜕 𝐷 ( 𝐷 , 𝜀 )( 𝐷 − 𝐷 ) . Next we define a function Φ ( 𝑎, 𝐷 ) as follows: Φ ( 𝑎, 𝐷 ) : = | 𝑎 | ( 𝑞 ( 𝐷 )) 𝑑 + 𝜆 𝐷 + 𝑞 ( 𝐷 ) . Note that Φ differs a bit from 𝜙 : it has one more argument and 𝜒 is replaced by its linearization in the denominator. While thismaneuver 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 verifythat the following decomposition holds: 𝜕 𝜙𝜕 𝑎 𝜕 𝑎 ⊤ ( 𝐽 ) = 𝑀 + ( 𝐽 ) + 𝑀 ± ( 𝐽 ) , (9)where 𝑀 + : = (cid:16) 𝐼 𝜕 𝐷𝜕 𝑎 (cid:17) (cid:32) 𝜕 Φ 𝜕 𝑎 𝜕 𝑎 ⊤ 𝜕 Φ 𝜕 𝑎 𝜕 𝐷𝜕 Φ 𝜕 𝐷 𝜕 𝑎 ⊤ 𝜕 Φ 𝜕 𝐷 (cid:33) (cid:16) 𝐼 𝜕 𝐷𝜕 𝑎 (cid:17) ⊤ , and 𝑀 ± : = 𝜕 Φ 𝜕 𝐷 𝜕 𝐷𝜕 𝑎 𝜕 𝑎 ⊤ − 𝜒 ′′ 𝜒 (cid:18) 𝑑 𝑓 𝜀 + 𝜆𝑔 𝜀 (cid:19) 𝜕 𝐷𝜕 𝑎 𝜕 𝐷𝜕 𝑎 ⊤ . The easiest way to check that the equality (9) holds is to note thatat the point 𝐽 we have 𝑞 = 𝜒 , 𝑞 ′ = 𝜒 ′ , and therefore we have 𝜕 𝜙 ( 𝑎 ) 𝜕 𝑎 = 𝜕 Φ ( 𝑎, 𝐷 ( 𝑎 )) 𝜕 𝑎 + 𝜕 Φ ( 𝑎, 𝐷 ( 𝑎 )) 𝜕 𝐷 𝜕 𝐷 ( 𝑎 ) 𝜕 𝑎 . To calculate the Hessian 𝜕 𝜙𝜕 𝑎 𝜕 𝑎 ⊤ ( 𝐽 ) , it suffices to differentiate thisexpression one more time and add the terms in 𝜒 ′′ that were zeroedout by the linearization.To sum up, in our computations, for each simplex we approximatethe Hessian matrix 𝜕 𝜙𝜕 𝑎 𝜕 𝑎 ⊤ by the 𝑑 × 𝑑 matrix 𝑀 + and we neglectthe term 𝑀 ± . Thanks to the convexity of Φ , it is trivial to verify that Technical Report, 2021. oldover-free maps in 50 lines of code • 9 for any choice of 𝐽 the matrix 𝑀 + is definite positive. Then we usethe chain rule over 𝑀 + to get the Hessian matrix with respect toour variables 𝑈 , and we assemble a ( 𝑑 𝑉 ) × ( 𝑑 𝑉 ) approximation 𝐻 + of the Hessian matrix for the energy function 𝐹 ( 𝑈 , 𝜀 ) . Matrix 𝐻 + is positive definite provided that at least 𝑑 mesh vertices arefixed. If less than 𝑑 points are fixed, rigid body transformations areallowed. 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 𝐻 + 𝑖𝑖 arealways positive definite. Refer to Appendix A for further details onthe finite element assembly procedure. 𝜀 𝑘 In this section we provide a strategy for the choice of the regu-larization parameter 𝜀 𝑘 at each iteration. Namely, we prove thatan idealized minimization algorithm reaches the admissible setmin det 𝐽 > Let us suppose that the admissible set is not empty,namely there exists a mesh 𝑈 ∗ satisfying 𝐹 ( 𝑈 ∗ , ) < +∞ . We alsosuppose that we have a minimization algorithm satisfying one of thefollowing efficiency conditions for some < 𝜎 < : • Either the essential descent condition holds 𝐹 ( 𝑈 𝑘 + , 𝜀 𝑘 ) ≤ ( − 𝜎 ) 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 ) , (10) • or the vector 𝑈 𝑘 satisfies the quasi-minimality condition: min 𝑈 𝐹 ( 𝑈 , 𝜀 𝑘 ) > ( − 𝜎 ) 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 ) . (11) Then the admissible set is reachable by solving a finite number ofminimization problems in 𝑈 with 𝜀 𝑘 fixed for each problem. In otherwords, under a proper choice of the regularization parameter sequence 𝜀 𝑘 , 𝑘 = . . . 𝐾 , we obtain 𝐹 ( 𝑈 𝐾 , ) < +∞ . Proof. The main idea is to expose an explicit way to build adecreasing sequence { 𝜀 𝑘 } ∞ 𝑘 = such that the sequence { 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 )} ∞ 𝑘 = is bounded from above. Then we can prove that the admissible setis reachable in a finite number of steps by a simple reductio adabsurdum argument.First of all, the function 𝐹 ( 𝑈 , 𝜀 ) can be rewritten as follows 𝐹 ( 𝑈 , 𝜀 ) = ∑︁ 𝑖 𝛼 𝑖 𝜓 𝑖 ( 𝑈 ) 𝜒 ( 𝐷 𝑖 , 𝜀 ) , (12)where 𝐷 𝑖 = 𝐷 𝑖 ( 𝑈 ) denotes the Jacobian determinant for 𝑖 -th simplexof the mesh ( 𝐷 𝑖 = det 𝐽 𝑖 ), and 𝛼 𝑖 > 𝜓 𝑖 ( 𝑈 , 𝜀 ) : = 𝜒 ( 𝐷 𝑖 , 𝜀 ) − 𝑑 tr 𝐽 ⊤ 𝑖 𝐽 𝑖 + 𝜆 ( 𝐷 𝑖 + ) defined according to Eq.(3) are positive and bounded from below as 𝛼 𝑖 𝜓 𝑖 ( 𝑈 , 𝜀 ) ≥ 𝜆 min 𝑖 𝛼 𝑖 . Note also that 𝜓 𝑖 ( 𝑈 , 𝜀 ) are increasing functions of 𝜀 .Our goal is to build a decreasing sequence { 𝜀 𝑘 } ∞ 𝑘 = such thatthe sequence { 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 )} ∞ 𝑘 = is bounded from above. We split theconstruction into two parts: first we suppose that at some iteration 𝑘 the essential condition (10) is satisfied, and then we explore thecase (11). Suppose that the condition (10) holds at iteration 𝑘 . In order toguarantee that the function does not increase, it suffices to establishthe following inequality: ( − 𝜎 ) 𝐹 ( 𝑈 𝑘 + , 𝜀 𝑘 + ) ≤ 𝐹 ( 𝑈 𝑘 + , 𝜀 𝑘 ) . (13)By noting that 𝜓 𝑖 ( 𝑈 𝑘 + , 𝜀 𝑘 + ) ≤ 𝜓 𝑖 ( 𝑈 𝑘 + , 𝜀 𝑘 ) , Ineq. (13) is impliedif the following condition holds: ∀ 𝑖 : ( − 𝜎 ) 𝜒 ( 𝐷 𝑘 + 𝑖 , 𝜀 𝑘 ) ≤ 𝜒 ( 𝐷 𝑘 + 𝑖 , 𝜀 𝑘 + ) (14)where 𝐷 𝑘 + 𝑖 : = 𝐷 𝑖 ( 𝑈 𝑘 + ) denotes the Jacobian determinant of sim-plex 𝑖 at iteration 𝑘 + 𝜀 𝑘 + such that Ineq. (14)is satisfied. To do so, first note that 𝜒 is convex with respect to thesecond argument, hence 𝜒 ( 𝐷 𝑘 + 𝑖 , 𝜀 𝑘 + ) − 𝜒 ( 𝐷 𝑘 + 𝑖 , 𝜀 𝑘 ) ≥ ( 𝜀 𝑘 + − 𝜀 𝑘 ) 𝜕𝜕𝜀 𝜒 ( 𝐷 𝑘 + 𝑖 , 𝜀 𝑘 ) . Thus, since 𝜒 and 𝜕𝜕𝜀 𝜒 ( 𝐷, 𝜀 ) are both positive, Ineq. (14) holds pro-vided that ∀ 𝑖 : 𝜀 𝑘 + ≥ 𝜀 𝑘 − 𝜎 𝜒 ( 𝐷 𝑘 + 𝑖 , 𝜀 𝑘 ) 𝜕𝜕𝜀 𝜒 ( 𝐷 𝑘 + 𝑖 , 𝜀 𝑘 ) (15)In its turn, since 𝜒 ( 𝐷,𝜀 ) 𝜕𝜕𝜀 𝜒 ( 𝐷,𝜀 ) is an increasing function of 𝐷 , condi-tion (15) can be simplified to 𝜀 𝑘 + ≥ 𝜀 𝑘 − 𝜎 𝜒 ( 𝐷 𝑘 + − , 𝜀 𝑘 ) 𝜕𝜕𝜀 𝜒 ( 𝐷 𝑘 + − , 𝜀 𝑘 ) , where 𝐷 𝑘 + − : = min 𝑖 𝐷 𝑘 + 𝑖 is the minimum value of the Jacobiandeterminant over all simplices.Hence, if 𝑈 𝑘 + is an approximate solution of the minimizationproblem arg min 𝑈 𝐹 ( 𝑈 , 𝜀 𝑘 ) with fixed parameter 𝜀 𝑘 , we can use thefollowing update rule for 𝜀 𝑘 + : 𝜀 𝑘 + = 𝜀 𝑘 − 𝜎 𝜒 ( 𝐷 𝑘 + − ,𝜀 𝑘 ) 𝜕𝜕𝜀 𝜒 ( 𝐷 𝑘 + − ,𝜀 𝑘 ) , if 𝐷 𝑘 + − < ,𝜀 𝑘 − 𝜎 𝜒 ( ,𝜀 𝑘 ) 𝜕𝜕𝜀 𝜒 ( ,𝜀 𝑘 ) , if 𝐷 𝑘 + − ≥ , (16)This update rule guarantees that Ineq. (13) is satisfied; coupled withthe assumption (10) of the theorem, this implies the required non-growth property of the function values sequence: 𝐹 ( 𝑈 𝑘 + , 𝜀 𝑘 + ) ≤ 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 ) . (17) Consider now the case where condition (11) holds at iteration 𝑘 . Notethat 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 𝐹 ( 𝑈 𝑘 + , 𝜀 𝑘 + ) ≤ ( − 𝜎 ) 𝐹 ( 𝑈 𝑘 + , 𝜀 𝑘 ) ≤ ( − 𝜎 ) 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 ) << ( − 𝜎 ) min 𝑈 𝐹 ( 𝑈 , 𝜀 𝑘 ) < ( − 𝜎 ) min 𝑈 𝐹 ( 𝑈 , ) . Here the last inequality provides a global bound on the functionvalues sequence, and it is based on the observation 𝜕𝜕𝜀 𝜒 ( 𝐷, 𝜀 ) > Technical Report, 2021.
To sum up, we have shown a way to build a sequence { 𝜀 𝑘 } ∞ 𝑘 = suchthat the sequence { 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 )} ∞ 𝑘 = is bounded from above. Now letus prove that the update rule (16) allows to reach the admissible setin 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 infinitedecreasing sequence { 𝜀 𝑘 } ∞ 𝑘 = built using the update rule (16), i.e. 𝐷 𝑘 + − < ∀ 𝑘 ≥ 𝜕𝜕𝜀 𝜒 ( 𝐷, 𝜀 ) ≤ and { 𝜀 𝑘 } ∞ 𝑘 = form adecreasing sequence of non-negative values, we can rewrite theupdate rule (16) for some 𝐾 > 𝜀 − 𝜀 𝐾 ≥ 𝜎 𝐾 − ∑︁ 𝑘 = 𝜒 ( 𝐷 𝑘 + − , 𝜀 𝑘 ) ≥ 𝜎𝐾 min ≤ 𝑘 < 𝐾 𝜒 ( 𝐷 𝑘 + − , 𝜀 𝑘 ) , with an immediate consequence that for an arbitrarily large 𝐾 wehave the following inequality:max ≤ 𝑘 < 𝐾 𝜒 ( 𝐷 𝑘 + − , 𝜀 𝑘 ) ≥ 𝜎𝐾𝜀 . Since all terms 𝛼 𝑖 𝜓 𝑖 ( 𝑈 ) in (12) are bounded from below, the resultingestimate contradicts the boundedness of 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 ) , thus concludingour proof. □ Remark 1.
An important corollary of Th. 1 is that, provided thatthe admissible set is not empty, there exists an iteration 𝐾 < ∞ suchthat the global minimum of the function 𝐹 ( 𝑈 , 𝜀 𝐾 ) belongs to the ad-missible set. The proof is rather obvious: suppose we have an idealizedminimizer such that 𝑈 𝑘 + = arg min 𝑈 𝐹 ( 𝑈 , 𝜀 𝑘 ) . This minimizer alwayssatisfies the conditions of Th. 1, therefore it can untangle the mesh ina finite number of steps. Remark 2.
For our choice of the regularization function 𝜒 , theupdate rule (16) can be instatiated as follows: 𝜀 𝑘 + = (cid:18) − 𝜎 √ ( 𝐷 𝑘 + − ) +( 𝜀 𝑘 ) | 𝐷 𝑘 + − |+ √ ( 𝐷 𝑘 + − ) +( 𝜀 𝑘 ) (cid:19) 𝜀 𝑘 , if 𝐷 𝑘 + − < , ( − 𝜎 ) 𝜀 𝑘 , if 𝐷 𝑘 + − ≥ , (18) In practice the global estimate 𝜎 is not known in advance, and theoptimization routine may be far from the ideal. For each minimizationstep we compute the local descent coefficient 𝜎 𝑘 : 𝜎 𝑘 : = − 𝐹 ( 𝑈 𝑘 + , 𝜀 𝑘 ) 𝐹 ( 𝑈 𝑘 , 𝜀 𝑘 ) . When 𝜎 𝑘 ≥ 𝜎 one can use the update rule (18) using the local value 𝜎 𝑘 guaranteeing that Ineq. (17) holds. In the case 𝜎 𝑘 < 𝜎 one shouldcheck that condition (11) holds for prescribed 𝜎 . If positive, we canassign 𝜎 𝑘 = 𝜎 and use update rule (18) . If one cannot assure (11) , itmeans that minimization procedure for 𝜀 𝑘 failed and theorem cannotbe applied (it does not mean that Alg. 1 will not reach the admissibleset!). In numerical experiments we use value 𝜎 = . Producing maps without reverted elements is a major challengein geometry processing. Inspired by untangling solutions in com-putational physics, our solution outperforms the state of the artboth in terms of robustness and distortion. It is easy to use since weprovide a simple implementation that is free of commercial prod-uct dependency (compiler, library, etc.). Moreover, the energy isestimated independently on each triangle / tetrahedra, making it agood candidate to be adapted to more difficult settings includingfree boundary and global parameterization.
REFERENCES
Noam Aigerman and Yaron Lipman. 2013. Injective and Bounded Distortion Mappingsin 3D.
ACM Trans. Graph.
32, 4, Article 106 (July 2013), 14 pages. https://doi.org/10.1145/2461912.2461931John M Ball. 1976. Convexity conditions and existence theorems in nonlinear elasticity.
Archive for rational mechanics and Analysis
63, 4 (1976), 337–403.J. M. Ball. 1981. Global invertibility of Sobolev functions and the interpenetration ofmatter.
Proceedings of the Royal Society of Edinburgh: Section A Mathematics
88, 3-4(1981), 315–328. https://doi.org/10.1017/S030821050002014XDavid 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. https://doi.org/10.1145/2461912.2462014J.U Brackbill and J.S Saltzman. 1982. Adaptive zoning for singular problems in twodimensions.
J. Comput. Phys.
46, 3 (1982), 342 – 368. https://doi.org/10.1016/0021-9991(82)90020-1Marcel Campen, Cláudio T. Silva, and Denis Zorin. 2016. Bijective Maps from SimplicialFoliations.
ACM Trans. Graph.
35, 4, Article 74 (July 2016), 15 pages. https://doi.org/10.1145/2897824.2925890AA Charakhch’yan and SA Ivanenko. 1997. A variational form of the Winslow gridgenerator.
J. Comput. Phys.
Memo, LawrenceLivermore National Lab
Finite Elements in Analysis and Design
ACM Trans. Graph.
39, 4, Article 120 (July 2020), 17 pages. https://doi.org/10.1145/3386569.3392484Matthias Eck, Tony DeRose, Tom Duchamp, Hugues Hoppe, Michael Lounsbery, andWerner Stuetzle. 1995. Multiresolution Analysis of Arbitrary Meshes. In
Proceedingsof the 22nd Annual Conference on Computer Graphics and Interactive Techniques (SIG-GRAPH ’95) . Association for Computing Machinery, New York, NY, USA, 173–182.https://doi.org/10.1145/218380.218440José Marıa Escobar, Eduardo Rodrıguez, Rafael Montenegro, Gustavo Montero, andJosé Marıa González-Yuste. 2003. Simultaneous untangling and smoothing of tetra-hedral meshes.
Computer Methods in Applied Mechanics and Engineering
Comput. Aided Geom. Des.
14, 3 (April 1997), 231–250. https://doi.org/10.1016/S0167-8396(96)00031-3Lori A Freitag and Paul Plassmann. 2000. Local optimization-based simplicial meshuntangling and improvement.
Internat. J. Numer. Methods Engrg.
49, 1-2 (2000),109–125.VA Garanzha. 2000. The barrier method for constructing quasi-isometric grids.
Com-putational Mathematics and Mathematical Physics
40 (2000), 1617–1637.V.A. Garanzha, L.N. Kudryavtseva, and S.V. Utyuzhnikov. 2014. Variational method foruntangling and optimization of spatial meshes.
J. Comput. Appl. Math.
269 (2014),24 – 41. https://doi.org/10.1016/j.cam.2014.03.006J. Gregson, A. Sheffer, and E. Zhang. 2011. All-Hex Mesh Generation via VolumetricPolyCube Deformation.
Computer Graphics Forum (Special Issue of Symposium onGeometry Processing 2011)
30, 5 (2011).Kai Hormann, Konrad Polthier, and Alia Sheffer. 2008. Mesh Parameterization: Theoryand Practice. In
ACM SIGGRAPH ASIA 2008 Courses (Singapore) (SIGGRAPH Asia’08) . Association for Computing Machinery, New York, NY, USA, Article 12, 87 pages.https://doi.org/10.1145/1508044.1508091SA Ivanenko. 1988. Construction of nondegenerate grids.
Zh. Vychisl. Mat. Mat. Fiz
Computer methods in applied mechanics andengineering
66, 3 (1988), 323–338.Technical Report, 2021. oldover-free maps in 50 lines of code • 11
Fig. 13. On each simplex the map (cid:174) 𝑢 ( (cid:174) 𝑥 ) is affine and is entirely defined bythe position of the vertices of the domain simplex { (cid:174) 𝑥 𝑖 } and its image {(cid:174) 𝑢 𝑖 } . Zhongshi Jiang, Scott Schaefer, and Daniele Panozzo. 2017. Simplicial Complex Aug-mentation Framework for Bijective Maps.
ACM Trans. Graph.
36, 6, Article 186 (Nov.2017), 9 pages. https://doi.org/10.1145/3130800.3130895Zhongshi Jiang, Teseo Schneider, Denis Zorin, and Daniele Panozzo. 2020. BijectiveProjection in a Shell.
ACM Trans. Graph.
39, 6, Article 247 (Nov. 2020), 18 pages.https://doi.org/10.1145/3414685.3417769Patrick M Knupp. 2001. Hexahedral and tetrahedral mesh untangling.
Engineering withComputers
17, 3 (2001), 261–268.Minchen Li, Zachary Ferguson, Teseo Schneider, Timothy Langlois, Denis Zorin, DanielePanozzo, Chenfanfu Jiang, and Danny M. Kaufman. 2020. Incremental PotentialContact: Intersection-and Inversion-Free, Large-Deformation Dynamics.
ACMTrans. Graph.
39, 4, Article 49 (July 2020), 20 pages. https://doi.org/10.1145/3386569.3392425Yaron Lipman. 2012. Bounded Distortion Mapping Spaces for Triangular Meshes.
ACMTrans. Graph.
31, 4, Article 108 (July 2012), 13 pages. https://doi.org/10.1145/2185520.2185604Matthias Nieser, Ulrich Reitebuch, and Konrad Polthier. 2011. CubeCover - Parameteriza-tion of 3D Volumes.
Computer Graphics Forum (2011). https://doi.org/10.1111/j.1467-8659.2011.02014.xMichael Rabinovich, Roi Poranne, Daniele Panozzo, and Olga Sorkine-Hornung. 2017.Scalable Locally Injective Mappings.
ACM Trans. Graph.
36, 2, Article 16 (April2017), 16 pages. https://doi.org/10.1145/2983621Martin Rumpf. 1996. A variational approach to optimal meshes.
Numer. Math.
72, 4(1996), 523–540.Christian Schüller, Ladislav Kavan, Daniele Panozzo, and Olga Sorkine-Hornung. 2013.Locally Injective Mappings.
Computer Graphics Forum (proceedings of Symposiumon Geometry Processing)
32, 5 (2013).Hanxiao Shen, Zhongshi Jiang, Denis Zorin, and Daniele Panozzo. 2019. ProgressiveEmbedding.
ACM Trans. Graph.
38, 4, Article 32 (July 2019), 13 pages. https://doi.org/10.1145/3306346.3323012Jason Smith and Scott Schaefer. 2015. Bijective Parameterization with Free Boundaries.
ACM Trans. Graph.
34, 4, Article 70 (July 2015), 9 pages. https://doi.org/10.1145/2766947Thomas Toulorge, Christophe Geuzaine, Jean-François Remacle, and Jonathan Lam-brechts. 2013. Robust untangling of curvilinear meshes.
J. Comput. Phys.
254 (2013),8–26.W. T. Tutte. 1963. How to Draw a Graph.
Proceedings of the London Mathemat-ical Society s3-13, 1 (01 1963), 743–767. https://doi.org/10.1112/plms/s3-13.1.743 arXiv:https://academic.oup.com/plms/article-pdf/s3-13/1/743/4385170/s3-13-1-743.pdfAlan M Winslow. 1966. Numerical solution of the quasilinear Poisson equation in anonuniform triangle mesh.
Journal of computational physics
1, 2 (1966), 149–172.C. Zhu, R. H. Byrd, and J. Nocedal. 1997. L-BFGS-B: Algorithm 778: L-BFGS-B, FORTRANroutines for large scale bound constrained optimization.
ACM Trans. Math. Software
23, 4 (1997), 550–560.
A COMPREHENSIVE DESIGN FORMULÆ
Given a map (cid:174) 𝑢 , let us denote by (cid:174) 𝑎 𝑖 , 𝑖 = , ( , ) the tangent basis,i.e. vectors forming the columns of the Jacobian matrix 𝐽 . For ex-ample, in 2D we have (cid:174) 𝑎 : = (cid:16) 𝜕 𝑢𝜕 𝑥 𝜕 𝑣𝜕 𝑥 (cid:17) ⊤ and (cid:174) 𝑎 : = (cid:16) 𝜕 𝑢𝜕 𝑦 𝜕 𝑣𝜕 𝑦 (cid:17) ⊤ .Let us denote by (cid:174) 𝑏 𝑖 the dual basis, i.e. vectors chosen in the waythat (cid:174) 𝑎 ⊤ 𝑖 (cid:174) 𝑏 𝑗 = 𝛿 𝑖 𝑗 det 𝐽 for all indices 𝑖, 𝑗 . In particular, for the 2Dsettings the dual basis can be written as (cid:174) 𝑏 : = (cid:16) 𝜕 𝑣𝜕 𝑦 − 𝜕 𝑢𝜕 𝑦 (cid:17) ⊤ and (cid:174) 𝑏 : = (cid:16) − 𝜕 𝑣𝜕 𝑥 𝜕 𝑢𝜕 𝑥 (cid:17) ⊤ . In the 3D case (cid:174) 𝑏 𝑘 = (cid:174) 𝑎 𝑖 × (cid:174) 𝑎 𝑗 , where 𝑖, 𝑗, 𝑘 iscyclic permutation from 1 , ,
3. It is a handy choice of variables, inparticular, tr 𝐽 ⊤ 𝐽 = (cid:205) 𝑖 | (cid:174) 𝑎 𝑖 | and 𝜕 det 𝐽𝜕 (cid:174) 𝑎 𝑖 = (cid:174) 𝑏 𝑖 . For further simplifi-cation of notations we will use 𝜒 for 𝜒 ( 𝐷, 𝜀 ) , 𝜒 ′ for 𝜕 𝜒 ( 𝐷,𝜀 ) 𝜕 𝐷 and 𝑎 ⊤ = ( (cid:174) 𝑎 ⊤ . . . (cid:174) 𝑎 ⊤ 𝑑 ) , 𝑏 ⊤ = ( (cid:174) 𝑏 ⊤ . . . (cid:174) 𝑏 ⊤ 𝑑 ) . A.1 Gradient
In order to derive expressions for the gradient and the Hessianmatrix of 𝐹 , we write down explicitly the Jacobian matrix 𝐽 for theaffine map of a simplex 𝑇 with vertices (cid:174) 𝑢 , (cid:174) 𝑢 , . . . , (cid:174) 𝑢 𝑑 : 𝐽 = ( (cid:174) 𝑎 . . . (cid:174) 𝑎 𝑑 ) = ( (cid:174) 𝑢 − (cid:174) 𝑢 (cid:174) 𝑢 − (cid:174) 𝑢 . . . (cid:174) 𝑢 𝑑 − (cid:174) 𝑢 ) 𝑆 − == ( (cid:174) 𝑢 . . . (cid:174) 𝑢 𝑑 ) 𝑍, where 𝑆 : = ( (cid:174) 𝑥 − (cid:174) 𝑥 (cid:174) 𝑥 − (cid:174) 𝑥 . . . (cid:174) 𝑥 𝑑 − (cid:174) 𝑥 ) , det 𝑆 > (cid:174) 𝑥 𝑖 are vertices of “ideal” or “target” shape forthe image of the simplex 𝑇 , and 𝑍 is a ( 𝑑 + ) × 𝑑 matrix defined as 𝑍 : = { 𝑧 𝑖 𝑗 } : = (cid:18) − . . . − 𝐼 (cid:19) 𝑆 − Since the Jacobian matrix is a linear function of (cid:174) 𝑢 𝑖 , we have 𝜕 (cid:174) 𝑎 𝑖 𝜕 (cid:174) 𝑢 ⊤ 𝑗 = 𝑧 𝑗𝑖 𝐼, 𝑖 = , . . . , 𝑑, 𝑗 = , . . . , 𝑑. The additive contribution to gradient of 𝐹 from the simplex 𝑇 canbe written using correspondence of local indices 0 − 𝑑 and globalindices 𝑔 − 𝑔 𝑑 in the list of vertices: (∇ 𝐹 ) 𝑔 𝑗 + = det 𝑆𝑑 ! 𝑑 ∑︁ 𝑖 = 𝜕 (cid:174) 𝑎 ⊤ 𝑖 𝜕 (cid:174) 𝑢 𝑗 𝜕 𝜙𝜕 (cid:174) 𝑎 𝑖 == det 𝑆𝑑 ! 𝑑 ∑︁ 𝑖 = 𝑧 𝑗𝑖 𝜕 𝜙𝜕 (cid:174) 𝑎 𝑖 , 𝑗 = , . . . , 𝑑, where function 𝜙 ( 𝐽 ) is defined in § 4.1. Let us provide an explicitexpression for 𝜕 𝜙𝜕 (cid:174) 𝑎 𝑖 : 𝜕 𝜙𝜕 (cid:174) 𝑎 𝑖 = 𝜒 𝑑 (cid:174) 𝑎 𝑖 − 𝜒 (cid:18) 𝑑 𝑓 𝜀 𝜒 ′ − 𝜆 det 𝐽 + 𝜆𝑔 𝜀 𝜒 ′ (cid:19) (cid:174) 𝑏 𝑖 A.2 Hessian
The blocks of the nonnegative definite part of Hessian matrix of 𝐹 can be updated using the following general formula 𝐻 + 𝑔 𝑗 𝑔 𝑖 + = det 𝑆𝑑 ! ∑︁ 𝑚,𝑙 𝜕 (cid:174) 𝑎 ⊤ 𝑚 𝜕 (cid:174) 𝑢 𝑗 𝑀 + 𝑚𝑙 𝜕 (cid:174) 𝑎 𝑙 𝜕 (cid:174) 𝑢 ⊤ 𝑖 , Technical Report, 2021. where 𝑀 + 𝑚𝑙 denotes a 𝑑 × 𝑑 block of 𝑑 × 𝑑 positive definite matrix 𝑀 + defined in Eq. (9). Let us provide an explicit expression for thematrix: 𝑀 + = (cid:0) 𝐼 𝑏 (cid:1) (cid:32) 𝜕 Φ 𝜕 𝑎 𝜕 𝑎 ⊤ 𝜕 Φ 𝜕 𝑎 𝜕 𝐷𝜕 Φ 𝜕 𝐷 𝜕 𝑎 ⊤ 𝜕 Φ 𝜕 𝐷 (cid:33) (cid:0) 𝐼 𝑏 (cid:1) ⊤ , where 𝜕 Φ 𝜕 𝑎 𝜕 𝑎 ⊤ = 𝜒 𝑑 𝐼𝜕 Φ 𝜕 𝐷 = 𝑑 (cid:18) + 𝑑 (cid:19) | 𝑎 | 𝜒 ′ 𝜒 + 𝑑 + 𝜆 (cid:32) 𝜒 − 𝐷 𝜒 ′ 𝜒 + ( + 𝐷 ) 𝜒 ′ 𝜒 (cid:33) 𝜕 Φ 𝜕 𝑎 𝜕 𝐷 = − 𝑑 𝜒 ′ 𝜒 + 𝑑 𝑎. Obviously the leading 𝑑 × 𝑑 blocks of the matrix 𝐻 + are strictlypositive definite and can be used to build Newton-type minimizationalgorithm. B CONVEXITY OF Φ Recall that the function Φ is defined as follows (§ 4.1): Φ ( 𝑎, 𝐷 ) : = | 𝑎 | 𝑞 𝑑 + 𝜆 𝐷 + 𝑞 , where 𝑎 ∈ R 𝑑 is the (columnwise) flattening of the Jacobian matrix 𝐽 , i.e. the vector composed of the elements of 𝐽 . Lemma 1. ∇∇ ⊤ Φ > ( 𝑑 + ) × ( 𝑑 + ) Hessian matrix of Φ can be written in the 2 × ∇∇ ⊤ Φ = 𝑃 + 𝜆𝑄, where 𝑃 : = (cid:169)(cid:173)(cid:173)(cid:171) 𝑞 𝑑 𝐼 − 𝑑 𝑞 ′ 𝑎𝑞 + 𝑑 − 𝑑 𝑞 ′ 𝑎 ⊤ 𝑞 + 𝑑 𝑑 ( + 𝑑 ) | 𝑎 | 𝑞 ′ 𝑞 + 𝑑 (cid:170)(cid:174)(cid:174)(cid:172) and 𝑄 : = (cid:32) 𝑞 − 𝐷 𝑞 ′ 𝑞 + ( + 𝐷 ) 𝑞 ′ 𝑞 (cid:33) . It is trivial to verify that 𝑄 ≥
0, since 𝑄 is a strictly positivequadratic function of argument 𝐷 . Since the leading blocks of thematrix 𝑃 are positive definite and the Schur complement 𝑃 − 𝑃 𝑃 − 𝑃 = | 𝑎 | 𝑞 + 𝑑 𝑑 (cid:18) − 𝑑 (cid:19) ≥ Φ is established. □ Remark 3.
Note that we have just proved the convexity of thefunction Φ . As an immediate consequence we obtain polyconvexity ofthe functional (1) , because it is a particular case of our functional (4) with 𝜒 = 𝑞 ..