Constrained willmore surfaces

Smooth curves and surfaces can be characterized as minimizers of squared curvature bending energies subject to constraints. In the univariate case with an isometry (length) constraint this leads to classic non-linear splines. For surfaces, isometry is too rigid a constraint and instead one asks for minimizers of the Willmore (squared mean curvature) energy subject to a conformality constraint. We present an efficient algorithm for (conformally) constrained Willmore surfaces using triangle meshes of arbitrary topology with or without boundary. Our conformal class constraint is based on the discrete notion of conformal equivalence of triangle meshes. The resulting non-linear constrained optimization problem can be solved efficiently using the competitive gradient descent method together with appropriate Sobolev metrics. The surfaces can be represented either through point positions or differential coordinates. The latter enable the realization of abstract metric surfaces without an initial immersion. A versatile toolkit for extrinsic conformal geometry processing, suitable for the construction and manipulation of smooth surfaces, results through the inclusion of additional point, area, and volume constraints.


INTRODUCTION
A draft person's physical spline, i.e., a thin elastic rod, has long been a model for smooth curves. Subject to point constraints, the spline minimizes bending, i.e., the squared curvature integrated along the curve, while maintaining its length. This captures both physical (least bending) and aesthetic ("as-round-as-possible") aspects. Top: Tori which minimize the Willmore energy subject to a conformality constraint. The inset parallelipiped visualizes the fundamental domain of each torus. Its "tilt" controls the twist of the torus. Without constraining the conformal class the minimum in all four cases would be achieved at the standard Cliord torus. Boom: A more complex shape with point constraints showing the minimizer of the Willmore energy without (le) and with (right) the conformal class constraint.
In the case of surfaces 5 : " ! R 3 the bending energy is given by the Willmore energy W(5 ) = π " 2 3 , which measures the squared ! 2 -norm of the mean curvature of the surface. As a bending energy it is of great interest in physical modeling ranging from bio-physics [Canham 1970;Helfrich 1973] to computer animation [Bridson et al. 2003;Grinspun et al. 2003]. Based on its "as-round-as-possible" aspect, the Willmore energy also nds application in surface modeling and form nding [Joshi and Séquin 2007;Vaxman et al. 2015Vaxman et al. , 2018Gruber and Aulisa 2020].
Unfortunately without any additional constraints the surface "material" underlying the Willmore energy is too exible. For example, all surfaces of genus 1 have their Willmore minimum at the Cliord torus [Marques and Neves 2014]. Or consider the case of surfaces which are topological spheres, constraining any 4 points always yields the round sphere through those 4 points as minimizer. What is missing is the analog of the isometry constraint from the univariate setting. Requiring the deformation of a surface to be isometric though leaves it too rigid. Instead the natural choice is a conformality constraint [Garsia 1961;Rüedy 1971], making the study of minimizers of the conformally constrained Willmore functional an active research area [Bohle et al. 2008;Rivière 2008;Schätzle 2013;Kuwert and Schätzle 2013;Heller and Ndiaye 2019]. Instead of an arbitrarily exible material the resulting surfaces are now made of an isotropic auxetic (negative Poisson ratio) material [Konaković et al. 2016]. Limiting the exibility of the surface in this manner enforces intrinsic (conformal class) constraints on the extrinsic appearance. In the case of the torus this xes the aspect ratio and skewness of the fundamental domain and the minimizer reects this (Figure 1, top as well as Figure 2). Similarly, for the spherical shape with 4 point constraints the surface deforms in a more intuitive way following the control points ( Figure 19). Figure 1 (bottom) shows a more complex example comparing the minimizer without and with conformal constraint.
More control over the shape is not the only practical benet of the conformal class constraint. Conformality also ensures that the deformation is locally a similarity. In the case of triangle meshes this preserves mesh quality. Degeneration of meshes during optimization is otherwise a common problem [Bobenko and Schröder 2005] (see also Figure 13) and has motivated methods based on conformal deformations [Crane et al. 2011[Crane et al. , 2013 or the addition of conformal penalty forces [Barrett et al. 2016;Gruber and Aulisa 2020]. Neither of these approaches though are entirely satisfactory since they exhibit numerical drift (Figure 12) in the conformal class resp. modify the energy being minimized.
In this paper we take a dierent approach: we employ an exact discrete conformality constraint for extrinsic variational conformal geometry processing of oriented discrete surfaces (simplicial 2-manifolds of arbitrary genus with or without boundary). This constraint can maintain the conformal class of a surface exactly, or be used to achieve a desired target conformal class. It is based on the discrete notion of conformal equivalence of triangle meshes (CETM) [Springborn et al. 2008;Pinkall and Springborn 2021;Gillespie et al. 2021;Campen et al. 2021] and is independent of the precise form of the energy.
For ecient numerical minimization we use gradient ow with an appropriate Sobolev metric [Renka and Neuberger 1995;Eckstein et al. 2007;Schumacher 2017]. Constraints are incorporated via Lagrange multipliers and the associated min/max problem is solved with the recently proposed competitive gradient descent (CGD) method [Schäfer and Anandkumar 2019]. Of interest beyond the conformality constraint are constraints on point positions, surface area, and enclosed volume (see for example Gruber and Aulisa [2020]) and we demonstrate examples of these in various combinations (see Figures 1, 3 and 19). Here we see how the same twists are realized by cylinders with boundary constraints (positions and tangents) and tori (see also Figure 1, top). The conformality of the deformation is visualized here with a checkerboard texture. Fig. 3. A sphere minimizing the Willmore energy subject to constraints on surface area, enclosed volume, and conformal class. Physically such shapes model vesicles, for example red blood cells, which are enclosed by a thin-shell of high membrane and low bending stiness [Canham 1970;Helfrich 1973;Gruber and Aulisa 2020]. Here, the conformality constraint only serves to preserve mesh quality.
While triangle mesh optimization can be parameterized by vertex positions, in particular when starting with a given immersion, we also give a formulation in terms of dierential coordinates [Sorkine 2006;Xu and Zhou 2009]. Aside from making our approach compatible with the many existing geometry processing algorithms which use dierential coordinates, it also enables us to nd conformal immersions for abstract metric surfaces, valuable in mathematical visualization. For triangle meshes such dierential coordinates are piecewise maps, constant per triangle, and we allow them to vary independently [Custers and Vaxman 2020]. Surprisingly, optimization over such triangle elds is well dened even if triangles do not "glue" together and provides a richer search space for the optimization, often allowing more rapid progress towards convergence. Recovery of an actual surface is then ensured through the inclusion of a (linear) integrability constraint. Dierentials as variables furthermore allow us to interpret the Lagrange multipliers associated with the conformal class constraint as quadratic dierentials [Weber et al. 2012] and, in a mechanical analogy, as stress tensors with corresponding forces acting along edges to counteract anisotropic distortion.

DISCRETE FORMULATION
In this section, we x our notation and dene the basic geometric objects needed to formulate variational problems for surfaces.

Discrete Surface
Topology. Throughout we describe the underlying surface by an abstract triangulation M = (V, E, F). The vertices, edges, and faces are denoted by indices 8 2 V, pairs 89 2 E, and triples 89: 2 F, respectively. The dual mesh is denoted by M ⇤ = (V ⇤ , E ⇤ , F ⇤ ) and we refer to dual cells by their corresponding primal cell, that is, we identify V ⇤ F, E ⇤ E, and F ⇤ V with no further comment. We assume that M is a closed oriented surface of arbitrary topology. Treating surfaces with boundary is relegated to Appendix D.
Dierential forms. We use discrete dierential forms from Discrete Exterior Calculus [Desbrun et al. 2008]. 0-, 1-, and 2-forms are assigned to vertices, edges, and facets respectively and represent integral quantities over their respective subdomain. Recall that forms on M ⇤ can be paired with forms on M and integrated: for where the sum is taken over all oriented :-cells f in M and the angle brackets indicate that the product used in the denition is the inner product on R < . This pairing identies the dual space of ⌦ : (M; R < ) with ⌦ 2 : (M ⇤ ; R < ) and the pairing U |V is interpreted as the integral Ø hU^Vi over M.
Geometry. The intrinsic geometry of the surface is specied by positive edge lengths ✓ : E ! R >0 satisfying the triangle inequality in each face. We call such an ✓ a discrete metric on M; this data endows each face with the geometry of a Euclidean triangle. Vertex positions 5 : V ! R 3 , interpolated linearly in the faces, provide a geometric realization of M in R 3 . We say that 5 is an immersion if every vertex star is embedded. An immersion induces a discrete metric given by the edge lengths of the realization ✓ 89 = |5 9 5 8 |.

Geometric Energies
Consider a discrete surface energy E(5 ). A variation of vertex positions5 : V ! R 3 induces an innitesimal change of energẙ We assume that the energy is dierentiable in the sense that all these partial derivatives exist. In this case, we dene the gradient 2-form of E, denoted as grad E(5 ) 2 ⌦ 2 (M ⇤ ; R 3 ), by the requirement hgrad E(5 ) 8 ,5 8 i for all surface variations5 . The components of the gradient 2-form are the partial derivatives with respect to the variables. These must be derived for an implementation. Critical points of E satisfy the Euler-Lagrange equations grad E(5 ) = 0. (1) From now on, we write grad E for the gradient 2-form if the surface 5 is clear from context.
Area. The area of a discrete surface is dened by where 89: is the area of the triangle 89: and L : 3⇤3 is the positive-denite cotan Laplacian. The gradient 2-form is grad A = L5 . The quantity 1 2 L5 is the integrated discrete mean curvature vector, which we informally denote by HN 3 .
Willmore energy. Measuring the ! 2 -norm of the mean curvature gives rise to the Willmore energy, which can be discretized most simply as where ⇤ is the Hodge-star on dual 2-forms, dened as the inverse of the diagonal matrix of vertex dual areas (Appendix D).

CONFORMAL CONSTRAINTS
In this section, we derive the Euler-Lagrange equation for variational problems with a conformal constraint. We begin with preliminaries, recalling the notion of discrete conformal equivalence and the role of length cross ratios and their logarithms in parameterizing the set of discretely conformally equivalent metrics, before examining the overall structure of conformally constrained critical points of a geometric energy. Finally we give concrete expressions for the discrete Euler-Lagrange equations and see that the associated Lagrange multipliers are quadratic dierentials which can be related to earlier work and interpreted as conformal stress tensors. Discrete Conformal Equivalence. In the smooth setting two metrics 6 and6 are conformally equivalent if there exists a function D : " ! R such that6 = 4 2D 6. In the discrete setting we follow [Luo 2004;Springborn et al. 2008;Bobenko et al. 2015;Pinkall and Springborn 2021], calling two discrete metrics ✓,✓ conformally equivalent if there exist vertex scale factors D : V ! R such that for all 89 2 Ẽ This denition has been successful in solving a number of conformal variational problems involving edge lengths [Chow and Luo 2003;Jin et al. 2007;Glickenstein 2011;Springborn 2017;Ge 2018]. The solutions reproduce structures found in the smooth setting, and exhibit convergence under renement [Luo et al. 2020] as shown by theoretical results, as well as by numerical experiments. Importantly, every innitesimal conformal deformation can be realized by innitesimal extrinsic deformations of vertex positions [Lam and Pinkall 2017], except in the case of isothermic surfaces.
Clearly, conformally equivalent metrics have equal length cross ratios. Thus c can be considered as a map dened on the set C(M) of conformal equivalence classes of discrete metrics on M. Moreover, the length cross ratios completely characterize the conformal class: Thus log c = C_ for a linear map C 2 R E⇥E and in terms of A and C the equivalence in Proposition 3.1 is expressed as Consequently, the map C(M) 3 [✓] 7 ! log c 2 im C denes a global chart for C(M).

Constrained Euler-Lagrange Equation
Here we characterize the critical surfaces of conformally constrained variational problems. We begin with a high level view before giving detailed discrete expressions. An innitesimal conformal deformation is a variation5 such that (log c)˚= C_ = 0. By denition, 5 is then critical for a surface energy where ⇡_(5 ) =_ denotes the dierential of _ in terms of surface variations5 . This can also be expressed equivalently as We derive the Euler-Lagrange equations by describing the orthogonal complement of ker(C ⇡_).
Using the fact that the kernel of a linear operator is the orthogonal complement of the image of its adjoint, we have By taking the orthogonal complement of Equation 4 we also have im C ⇤ = ker A ⇤ . Furthermore, if we identify (R E ) ⇤ with R E using the standard inner product then So far, we have shown that 5 is a conformally constrained critical point if and only if there exists @ 2 ker A ⇤ satisfying grad E = ⇡_ ⇤ @. It remains to compute ⇡_ ⇤ @: for a surface variation5 : V ! R 3 we have Consequently, and so we have the Euler-Lagrange equations below.
T 3.2. A discrete immersion of a closed surface 5 : V ! R 3 is critical for a surface energy E under all innitesimal conformal variations if and only if there exists @ : E ! R satisfying ' 89 @ 89 = 0 (6) for all vertices 8 2 V.
Equation 7 can be interpreted as the condition that grad E is orthogonal to the space of discrete conformal immersions, and so we call @ the Lagrange multiplier of the conformally constrained critical surface.

Discrete adratic Dierentials
We call the conformal Lagrange multipliers discrete quadratic differentials, and we write Q(M) : ker A ⇤ . This notion has appeared in the literature [Lam and Pinkall 2016Pinkall , 2017 in the context of discrete minimal surfaces, circle packings, and discrete harmonic functions. In the smooth theory, holomorphic quadratic dierentials also appear as Lagrange multipliers for conformally constrained variational problems. This discrete notion plays an important role in our numerical optimization-below, we interpret @ as a surface stress and then discretize a metric on the space of quadratic dierentials.
A Lagrange multiplier @ can be interpreted as tension along the edges balancing the energy force at the vertices. Dene g 2 Notice that g is a dual 1-form, so 3g is the dual 2-form obtained by totaling the contributions g 89 around a vertex 8. Pairing the Euler-Lagrange equations with an arbitrary surface variation5 gives Each summand hg 89 , 35 89 i measures the (virtual) work done in response to the edge strain 35 89 , and so we can interpret g 89 as an edge-stress. We also see that Euler-Lagrange equations can be expressed as the force balancing 3g = grad E. This shows that @ 89 generates the force g 89 at the vertex 8, and so if @ 89 > 0 then there is a compressive stress on the edge; conversely, if @ 89 < 0 then there is a tensile stress. In this context we think of @ as a force per unit length squared, which diers from the traditional notion of stress as a force per cross-sectional area. Figure 6 visualizes the forces that the Lagrange multiplier induces on the vertices at a conformally constrained Willmore resp. area minimizing surface. Fig. 6. We visualize the conformal tension force at a constrained Willmore surface (Le) and at a constrained minimal surface (Right). Notice that the cylinder is constrained minimal since the induced force counteracts the mean curvature normal, which tries to squeeze the cylinder.
Measuring Innitesimal Conformal Distortion. The interpretation of the Lagrange multipliers as a conformal tension force provides a natural metric on Q(M), which is needed for our numerical optimization. Since @ induces the conformal tension force g (Equation 8) we can dene the ! 2 -norm of @ by the ! 2 -norm of g: where F ⇤ 89 = ✓ 89 /✓ ⇤ 89 are dual edge weights, ✓ 89 are primal edge lengths, and ✓ ⇤ 89 are dual edge lengths. This gives us a diagonal mass matrix that we call I Q 2 R E⇥E . One can dene other ! 2 -metrics on Q(M) by choosing a metric dierent from the one induced by the immersion. Two natural choices are given by the uniformized metric of constant curvature, which discretizes the Weil-Petersson metric, or using the Bergman metric [Crane et al. 2011].

OPTIMIZATION
In this section we develop ows to minimize a geometric energy E(5 ), the principal example of which is gradient descent. To describe a variety of dierent descent methods in a common framework we derive the gradient descent update step as the minimization of a quadratically regularized, local linear approximation of the given E: Here the rst two terms are from the Taylor series expansion of E at 5 in the direction5 and the minimization is over5 . The step size C > 0 is encoded in the quadratic regularization, and describes our trust in the local linear approximation. If the inner product is described by the mass matrix I, i.e., if 5 ,5 =5 > I5 , then the minimizer of the local approximation is the solution to the linear system So we see that the approximation in Equation 10 results in the usual denition of gradient descent. Of course, this gradient descent will not preserve any of the desired constraints, in particular, this ow will change the discrete conformal class.
Sobolev gradients. Choosing a dierent metric (inner product) in the regularizer above denes a dierent notion of gradient. We can then design ecient gradient ows by carefully choosing the metric so that Equation 10 closely approximates E(5 +5 ). For example, replacing the metric by the Hessian of E results in Newton's method; the Hessian, however will generally not be positive-denite and so does not actually dene an inner product. For many geometric energies an appropriate metric is given by the Sobolev metric matching the energy space of E [Schumacher 2017;Yu et al. 2020].
Since the bi-Laplacian energy is the linearization of the Willmore energy, the appropriate inner product is the Sobolev 2 -metric: The existence and analysis of this Willmore gradient ow in the smooth setting is found in [Schumacher 2017], and the numerical results show that using Sobolev gradients makes Willmore minimization remarkably easy (see the comparison in Figure 7).

Projected gradient flow
To preserve the conformal class we consider the projected gradient descent dened by restricting our local approximation to the set of innitesimally conformal variations. Recall that5 : V ! R 3 is innitesimally conformal if C ⇡_(5 ) = 0. Minimizing Equation 10 subject to this linear constraint results in the saddle point system ✓ 1 Note that we use˚as a variable to parameterize@ = C ⇤˚t o avoid the constraint on@ (Equation 6). This is projected gradient descent, and the variables˚are the linearized Lagrange multipliers which ensure that5 is an innitesimal conformal deformation. Since the conformal constraint is nonlinear, taking an explicit step in the direction5 will not exactly preserve the conformal class and so one still needs to project the surface back to the space of conformal immersions [Schumacher 2017].

Competitive gradient flow
As we just saw, evolving the surface by a projected gradient step and a constraint projection is inecient, due to the nonlinear dependence of the discrete conformal class on the vertex positions. Put another way, restricting the evolution to the space of conformal immersions means that the projected gradient descent takes very small steps along the nonlinear subspace in R 3V of conformal immersions.
Rather than choosing the Lagrange multipliers˚to make5 a conformal variation, as in (PGD), we will update both the surface and the Lagrange multipliers according to the competitive gradient descent (CGD) of [Schäfer and Anandkumar 2019]. This approach aims to nd saddle points of the minimax problem on the Lagrangian min 5 max L(5 ,`) : E(5 ) + C ⇤`| _(5 ), which correspond to constrained critical points of E. We extend this evolution to the Riemannian setting and we will see how it naturally generalizes the projected gradient ow even when not on the prescribed constraint manifold.
Analogous to the derivation of gradient descent through a linear approximation of E, we consider a local bilinear approximation of L in terms of the surface (5 ) and the Lagrange multiplier (˚) variations: Notice that the bilinear (third) term is present in both equations and couples them. This particularly simple form allows us to explicitly obtain the updates5 and˚by solving the generalized saddle point system ✓ 1 where X : C ⇡_ is the linearization of the constraints, and where I 5 and I`are the mass matrices in the space of surfaces and Lagrange multipliers, respectively. For the Lagrange multipliers, we use the ! 2 -inner product of the induced conformal tension, dened in Equation 9. Notice that the right hand side is the same as the one that occurs in (PGD) when 5 lies on the constraint manifold and`is orthogonal to the constraint manifold-this shows that the competitive gradient descent provides a minimal modication of projected gradient descent, which can however also explore the space of non-conformal immersions (see Figure 7 for a comparison).
In the same way that the Sobolev gradient descent denes a quasi-Newton method for grad E(5 ) = 0, we see that the Sobolev CGD also denes a quasi-Newton method for the constrained Euler-Lagrange equation from Theorem 3.2. Since this simple modication to PGD makes sure the ow points back to the constraint manifold, it is no longer necessary to explicitly project the surface back to a conformal immersion. This generalization has an additional practical benet over projected methods: one can now prescribe the conformal class to be dierent than that of the initial surface (see Figure 2).

RELAXING INTEGRABILITY
So far we have considered vertex positions as the variables of our constrained optimization. Since dierential coordinates are used extensively in geometry processing [Sorkine 2006;Xu and Zhou 2009], it is worthwhile to see how our conformality constraint and the attendant constrained minimization problems can be parameterized with dierentials. As we will see in this section, doing so reveals the underlying connection between the discrete Lagrange multipliers and their smooth counterparts. It also provides additional computational benets on challenging examples. To make this concrete we introduce triangle elds, which are natural extensions of dierential coordinates beyond just parameterizing surface deformations, and study discrete conformal variational problems at the level of dierential coordinates (see Figure 8).

Triangle Fields and the Lens Complex
The dierential 35 of a piecewise linear immersion 5 : M ! R 3 is a piecewise constant R 3 -valued 1-form on M, i.e. In particular, also the lengths of corresponding halfedges may dier. Therefore it is convenient to consider all objects as living on a modied cell complex, whose discrete dierential forms naturally contain the discrete forms of M as well as the piecewise constant 1-forms ⌦ 1 F (M; R 3 ). The idea is closely related to the halfedge forms introduced in [Custers and Vaxman 2020].
Lens Complex. The lens complex e M of M is obtained by gluing digons in between opposite halfedges in M (Figure 10). Hence we have the following immediate correspondences, Whenever clear from context we label the cells of e M by the corresponding elements in M-for example, when working over the lens complex we choose to write V to indicate the vertex set of e M. Under these identications we immediately get the equalities: Moreover, note that each edge in e E comes with a canonical orientation given by the orientation of the corresponding half edge in H. Hence ⌦ 1 ( e M; R 3 ) R 3H -with no further conditions imposed. Among the elements of ⌦ 1 ( e M; R 3 ) the 1-forms l on M are distinguished by the condition that l 89 = l 98 for all 89 2 H, i.e. l is closed on digons: A similar condition characterizes the piecewise constant 1-forms among the 1-forms on e M: Since the edge vectors of a triangle sum up to zero, every l 2 ⌦ 1 F (M; R 3 ) denes a 1-form on e M which is closed on triangles F. Conversely, every 1-form l 2 ⌦ 1 ( e M; R 3 ) which is closed on triangles F denes a piecewise constant 1-form on M. Hence TM df ω Fig. 9. R 3 -valued discrete piecewise constant 1-forms identify intrinsic tangent spaces (shown as per-triangle copies of the 2D plane, le) with subspaces of R 3 spanned by the embedded triangles (middle, right). Discrete 1-forms 3 5 that arise as dierentials of R 3 -valued zero forms 5 "stitch up" to discrete immersions (middle). However, using this tangent space identification we can define geometric quantities also for more general, not necessarily integrable candidate dierentials l (triangle fields, right).
! 2 inner product. To express orthogonality between forms as it appears, for example, in the Hodge decomposition we will use the ! 2 inner product on ⌦ 1 F (M; R 3 ) l,l = π " hl^⇤li.
Here star denotes the Hodge star with respect to the triangle eld metric and the integral of the piecewise constant 2-form hl^⇤li is given as a sum of its integrals over triangles. Sincel 2 ⌦ 1 F (M; R 3 ) is constant (and hence closed) on each triangle 89:, its integral over 89: can be computed by the usual cotangent formula: π 89: hl 89:^⇤l89: i = 1 2 cot \ 89 |l 89 | 2 + cot \ 9: |l 9: | 2 + cot \ :8 |l :8 | 2 , where \ 89 , \ 9: , \ :8 2 (0, c/2) denote the interior angles in 89: opposite the halfedges 89, 9:, :8. In particular, (⇤U) 89 : where on the right hand sidel is regarded as a 1-form 89 7 !l 89 on the lens complex. Integrability. A triangle eld l is locally the dierential of a surface if it is closed. Since l is already closed on the triangles, the only obstruction to local integrability is whether l is closed on the lenses. For a lens 89 2 E ⇢ e F (3l) 89 = (l 89 + l 98 ).
If l is closed on the lenses then l denes an R 3 -valued discrete 1form on M. The DEC Hodge decomposition shows that the remaining obstruction to global integrability is the harmonic component of l.
We conclude that l is the dierential of a surface if and only if 8=1 is a basis of harmonic 1-forms on M.

Conformal Lagrange Multipliers
Here we derive the conformal Lagrange multipliers for triangle eld energies. We explain how they can be interpreted as conformal stress tensors (quadratic dierentials) and how this discretization is related to other notions in geometry processing.
Length Multi-Ratios. A triangle eld l does not dene a metric on M but induces edge lengths ✓ on the lens complex e M which satisfy the triangle inequality in each face of M with no compatibility condition over the lenses-we call such a collection of halfedge lengths a triangle eld metric on M.
The notion of conformal equivalence extends to triangle eld metrics in the obvious way-two triangle eld metrics ✓ and✓ are called equivalent if there is a function D : V ! R such that ✓ 89 = 4 (D 8 +D 9 )/2 ✓ 89 for all 89 2 H. Following [Bobenko et al. 2016], where discrete conformal equivalence of polyhedral surfaces is studied, we characterize discrete conformal equivalence over the lens-complex by logarithmic length multi-ratios log m : e E ! R log m(✓) 89 = log ✓ ✓ 8; ✓ 9: ✓ 98 ✓ ; 9 ✓ :8 ✓ 89 ◆ = (_ 8; _ ; 9 + _ 9: _ :8 ) + (_ 98 _ 8 9 ) where :, ; are the vertices across the edge 8 9 as before. Now, if we dene-analogously to Section 3-the averaging and multi-ratio matrices, e A 2 R e E⇥V and e C 2 R e E⇥ e E , by then the relationship im e A = ker e C (Lemma A.1) directly implies: P 5.1. Two triangle eld metrics are discretely conformally equivalent if and only if their length multi-ratios are the same.
Euler-Lagrange Equations. Given that im e A = ker e C, the derivation of the Euler-Lagrange equations in the triangle eld setting follows verbatim the arguments of Section 3. Carrying this out one ends up with the following characterization: A triangle eld l is a critical point of a triangle eld energy e E under innitesimal conformal deformations if and only if grad e E = ⇡_ ⇤ @ for @ 2 ker e A ⇤ @ : e E ! R : where grad e E is the ! 2 -gradient of the energy.
To relate this to the surface case we now consider innitesimal conformal triangle eld variations that preserve integrability. By the DEC Hodge decomposition on the lens complex we know that a variationl 2 ⌦ where grad e E is the ! 2 -gradient of the energy.
that acts on deformations of the extrinsic triangle eld. Noticing that this is a tangential R 3 -valued 1-form, we can actually write e g 89: = 1 89: l 89: ( 89: where ( 89: : ) 89: M ! ) 89: M is the conformal stress tensor ( 89: : @ 89 % 89 + @ 9: % 9: + @ :8 % :8 . An innitesimal deformation of the intrinsic geometry of the faces is described by a self-adjoint linear map ⇡ 89: : ) 89: M ! ) 89: M per face 89: 2 F, and as usual in mechanics the stress tensor gives rise to the energy change h( 89: , ⇡ 89: i when the strain ⇡ 89: is attempted. Fig. 11. A conformal stress tensor can be visualized through its principal stress directions. In each face 89: we compute the eigenspace of ( 89: with the largest eigenvalue which defines a line field (aka the horizontal foliation of the quadratic dierential @) and we visualize its integral curves.
Since tr % 89 = 1, the condition that @ is a conformal Lagrange multiplier (Equations 15 and 17) states that the conformal stress tensor is trace-free when averaged around a vertex. To better understand what this means one can compare this to the case when ( 89: is trace-free inside every face: in this case, ( 89: does not resist scaling deformations and only resists deformations that change the corner angles. The trace-free condition averaged to the vertices means that conformal stress tensor does not measure resistance to isotropic scaling, and only resistance to anisotropic distortion. In this context, the notion of anisotropic distortion is described precisely through discrete conformal equivalence.

RESULTS
We begin this section with a discussion of our implementation and its numerical properties before turning to examples of our method in action.
Our algorithm was implemented in C++ and we performed numerical experiments on a MacBookPro class computer (2.8GHz Intel Core i7-7700HQ, 8GB of RAM). To solve the saddle point systems from Section 4 we precomputed symbolic factorizations via PAR-DISO [Alappat et al. 2020;Bollhöfer et al. 2020Bollhöfer et al. , 2019 and then updated the numerical factorizations and applied backsubstitution for each subsequent problem. We also precomputed the numerical factorization of the cotan Laplacian.
The iterations of our optimization are dominated by the solution of a sparse linear system of size 3|V| + |E|. For an overall sense of the runtime we note that for a mesh of 125k triangles, preprocessing takes about 6 seconds; the following iterations take approximately 3 seconds-computing geometric quantities and derivatives takes about 600ms, and constructing and solving the saddle point problem takes an additional 3 seconds.
initial surface [Crane et al. 2013] CETM Constrained Minimizer (ours) quasi-conformal error Fig. 12. Starting with a twisted torus (upper le) our method finds a conformally equivalent Willmore minimizer with very low quasi-conformal distortion (lower le). Using the method of [Crane et al. 2013], based on conformal deformations, numerical dri leads to accumulating quasi-conformal error (lower right). Note in particular the "neck" region where parametric squares are distorted into high aspect ratio rectangles.
A distinction of our approach relative to earlier works is the use of discrete conformal equivalence as the basis for the conformality constraint. This ensures that the conformal class is preserved up to double precision accuracy. Practically this implies that our approach preserves mesh quality at convergence (and in practice throughout the ow). In Figure 13 we visualize the quasiconformal error Q of Willmore minimizers comparing presence resp. absence of the conformality constraint. There is no quasiconformal error when Q = 1. Without any conformality constraint the element quality can be severely degraded since the reparameterization degree of freedom is wholly uncontrolled. This deterioration of the mesh is a well known issue [Bobenko and Schröder 2005;Barrett et al. 2016;Gruber and Aulisa 2020]. Simply discretizing a smooth notion of conformal deformations [Crane et al. 2011[Crane et al. , 2013 is still insucient due to numerical discretization errors accumulating along the ow leading to signicant drift away from the initial conformal class (Figure 12). Such discretization errors depend heavily on the aspect ratios of triangles which our method is quite insensitive to. Figure 14 illustrates this by showing that two quite dierent triangulations of a given shape yield conformally constrained minimizers of the same shape. Fig. 13. Although Willmore minimization with no conformal class constraint (boom row) yields smooth surfaces, there can be severe distortion in the mesh quality. Fixing the discrete conformal class (top row) yields low quasiconformal error as expected. Recall that a conformal deformation is locally a similarity hence tends to preserve element quality.
Our choice of competitive gradient descent as a method for constrained optimization was motivated by comparing it with standard constrained optimization methods (Figure 7). In particular, we considered the projected gradient descent (Section 4.1) and an augmented Lagrangian method [Bertsekas 1996;Nocedal and Wright 2006]-the projected and competitive gradients were dened with respect to both the ! 2 and 2 inner products, and updates were computed via explicit forward steps. Overall, our strategy was both more ecient and more robust to the quality of the initial surface (or triangle eld). As our approach can be seen as a relaxation of the projected gradient descent, we found that the projected gradient descent performed quite well in general. However, it was comparatively slow since it enforced constraint satisfaction every iteration. Not only did it increase the per-iteration computation time, but the evolution would get "stuck" in regions where the projected Willmore gradient is nearly zero- Figure 7 shows an example where our method used non-conformal deformations (highlighted box in gure) to resolve these challenging situations. The Augmented Lagrangian method was the least consistent and it would often get stuck along the way or in undesirable local minima (inset); this method had slow convergence and its performance was very sensitive to the choice of penalty parameters.
Since we also consider a non-convex optimization problem, there are situations where our approach has diculty converging. In particular, initializing the optimization from a nearly degenerate and non-conformal triangle eld results in seemingly random evolution. In practice, we only ran into these problems when trying to dramatically change the prescribed discrete conformal class of a non-integrable triangle eld. Fig. 14. Length cross ratios are a discrete notion of conformal equivalence and remain invariant under conformal change. Consequently our method is insensitive to the quality of the underlying triangulation. We demonstrate this here by starting with two dierent triangulations of a given input shape, resulting in conformally constrained Willmore minimizers of the same shape.

Constrained Willmore Surfaces
Finding optimal shapes has an extensive history in mathematics, with motivation that goes beyond aesthetics. This search often reveals an understanding of deep mathematical concepts and the physical phenomena governed by these ideas. Willmore surfaces arose from the desire to combine minimal surface theory with (Möbius) conformal geometry of the ambient space [Blaschke 1929]. Constrained Willmore surfaces then arise naturally by additionally constraining the intrinsic conformal geometry of the surface-they are the subject of an active mathematical research eld Besides general results about existence and uniqueness of constrained minimizers, the existing theory and examples almost exclusively concern tori. Neither explicit constructions nor numerical visualizations of constrained Willmore surfaces have been available in the general case. By minimizing the Willmore energy in a xed conformal class, we provide a new numerical framework to experimentally study constrained Willmore surfaces of any genus.
sphere. We then expect our Willmore minimizing algorithm to produce round spheres for any genus zero input. In Figure 15 we start with a Morin surface, a fourfold rotationally symmetric surface that appears as the half-way model of a sphere eversion [Francis 1987]. The surface was slightly perturbed to break the symmetry, and given as input to our algorithm, which converges to a round sphere. Figure 3 gives an example of a sphere minimizing the Willmore energy under conformal class, area, and volume constraints.
6.1.2 Genus 1 Surfaces (Tori). Several classes of explicit examples of constrained Willmore tori have been constructed [Pinkall 1985;Heller 2013], but there is still little known about minimizers in conformal classes far away from the square torus (Cliord torus).
Using triangle elds, we nd constrained Willmore tori that are local minimizers for a given conformal class-this also yields a visualization of every abstract conformal type as a surface in space. To nd these constrained minimizers via surface ows it is necessary to begin with an initial immersion of the abstract conformal class, which is just part of the original problem we're trying to solve.
On the other hand, triangle elds of every conformal class are easily constructed by simply taking a triangulation of the fundamental domain. Competitive gradient descent of the triangle eld gives a simple algorithm to compute constrained Willmore tori directly from a description of the conformal class. We build on the spinorial framework of [Chern et al. 2018] to ensure that our minimizers do not have pinch points. Technical and implementation details are given in Appendices C and D. Convergence to a critical point is measured by the residual norm of the Euler-Lagrange equations-we stop our optimization once this error is below 10 6 .
When starting from an initial immersion that exhibits certain symmetries (Figure 17), our algorithm may rst ow towards a critical point that exhibits those same symmetries. We can use this feature to produce symmetric critical points of the Willmore energy that are not minimizers. 6.1.3 Higher Genus Surfaces. Since we consider general conformal variational problems we are able to produce numerical realizations of general constrained Willmore minimizers of higher genus for the very rst time. A gallery of such results, produced by giving various high genus surfaces as input, is shown in Figure 5.

Conformal Surface Modeling
Conformal variational problems provide an extensible surface modeling framework where many important geometric features can be easily controlled. For example, since conformal maps preserve textures our conformal Willmore minimization provides an excellent replacement for standard fairing methods which can distort textures even with a very small amount of smoothing. Below, we consider several illustrative examples of the conformal constraint applied to problems in geometry processing.

Point Constraints.
In the presence of conformality, xing positions of some points provides control of surface geometry. In Figure 1 (bottom), we see that xing vertex positions while minimizing the Willmore energy without the conformal constraint results in a surface that does not resemble the initial geometry at all. On the other hand requiring the deformation to be conformal is enough for few point constraints to encode the overall geometry of the initial surface. Figure 19 helps explain this behavior in the simple case of a ours Moebius subdivision ours-input original Fig. 18. Our Willmore-minimizing flows (right) produce interpolatory surfaces that are locally as sphere-like as possible. On the le, a result of canonical Möbius subdivision [Vaxman et al. 2018], which specifically targets sphere-reproducing interpolation, albeit using a subdivision approach.
sphere. Changing the cross-ratio of any four point constraints leads to constrained Willmore surfaces that are not Möbius equivalent. Instead conformally constrained Willmore surfaces with additional point constraints provide as-spherical-as-possible interpolants of the points. This property has been recognized as desirable in architectural geometry processing. Subdivision schemes based on Möbius geometry were developed to provide approximations of such sphere preserving interpolants [Vaxman et al. 2018]. Figure 18 compares our conformally constrained Willmore surface with Möbius subdivision. Fig. 19. Four points selected on a round sphere fix the edge length cross-ratio on the two triangles formed by them. Moving one of the points asymmetrically changes the induced cross-ratio means, so no Möbius transformation can map between the old and the new point positions. This causes the shape to deviate from the original sphere.
6.2.2 Gradient Domain Surface Processing. Another way to encode the input shape into our framework is to modify the energy functional so that the minimizing surfaces balance smoothness and closeness to the input shape. We used the ! 2 -norm as a measure of similarity, resulting in the modied energy where Y is a smoothness parameter. These energies are commonplace in gradient domain image/surface processing [Mumford and Shah 1989] and computer graphics [Chuang et al. 2016] in the context of image and surface denoising. In contrast to using point constraints, minimizing Equation 21 in a xed conformal class gives only a "soft" constraint on the vertex positions ( Figure 20).

Constrained Area Minimal
Surfaces. The discrete conformal class constraint can be applied to other than curvature energies. A classical example is the area energy, whose minimizers are minimal surfaces modeling the behavior of soap lms. Applying a conformal point constraints Fig. 20. Willmore-regularized surface fairing which jointly optimizes for fidelity to the input, and minimal Willmore energy with the parameter Y controlling the regularizing eect of the Willmore energy.
class constraint in this setting yields a material which can shrink but only isotropically, in eect a type of isotropic auxetic soap lm. Recall that auxetic materials are characterized by a negative Poisson ratio and in the isotropic case best modeled by conformal maps [Konaković et al. 2016]. Figure 6 (cylinder) shows an example of such a conformally constrained minimal surface. The usual shrinkage towards the middle of the cylinder (with xed boundaries) is not present here since it would require anisotropic deformation. In fact the conformal forces counterbalance the mean curvature normal force. Figure 21 shows more complex examples of a cylinder with twist and a sphere with six caps cut o.
We leave exploration of such materials in the context of geometric and physical modeling to future work. Fig. 21. Minimal surfaces which minimize the area energy oen collapse due to their propensity for shrinkage. Once the conformality constraint is added, as done here, any area scaling must be isotropic. The resulting shapes, eectively made of an isotropic auxetic soap film material, are distinctly dierent.

CONCLUSION
We presented conformally constrained variational problems as a new framework for solving problems in mathematical visualization and geometry processing. The conformal constraint is discretized using the notion of discrete conformal equivalence, and we show that this constraint on the intrinsic parameterization has profound eects on the extrinsic shape of a surface. We used this framework to numerically compute constrained Willmore for surfaces of any genus and conformal type for the very rst time. The numerical optimization is based on a manifold reformulation of competitive gradient descent, which we found to be more robust than standard methods. The conformal constraint, in the presence of Willmore optimality, produces interesting surfaces with rich geometric properties that may be desirable for form nding purposes. Our algorithm is suciently general and can be adapted to other geometric energies or user-dened modeling constraints. Investigating the eect of the conformal constraint on other variational problems is an interesting avenue for future work. On the theory side, we look forward to deepen our understanding of the relationship between our results and isothermic triangulated surfaces, discrete holomorphic quadratic dierentials, Dirac operators, and Teichmüller spaces, which would nicely couple smooth and discrete theory.

ACKNOWLEDGMENTS
This work was supported in part by a National Science Foundation graduate research fellowship, the Kortschak Scholars program at Caltech, the Einstein Foundation Berlin, the Deutsche Forschungs Gemeinschaft (SFB/TRR 109 "Discretization in Geometry and Physics"), and the Taiwanese National Center for Theoretical Sciences (Mathematics Division). Additional support was provided by the Information Science and Technology Initiative at Caltech and SideFX Software.

A THE LENS COMPLEX
This appendix details the discrete conformal equivalence over the lens complex e M.

A.1 Discrete Conformal Equivalence
Below, we prove that discrete conformal equivalence over the lens complex is characterized by logarithmic multi-ratios.
L A.1. The averaging and multi-ratio maps dened in Equation 14 satisfy im e A = ker e C.
P. Consider _ 2 im e A. Notice that _ 89 = _ 98 , and so _ descends to a well-dened function over the edges of M. and so _ 2 im e A. This shows that ker e C ✓ im e A and we conclude that im e A = ker e C, as desired. ⇤

B TRIANGLE FIELD EXTENSIONS
To extend geometric variational problems to the space of triangle elds we need to reformulate functions that are dened in terms of vertex positions in terms of the halfedge vectors. More precisely, we say that the triangle eld function e F (l) is an extension of a function F (5 ) that depends on vertex positions 5 if the relationship e F (35 ) = F (5 ) is satised for all discrete surfaces 5 : V ! R 3 . Since 35 is translationally invariant, an extension does not exist unless F (5 ) is translationally invariant. In this section we obtain triangle eld extensions via the Poisson reconstruction; Recall that the Poisson reconstruction, 5 (l), is the projection onto the space of integrable (exact) 1-forms, and it solves the Poisson equation L5 = 3 ⇤ l.
By applying the function to this projection we get an extension: e F (l) : F (5 (l)).
Now we compute the gradient in the space of triangle elds. An innitesimal deformationl induces a deformation5 of the Poisson reconstruction, and so the variation can be computed from the gradient of F : grad e F ,l =e F = grad F |5 .
Since Poisson equation is linear,5 is also the Poisson reconstruction ofl. This implies that a co-closed triangle eld variation causes no change in the reconstruction5 ⌘ 0. By considering co-closedl in the equation above we nd that grad e F is orthogonal to all co-closed 1-forms. By the Hodge decomposition, the gradient is exact and we can nd i : V ! R 3 satisfying grad e F = 3i. The gradient relationship, in terms of i, becomes grad e F ,l = 3i,l = 3 ⇤l |i where we used the self-adjointness of L and the innitesimal Poisson reconstruction L5 = 3 ⇤l. Since5 is arbitrary, we conclude that Li = grad F . Summarizing, the gradient of the extension can be computed by solving a Poisson equation with the gradient 2-form of F (5 ) on the right-hand side: ( Recall that a Poisson equation has a solution if and only if the right hand side has zero mean. The gradient 2-form grad F integrates to zero precisely when F is translationally invariant, and so Equation 22 can always be solved in our setup. We now give some examples of triangle eld extensions.

B.1 Point constraints
Since an integrable triangle eld only describes a surface up to translation, it does not make sense to talk about the vertex positions of the triangle eld 35 2 M. Nevertheless, it is still meaningful to consider the translationally invariant relative point positions, P (5 ) : 5 D 5 E , between distinct vertices D, E 2 V. It is straightforward to check that where X E is the unit impulse located at the vertex E 2 V. We use Equation 22 to nd that grad e P is the harmonic 1-form with source and sink at D and E, respectively-that is, grad e P = 3? where ? : V ! R is the harmonic function with source and sink at D and E, respectively L? = X D X E .
This shows that grad e P (l) = 3? does not depend on l, and so e P is linear. Explicitly, we have e P (l) = l, 3?.
Notice the similarity between these point constraints and the integrability condition in Equation 13: by considering the inner product with a harmonic 1-form dual to a homology generator W we can prescribe the periods Ø W l. On the other hand, considering the inner product with the harmonic 1-form with sources and sinks we can prescribe the dierence of point positions, which can also be expressed as the integral 5 D 5 E = Ø W D E l, where W D E is a curve from E to D. In both cases, inner products with harmonic 1-forms encode the integrals of the triangle eld over curves in the mesh.
Recall that the gradient 2-form is the discrete area vector: The gradient of our extension e V is then 3o 2 ⌦ 1 (M; R 3 ) where o is the solution of the Poisson equation 3 ⇤ 3o = grad V.

C STARTING FROM RANDOM TRIANGLE FIELDS
The integrability conditions of a triangle eld ensure that all faces can be glued together in a consistent way, but they do not ensure that the integrated surface comes from an immersion since there may be pinch points where the vertex star is not embedded (inset). Practically speaking, pinch points very rarely develop when we initialize the optimization with l 0 = 35 , where 35 is the dierential of a smooth surface. On the other hand, when no initial immersion is available (e.g.in Figure 16 and Figure 1) energy minimization of an arbitrary triangle eld l 0 tends to develop many pinch points. In [Chern et al. 2018] the authors introduce discrete spin structures and show that by parameterizing rotations with spinors one can eliminate pinch points through energy minimization. Similarly, by parameterizing triangle elds with spinors we are able to nd immersed minimizers of a much larger class of challenging variational problems-for example, computing constrained Willmore surfaces when only the discrete conformal class is given.

C.1 Discrete Spinors
We quickly review the necessary concepts from [Chern et al. 2018], relegating the technical details to their paper. Given a background metric ✓, a discrete spinor eld is a function k : F ! ( 3 ⇢ H that we interpret as dening a triangle eld by l k = k 3I k inside each face, where I 89: are complex local coordinates. A discrete spin connection is then a choice of square-root of parallel transport g 89 2 C across the edges 89 2 E. The remarkable feature of this framework is that pinch points are encoded in the smoothness of k : an integrable triangle eld l k = 35 has no pinch points if |k 98; g 89 k 89: | < p 2 for all edges 89 2 E. Consequently, the connection Dirichlet energy E g (k ) provides a bound on the number of pinch points.

C.2 Spinorial Triangle Fields
General triangle elds can be parameterized by a triangle eld metric ✓ : e E ! R and face spinors k : F ! ( 3 : where (I 89: ) are isometric local coordinates in each face with respect to the metric ✓.
Conformal Triangle Fields. By xing a reference metric ✓ 0 we can specialize the parameterization above to the space of conformal triangle elds using vertex scale factors D : V ! R: Since this parameterization only produces conformally equivalent triangle elds the discrete conformal constraint can safely be dropped from the optimization when using these coordinates.

C.3 Immersive Regularization
Finally, we get to the main point: when parameterizing triangle elds with spinors we can regularize our variational problem by adding a small multiple of the spin Dirichlet energy This simple regularization almost completely eliminates the emergence of any pinch points, and constrains the search space to immersed surfaces.

D IMPLEMENTATION
In this section we give the details needed to implement the ow from Section 4. We give an overview of the algorithm before describing explicit matrix expressions: The conformal constraint is computed through the logarithmic edge lengths, _ 89 = |5 9 5 8 |, and the cross-ratio map C 2 R E⇥E (Equation 3), which can be expressed as a sum of sparse matrices C 89: 2 R E⇥E dened per face, with non-zero elements as follows: C 89: = 89 9: :8 ! 0 1 1 89 1 0 1 9: