On a linear programming approach to the discrete Willmore boundary value problem and generalizations

On a linear programming approac h to the discrete Willmore b oundary v alue problem and generalizations Thomas Sc ho enemann 1 , Simon Masnou 2 , and Daniel Cremers 3 1 Departmen t of Mathematical Sciences, Lund Universit y , Sw eden 2 Institut Camille Jordan, Universit ´ e Claude-Bernard Lyon 1, CNRS, F rance 3 Departmen t of Computer Science, TU M ¨ unc hen, Germany Abstract. W e consider the problem of finding (possibly non connected) discrete surfaces spanning a finite set of discrete b oundary curv es in the three-dimensional space and minimizing (globally) a discrete energy in volving mean curv ature. Although w e consider a fairly general class of energies, our main focus is on the Willmore energy , i.e. the total squared mean curv ature. Most w orks in the literature hav e been devoted to the appro ximation of a surface ev olving by the Willmore flo w and, in particular, to the ap- pro ximation of the so-called Willmore surfaces, i.e., the critical points of the Willmore energy . Our purpose is to address the delicate task of ap- pro ximating glob al minimizers of the energy under boundary constraints. The main contribution of this work is to translate the nonlinear b ound- ary v alue problem into an integer linear program, using a natural for- m ulation inv olving pairs of elementary triangles c hosen in a pre-sp ecified dictionary and allowing self-in tersection. The reason for suc h strategy is the well-kno wn existence of algorithms that can compute glob al minimizers of a large class of linear optimiza- tion problems, ho wev er at a significant computational and memory cost. The case of integer linear programming is particularly delicate and usual strategies consist in relaxing the in tegral constraint x ∈ { 0 , 1 } in to x ∈ [0 , 1] which is easier to handle. Our work focuses essentially on the connection betw een the in teger linear program and its relaxation. W e prov e that: – One cannot guaran tee the total unimo dularit y of the constraint ma- trix, which is a sufficient condition for the global solution of the relaxed linear program to b e alwa ys in tegral, and therefore to b e a solution of the integer program as well; – F urthermore, there are actually experimental evidences that, in some cases, solving the relaxed problem yields a fractional solution. These facts pro ve that the problem cannot b e tackled with classical linear programming solvers, but only with pure in teger linear solvers. Neverthe- less, due to the very specific structure of the constraint matrix here, w e strongly b eliev e that it should b e p ossible in the future to design ad-hoc in teger solv ers that yield high-definition appro ximations to solutions of sev eral b oundary v alue problems in volving mean curv ature, in particular the Willmore b oundary v alue problem. 2 Thomas Schoenemann, Simon Masnou, and Daniel Cremers 1 In tro duction The Willmore energy of an immersed compact oriented surface f : Σ → R N with b oundary ∂ Σ is defined as W ( f ) = Z Σ | H | 2 dA + Z ∂ Σ κ ds where H is the mean curv ature vector on Σ , κ the geo desic curv ature on ∂ Σ , and dA , ds the induced area and length metrics on Σ , ∂ Σ . The Willmore en- ergy of surfaces with or without b oundary plays an imp ortan t role in geom- etry , elastic membranes theory , strings theory , and image pro cessing. Among the many concrete optimization problems where the Willmore functional ap- p ears, let us mention for instance the mo deling of biological membranes, the design of glasses, and the smo othing of meshed surfaces in computer graphics. The Willmore energy is the sub ject of a long-standing researc h not only due to its relev ance to some ph ysical situations but also due to its fundamental prop- ert y of b eing conformal inv ariant, which mak es it an in teresting substitute to the area functional in conformal geometry . Critical p oin ts of W with resp ect to in terior v ariations are called Willmore surfaces. They are solutions of the Euler- Lagrange equation δ W = 0 whose expression is particularly simple when N = 3: ∆H + 2 H ( H 2 − K ) = 0, b eing K the Gauss curv ature. It is kno wn since Blaschk e and Thomsen [23] that stereographic pro jections of compact minimal surfaces in S 3 ⊂ R 4 are alwa ys Willmore surfaces in R 3 . How ev er, Pink all exhibited in [22] an infinite series of compact embedded Willmore surfaces that are not stereo- graphic pro jections of compact embedded minimal surfaces in S 3 . Y et Kusner conjectured [17] that stereographic pro jections of Lawson’s g -holed tori in S 3 should b e global minimizers of W among all genus g surfaces. This conjecture is still op en, except of course for the case g = 0 where the round sphere is known to b e the unique global minimizer. The existence of smo oth surfaces that minimize the Willmore energy span- ning a given b oundary and a conormal field has b een pro ved by Sch¨ atzle in [27]. F ollowing the notations in [27], we consider a smo oth embedded closed oriented curv e Γ ⊂ R N together with a smooth unit normal field n Γ ∈ N Γ and w e denote as ± Γ and ± n Γ their p ossible orien tations. W e assume that there ex- ist oriented extensions of ± Γ , ± n Γ , that is, there are compact orien ted sur- faces Σ − , Σ + ⊂ R N with b oundary ∂ Σ ± = ± Γ and conormal vector field co Σ ± = ± n Γ on ∂ Σ ± . W e also assume that there exists a b ounded op en set B ⊃ Γ such that the set { Σ ± orien ted extensions of ( Γ , n Γ ) , Σ + connected , Σ + ∪ Σ − ⊂ B , W ( Σ + ∪ Σ − ) < 8 π } is not empt y . The condition on energy ensures that Σ + ∪ Σ − is an em b edding. It follows from [27], Corollary 1.2, that the Willmore b oundary problem as- so ciated with ( Γ , n Γ ) in B has a solution, i.e., there exists a compact, oriented, A linear programming approach to the discrete Willmore b oundary problem 3 connected, smo oth surface Σ ⊂ B with ∂ Σ = Γ , co Σ = n Γ on ∂ Σ , and W ( Σ ) = min { W ( ˜ Σ ) , ˜ Σ smo oth , ˜ Σ ⊂ B , ∂ ˜ Σ = Γ , co ˜ Σ = n Γ on ∂ ˜ Σ } There ha v e b een many con tributions to the n umerical sim ulation of Willmore surfaces in space dimension N = 3. Among them, Hsu, Kusner and Sulliv an ha ve tested exp erimen tally in [16] the v alidit y of Kusner’s conjecture: starting from a triangulated polyhedron in R 3 that is close to a La wson’s surface of gen us g , they let it evolv e by a discrete Willmore flow using Brakk e’s Surface Ev olver [6] and c heck that the solution obtained after conv ergence is W -stable. Recen t up dates that Brakke brough t to its program give no w the p ossibilit y to test the flo w with v arious discrete definitions of the mean curv ature. May er and Simonett [19] introduce a finite difference scheme to approximate axisymmet- ric solutions of the Willmore flow. Rusu [26] and Clarenz et al. [8] use a finite elemen ts appro ximation of the flow to compute the evolution of surfaces with or without b oundary . In b oth works, p osition and mean curv ature vector are tak en as indep enden t v ariables, which is also the case of the con tribution by V erdera et al. [32], where a triangulated surface with a hole in it is restored using the following approac h: by the coarea formula, the Willmore energy (actu- ally a generalization to other curv ature exp onen ts) is replaced with the energy of an implicit and smo oth representation of the surface, and the mean curv ature term is replaced by the divergence of an unkno wn field that aims to represent the normal field. Droske and Rumpf [9] prop ose a finite element approac h to the Willmore flow but replace the standard flo w equation b y its level set for- m ulation. The contribution of Dziuk [10] is tw ofold: it provides a finite elemen t appro ximation to the Willmore flow with or without b oundary conditions that can handle as well em b edded or immersed surfaces (turning the surface problem in to a quasi-planar problem), and a consistency result showing the conv ergence of b oth the discrete surface and the discrete Willmore energy to the contin uous surface and its energy when the approximated surface has enough regularity . Bob enk o and Schr¨ oder [4] use a difference strategy: they introduce a discrete notion of mean curv ature for triangulated surfaces computed from the circles circumscrib ed to eac h triangle that shares with the con tinuous definition a few prop erties, in particular the inv ariance with resp ect to the full M¨ obius group in R 3 . This discrete definition is vertex-based and a discrete flow can b e derived. Based also on several axiomatic constraints but using a finite elements frame- w ork, W ardetzky et al. [33] introduce an edge-based discrete Willmore energy for triangulated surfaces. Olischl¨ ager and Rumpf [21] introduce a tw o step time discretization of the Willmore flow that extends to the Willmore case, at least formally , the discrete time approximation of the mean curv ature motion due to Almgren, T aylor, and W ang [2], and Luckhaus and Sturzenheck er [18]. The strategy consists in using the mean curv ature flow to compute an appro ximation of the mean curv ature and plug it in a time discrete approximation of the Will- more flow. Grzib o vskis and Heintz [14], and Esedoglu et al. [11] discuss how 4th order flows can b e approximated by iterativ e conv olution with suitable kernels and thresholding. 4 Thomas Schoenemann, Simon Masnou, and Daniel Cremers While all the previous approaches yield approximations of critical p oints of the Willmore energy , our motiv ation in this pap er is to approximate global minimizers of the energy . This is an ob viously non trivial task due to the high nonlinearit y and noncon v exity of the energy . Y et, for the simpler area functional, Sulliv an [31] has shown with a calibration argumen t that the task of finding min- imal surfaces can b e turned into a linear problem. Even more, when a discrete solution is seeked among surfaces that are union of faces in a cubic grid parti- tion of R 3 , he prov ed that the minimization of the linear program is equiv alen t to solving a minim um-cost circulation netw ork flow problem, for which efficient co des hav e b een developed by Boyk o v and Kolmogoro v [5] after F ord and F ulk- erson [12]. Sulliv an [31] did not pro vide exp erimen ts in his pap er but this was done recently by Grady [13], with applications to the segmentation of medical images. The linear form ulation that we prop ose here is based on tw o k ey ideas: the concept of surface con tinuation constraints that has b een pioneered by Sulli- v an [31] and Grady [13], and the representation of a triangular surface using pairs of triangles. With this represen tation and a suitable definition of discrete mean curv ature, w e are able to turn into a linear formulation the task of mini- mizing discrete represen tations of any functional of the form W ϕ ( Σ ) = Z Σ ϕ ( x, n, H ) dA among discrete immersed surfaces with b oundary constraints: ∂ Σ = Γ , co ˜ Σ = n Γ on ∂ Σ . In the expression of W ϕ ( Σ ), x denotes the space v ariable, n the normal v ector field on Σ and H the mean curv ature vector. The linear problem w e obtain in volv es integer-v alued unknowns and do es not seem to admit any simple graph- based equiv alent. W e will therefore discuss whether classical strategies for linear optimization can b e used. The paper is organized as follows: in section 2 we discuss both the c hosen rep- resen tation of surfaces and the definition of discrete mean curv ature. In section 3 w e present a first possible approach yielding a quadratic energy . W e present in section 4 our linear formulation and discuss whether it can be tac kled b y classical linear optimization tec hniques. 2 Discrete framew ork 2.1 T riangular meshes from a set of pre-defined triangles The equiv alence shown b y Sulliv an b et ween finding minimal surfaces and solving a flow problem holds true for discrete surfaces defined as a connected set of cell faces in a cellular complex discrete representation of the space. W e will consider here p olyhedral surfaces defined as union of triangles with vertices in (a finite subset of ) the cubic lattice  Z 3 where  = 1 n is the resolution scale. Not all A linear programming approach to the discrete Willmore b oundary problem 5 p ossible triangles are allow ed but only those resp ecting a specified limit on the maximal edge length. W e assume that eac h triangle, as well as each triangle edge, is represented t wice, once for eac h orientation. W e let I denote the collection of orien ted triangles, N = |I | its cardinality , and M the n um b er of oriented triangle edges. The constrained boundary is given as a con tiguous orien ted set of triangle edges. The orientation of the boundary constrains the spanning surfaces since w e will allow only spanning triangles whose orientation is compatible. In this framework, one can represent a triangular mesh as a binary indicator v ector x = { 0 , 1 } N where 1 means that the resp ective triangle is presen t in the mesh, 0 that it is not. Obviously , not all binary indicator vectors can b e asso ciated with a triangular surface since the corresp onding triangles may not b e con tiguous. How ever, as discussed by Grady [13] and, in a slightly differen t setting, by Sulliv an [30,31], it is p ossible to write in a linear form the constraint that only binary vectors that corresp ond to surfaces spanning the given b oundary are considered. W e will see that using the same approach here turns the initial b oundary v alue problem in to a quadratic program. Another formulation will b e necessary to get a linear problem. 2.2 Admissible indicator vectors: a first attempt e4 e3 e2 e1 Fig. 1. Incidence of oriented triangles and edges. e 1 is p ositively inciden t to the orien ted triangle, e 2 and e 3 are negatively incident, and e 4 is not incident to the triangle. T o define the set of admissible indicator vectors, we first consider a relation- ship b et ween oriented triangles and oriented edges which is called incidenc e : a triangle is p ositive incident to an edge if the edge is one of its b orders and the t wo agree in orientation. It is negativ e incident if the edge is one of its b orders, but in the opp osite orientation. Otherwise it is not incident to the edge. F or example, the triangle in Figure 1 is p ositive incident to the edge e 1 , negativ e inciden t to e 2 and e 3 and not inciden t to e 4 . Being defined as ab o ve the set of N oriented triangles and their M oriented edges, w e in tro duce the matrix B = ( b ij ) i ∈{ 1 , ··· ,N } j ∈{ 1 , ··· ,M } whose elemen t b ij giv es accoun t of the incidence b et ween triangle i and edge j . More precisely b ij =      1 if edge i is an edge of triangle j with same orien tation − 1 if edge i is an edge of triangle j with opp osite orientation 0 otherwise 6 Thomas Schoenemann, Simon Masnou, and Daniel Cremers The kno wledge of which edges are present in the set of prescrib ed b oundary segmen ts is expressed as a vector r ∈ {− 1 , 0 , 1 } M with r j =                1 if the orien ted b oundary con tains the edge j with agreeing orien tation − 1 if the orien ted b oundary con tains the edge − j with opp osing orientation 0 otherwise With these notations set up w e can no w describ e the equation system defining that a vector x ∈ { 0 , 1 } N enco des an oriented triangular mesh with the pre- sp ecified oriented b oundary . This system has one equation for each edge. If the edge is not con tained in the given boundary , this equation expresses that, among all triangles indicated by x that contain the edge, there are as many triangles with same orientation as the edge as triangles with opp osite orientation. If the edge is con tained in the boundary with coherent orientation, there m ust be one more p ositiv e incident triangle than negativ e incident. If it is con tained with opp osite orien tation, there is one less positive than negativ e inciden t. Altogether the constrain t for edge j can b e expressed as the linear equation X i b ij x i = r j and the en tire system as B x = r . (1) So far, w e did not incorporate the conormal constrain t. Actually not all conormal constrain ts are p ossible, exactly lik e not all discrete curves can be spanned in our framework but only union of edges of dictionary triangles, i.e. the collection of triangles defined in the previous section that determine the p ossible surfaces. F or the conormal constraint, only the conormal vectors that are tangent to dic- tionary triangles sharing an edge with the b oundary curve are allow ed. Then the conormal constrain t can b e easily plugged into our formulation by simply imp osing the corresp onding triangles to b e part of the surface, see Figure 2, and b y defining accordingly a new b oundary indicator v ector ˜ r . Denoting as J the collection of those additional triangles, the complete con- strain t reads  B x = ˜ r x j = 1 , j ∈ J (2) W e discuss in the next section how discrete mean curv ature can b e ev aluated in this framew ork. 2.3 Discrete mean curv ature on triangular meshes The v arious definitions of discrete mean curv ature that hav e b een prop osed in the literature obviously dep end on the c hosen discrete representations of surfaces. A linear programming approach to the discrete Willmore b oundary problem 7 Fig. 2. The b oundary and conormal constraints can be imp osed b y pre-sp ecifying suit- able triangles to b e part of the surface. Presen ting and discussing all p ossible definitions is out of the scop e of the presen t pap er. The important thing to kno w is that there is no fully consistent definition: the p oin twise conv ergence of mean curv ature cannot b e guaranteed in general but only in sp ecific situations [15,20]. Among the many p ossible definitions, w e will use the edge-based one prop osed by Polthier [24] for it suits with our framew ork. Recalling that, in the smo oth case but also for generalized surfaces lik e v arifolds [29], the first v ariation of the area can b e written in terms of the mean curv ature, the definition due to Polthier of the mean curv ature vector at an in terior edge e of a simplicial surface reads H ( e ) = | e | cos θ e 2 N e (3) where | e | is the edge-length, θ e is the dihedral angle b et ween the tw o triangles adjacen t to e , and N e is the angle bisecting unit normal vector, i.e., the unit v ector collinear to the half sum of the tw o unit vectors normal to the adjacen t triangles (see figure 3). Remark that this formula is a discrete counterpart of the contin uous H = κ 1 + κ 2 dep ending on the principal curv atures, which is used in man y pap ers for simplicity as definition of mean curv ature. When the correct contin uous definition H = 1 2 ( κ 1 + κ 2 ) is used, the formulas ab ov e and hereafter should b e adapted. The justification of formula (3) b y Polthier [24,25] θ e N e e Fig. 3. The edge-based definition of a discrete mean curv ature vector due to Polth- ier [24] depends on the dihedral angle θ e and the angle bisecting unit normal vector N e . is as follows: it is exactly the gradien t at any p oin t m ∈ e of the area of the 8 Thomas Schoenemann, Simon Masnou, and Daniel Cremers t wo triangles T 1 and T 2 adjacen t to e , and this gradient does not dep end on the exact p osition of m . Indeed, one can sub divide T 1 , T 2 in four triangles T 0 i , i ∈ { 1 , · · · , 4 } ha ving m ∈ e as a vertex and such that T 1 = T 0 1 ∪ T 0 2 and T 2 = T 0 3 ∪ T 0 4 . The area of each triangle is half the pro duct of the opp osite edge’s length and the height. Therefore, if e i is the p ositively orien ted edge opp osite to m in the triangle T 0 i and J 1 , J 2 the rotations in the planes of T 1 , T 2 b y π 2 , the area gradients of T 0 i , i ∈ { 1 , · · · , 4 } at m are 1 2 J 1 e 1 , 1 2 J 1 e 2 , 1 2 J 2 e 3 , 1 2 J 2 e 4 . The sum is the total area gradient of T 1 ∪ T 2 at m and equals 1 2 ( J 1 e + J 2 e ), which coincides with the form ula ab o ve. As discussed by W ardetsky et al. using the Galerkin theory of appro ximation, this discrete mean curv ature is an integrated quantit y: it scales as λ when each space dimension is rescaled b y λ . A p oin t wise discrete mean curv ature rescaling as 1 λ is giv en by (see [33]) H pw ( e ) = 3 | e | A e cos θ e 2 N e , where A e denotes the total area of the tw o triangles adjacent to e . The factor 3 comes from the fact that, when the mean curv atures are summed up o ver all edges, the area of each triangle is counted three times, once for each edge. Then a discrete coun terpart of the energy Z Σ ϕ ( H ) dA is given by X edges e A e 3 ϕ ( 3 | e | A e cos θ e 2 N e ) . (4) In particular, the edge-based total squared mean curv ature is X edges e 3 | e | 2 A e (cos θ e 2 ) 2 . (5) 3 A quadratic program for the minimization of the discrete Willmore energy Ultimately w e are aiming at casting the optimization problem in a form that can b e handled by standard linear optimization softw are. Having in mind the frame- w ork described ab o ve where a discrete surface spanning the prescrib ed discrete b oundary is giv en as a collection of orien ted triangles satisfying equation (2) and c hosen among a pre-sp ecified collection of triangles, a somewhat natural direc- tion at first glance seems to b e solving a quadr atic pr o gr am . Lik e in section 2.1, let us indeed denote as ( x i ) the collection of binary v ariables asso ciated to the “dictionary” of triangles ( T i ) and define – e ij the common edge to t wo adjacent triangles T i and T j ; – θ ij the corresp onding dihedral angle; – N ij the angle bisecting unit normal; A linear programming approach to the discrete Willmore b oundary problem 9 – A ij the total area of b oth triangles. Then a con tinuous energy of the form Z Σ ϕ ( x, n, H ) dA can b e discretized as X i,j q ij x i x j (6) with q ij =        1 2 A ij 3 ϕ ( e ij , N ij , 3 | e ij | A ij cos θ ij 2 N ij ) if i 6 = j are adjacent ˜ ϕ ( T i , N i ) if i = j 0 otherwise where ˜ ϕ allows to incorp orate dep endences on each triangle T i ’s p osition and unit normal N i . In particular, the discrete Willmore energy is X i,j q w ij x i x j (7) with q w ij =    3 | e ij | 2 2 A ij (cos θ ij 2 ) 2 if i 6 = j are adjacent 0 otherwise Assuming that the maps ϕ and ˜ ϕ are p ositive-v alued, both energy matrices Q = ( q ij ) and Q w = ( q w ij ) are symmetric matrices in R + N × N , and the mini- mization of either (6) or (7) with b oundary constraints turns to b e the following quadratic program with linear and in tegrality constraints: min x h Q x, x i suc h that B x = r x i = 1 ∀ i ∈ J x ∈ { 0 , 1 } N . W e know of no solution to solve this problem efficiently due to the integralit y constrain t. What is worse, ev en the relaxed problem where one optimizes ov er x ∈ [0 , 1] N is very hard to solve: terms of the form x i x j with i 6 = j are indefinite, so (unless Q has a dominan t diagonal) the ob jective function is a non-con vex one. Moreo ver, a solution to the relaxed problem w ould not b e of practical use: already for the 2D-problem of optimizing curv ature energies o ver curves in the plane, the resp ectiv e quadratic program fav ors fractional solutions. The relax- ation would therefore not b e useful for solving the integer program. How ever, in this case Amini et al. [3] show ed that one can solve a linear program instead. This inspired us for the ma jor contribution of this work: to cast the problem as an in teger linear program. 10 Thomas Schoenemann, Simon Masnou, and Daniel Cremers 4 An in teger linear programming approac h 4.1 Augmen ted indicator vectors The k ey idea of the prop osed in teger linear program is to consider additional indicator vectors. Aside from the indicator v ariables x i for basic triangles, one no w also considers en tries x ij corresp onding to p airs of adjacen t triangles. Such a pair is called quadr angle in the follo wing. W e will denote ˆ x the augmen ted v ector ( x 1 , · · · , x N , · · · , x ij , · · · ) where i 6 = j run ov er all indices of adjacent triangles. The cost function can b e easily written in a linear form with this augmented v ector, i.e. it reads X w k ˆ x k with (see the notations of the previous section) w k =  q ii if ˆ x k = x i q ij if ˆ x k = x ij The ma jor problem to ov ercome is ho w to set up a system of constraints that guaran tees consistency of the augmented vector: the indicator v ariable x ij for the pair of triangles i and j should b e 1 if and only if b oth the v ariables x i and x j are 1. Otherwise it should b e 0. In addition, one again wan ts to optimize only o ver indicator vectors that corresp ond to a triangular mesh. T o enco de this in a linear constraint system, a couple of c hanges are necessary . First of all, we will now hav e a constraint for each pair of triangle and adjacent edge. Secondly , edges are no longer orien ted. Still, the set of pre-sp ecified indices J implies that the orientation of the b order is fixed - w e still require that for eac h edge of the b oundary an adjacent (oriented) triangle is fixed to constrain the conormal information. T o enco de the constrain t system w e in tro duce a modified notion of incidence. W e are no longer in terested in incidence of triangles and edges. Instead we now consider the incidence of both triangles and quadrangles to pairs of triangles and (adjacen t) edges. F or con venience, we define that triangles are p ositiv e inciden t to a pair of edge and triangle, whereas all quadrangles are negativ e incident. W e prop ose an incidence matrix where lines corresp ond to pairs (triangle, edge) and columns to either triangles or quadrangles. The entries of this incidence matrix are either the incidence of a pair (triangle, edge) with a triangle, defined as d ((triangle k , edge e ) , triangle i ) = ( 1 if i = k , e is an edge of triangle i 0 otherwise , or the incidence of a pair (triangle, edge) with a quadrangle, defined as d ((triangle k , edge e ) , quadrangle ij ) = ( − 1 if i = k or j = k and i, j share e 0 otherwise . A linear programming approach to the discrete Willmore b oundary problem 11 The columns of this incidence matrix are of t wo types: either with only 0’s and exactly three 1 (a column corresponding to a triangle T , whose three edges are found at lines ( T , e 1 ), ( T , e 2 ), ( T , e 3 )), or with only 0’s and exactly tw o ( − 1)’s (a column corresp onding to a quadrangle ( T 1 , T 2 ) that matches with lines ( T 1 , e 12 ) and ( T 2 , e 12 )). Again, b oth the conormal constraints and the b oundary edges can be imp osed b y imposing additional triangles indexed b y a collection J of indices. The general constrain t has the form X i d (( x k , e ) , x i ) + X i,j d (( x k , e ) , x ij ) = r 0 ( k,e ) , where the right-hand side dep ends whether the edge e is shared b y tw o triangles of the surface (and even several quadrangles in case of self-intersection), or b e- longs to the new b oundary indicated by the additional triangles. If e is an inner edge, then the sum must b e zero due to our definition of d , otherwise there is an adjacen t triangle, but no adjacent quadrangle, so the right-hand side should b e 1: r 0 ( k,e ) = ( 1 if k ∈ J , e is part of the mo dified b oundary 0 otherwise T o sum up, we get the following integer linear program: min ˆ x h w , ˆ x i (8) suc h that D ˆ x = r 0 ˆ x j = 1 ∀ j ∈ J ˆ x i ∈ { 0 , 1 } ∀ i ∈ { 1 , . . . , ˆ N } where ˆ N is the total num ber of entries in ˆ x , namely all triangles plus all pairs of adjacen t triangles. It is w orth noticing that such formulation allows triangle surfaces with self-in tersection. 4.2 On the linear programming relaxation Solving integer linear programs is an NP-complete problem, see e.g. [28, chapter 18.1]. This implies that, to the noticeable exception of a few particular prob- lems [28], no efficient solutions are known. As a consequence one often resorts to solving the corresp onding linear programming (LP) relaxation, i.e. one drops the in tegrality constraints. In our case this means to solve the problem: min ˆ x h w , ˆ x i (9) suc h that D ˆ x = r 0 ˆ x j = 1 ∀ j ∈ J 0 ≤ ˆ x i ≤ 1 ∀ i ∈ { 1 , . . . , ˆ N } 12 Thomas Schoenemann, Simon Masnou, and Daniel Cremers or, equiv alently , by suitably augmen ting D and r 0 in order to incorp orate the second constrain t ˆ x j = 1, ∀ j ∈ J : min ˆ x h w , ˆ x i (10) suc h that ˆ D ˆ x = ˆ r 0 ≤ ˆ x i ≤ 1 ∀ i ∈ { 1 , . . . , ˆ N } There are v arious algorithms for solving this problem, the most classical b eing the simplex algorithm and several interior p oin t algorithms. Let us now discuss the conditions under which these relaxed solutions are also solutions of the orig- inal integer linear program. Recalling the basics of LP-relaxation [28], the set of admissible solutions P = { ˆ x ∈ R ˆ N , ˆ D ˆ x = ˆ r , 0 ≤ x ≤ 1 } is a polyhedron, i.e. a finite intersection of half-spaces in R ˆ N . A classical result states that minimizing solutions for the linear ob jectiv e functions can b e seek ed among the extremal p oin ts of P only , i.e. its vertices. Denoting P e the integral en velope of P , that is the con vex env elop e of P ∩ Z ˆ N , another classical result states that P has in tegral vertices only (i.e. vertices with in tegral co ordinates) if and only if P = P e Since P = { ˆ x ∈ R ˆ N , ˆ D ˆ x = ˆ r , 0 ≤ ˆ x ≤ 1 } , according to Theorem 19.3 in [28], a sufficient condition for ha ving P = P e is the prop ert y of B b eing totally unimo dular, i.e. any square submatrix has determinant either 0, − 1 or 1. Under this condition, an y extremal p oin t of P that is a solution of min ˆ D ˆ x = ˆ r , ˆ x i ∈ [0 , 1] h w , ˆ x i has integral co ordinates therefore is a solution of the original in teger linear pro- gram min ˆ D ˆ x = ˆ r , ˆ x i ∈{ 0 , 1 } h w , ˆ x i . Theorem 19.3 in [28] mentions an interesting characterization of total unimo du- larit y due to Paul Camion [7]: a matrix is totally unimo dular if, and only if, the sum of the entries of every Eulerian square submatrix (i.e. with even rows and columns) is divisible b y four. Unfortunately , w e can pro ve that, as so on as the triangle space is ric h enough, the incidence matrix ˆ D do es not satisfy Camion’s criterion, therefore is not totally unimo dular, and neither are the matrices for ric her triangles spaces. As a consequence, there are choices of the triangle space for which the p olyhedron P = { ˆ x ∈ R ˆ N , ˆ D ˆ x = ˆ r , 0 ≤ ˆ x ≤ 1 } may ha ve not only in tegral vertices, or more precisely one cannot guaran tee this property thanks to total unimodularity . This is summarized in the follo wing theorem. Theorem 1. The incidenc e matrix asso ciate d with any triangle sp ac e wher e e ach triangle has a lar ge enough numb er of adjac ent neighb ors is not total ly unimo dular. A linear programming approach to the discrete Willmore b oundary problem 13 Pr o of. W e sho w in Figure 4 a configuration and, in T able 1, an asso ciated square submatrix of the incidence matrix. The sum of en tries o ver each line and the sum o ver each column are even, though the total sum of the matrix entries is not divisible b y four. By a result of Camion [7], the incidence matrix is not totally unimo dular which yields the conclusion according to [28][Thm 19.3]. Clearly , an y triangle space for which this configuration can o ccur is also asso ciated to an incidence matrix that is not totally unimo dular. It is worth noticing that the previous theorem do es not imply that the extremal p oin ts of the p olyhedron P are necessarily not all in tegral. It only states that this cannot b e guaranteed as usual by the criterion of total unimo dularit y . W e will discuss in the next section what additional informations ab out integralit y can be obtained from a few exp erimen ts that w e ha ve done using classical solv ers for addressing the relaxed linear problem. ... 24 T T 9 T 2 T 1 e 2 e e 3 1 T T ... T 10 13 3 T T T T 4 6 7 5 8 T T 14 Fig. 4. A configuration in a triangle space with sufficient resolution. The asso ciated incidence matrix is Eulerian (see text) but do es not satisfy Camion’s criterion, th us is not totally unimo dular. 14 Thomas Schoenemann, Simon Masnou, and Daniel Cremers T able 1. A square incidence matrix asso ciated with the configuration in Figure 4. It is Eulerian, i.e. the sum along each line and the sum along each column are even, but the total sum is not divisible by four. According to Camion [7], the matrix is not totally unimodular. T 1 T 5 T 9 T 1 , 2 T 2 , 3 T 3 , 4 T 4 , 5 T 2 , 5 T 5 , 6 T 1 , 6 T 1 , 7 T 7 , 8 T 2 , 8 T 1 , 9 T 9 , 10 T 10 , 11 T 1 , 11 T 1 , 12 T 12 , 13 T 9 , 13 T 9 , 5 T 5 , 14 T 14 , 15 T 9 , 15 T 9 , 16 T 16 , 17 T 5 , 17 P ( T 1 , e 1 ) 1 -1 -1 -1 -2 ( T 2 , e 1 ) -1 -1 -1 -1 -4 ( T 3 , e 1 ) -1 -1 -2 ( T 4 , e 1 ) -1 -1 -2 ( T 5 , e 1 ) 1 -1 -1 -1 -2 ( T 6 , e 1 ) -1 -1 -2 ( T 7 , e 1 ) -1 -1 -2 ( T 8 , e 1 ) -1 -1 -2 ( T 1 , e 2 ) 1 -1 -1 -1 -2 ( T 9 , e 2 ) 1 -1 -1 -1 -2 ( T 10 , e 2 ) -1 -1 -2 ( T 11 , e 2 ) -1 -1 -2 ( T 12 , e 2 ) -1 -1 -2 ( T 13 , e 2 ) -1 -1 -2 ( T 9 , e 3 ) 1 -1 -1 -1 -2 ( T 5 , e 3 ) 1 -1 -1 -1 -2 ( T 14 , e 3 ) -1 -1 -2 ( T 15 , e 3 ) -1 -1 -2 ( T 16 , e 3 ) -1 -1 -2 ( T 17 , e 3 ) -1 -1 -2 plus 7 lines ( T 18 , e 3 ) , · · · , ( T 24 , e 3 ) with only 0 entries to hav e a square matrix P 2 2 2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -42 A linear programming approach to the discrete Willmore b oundary problem 15 4.3 T esting the relaxed linear problem W e hav e tested the relaxed form ulation on a few examples at lo w-resolution using the dual simplex metho d implemen ted in the CLP solv er. The main reason for using low-resolution is that the num b er of triangles b ecomes significan tly imp ortan t as the resolution increases, and both the computational cost and the memory requiremen ts tend to b ecome large. Another reason for working at low-resolution is that there is no need to go high b efore finding a case of non-in tegrality . Indeed, consider the examples in figure 5: integral solutions are obtained when the resolution is very low (i.e. when there is no risk to hav e configurations like in figure 4). In the last configuration, how ever, the optimal solution of the relaxed problem has fractional entries. This confirms that our initial problem cannot b e addressed though the classical tec hniques of relaxation, and with usual LP solv ers. 4.4 On in teger linear programming Our results ab o ve indicate that, necessarily , in teger linear solvers [28,1] should b e used. These commonly start with solving the linear programming relaxations, then derive further v alid inequalities (called cuts ) and/or apply a branch-and- b ound scheme. Due to the small num ber of fractional v alues that we hav e ob- serv ed in our exp eriments, it is quite likely that the deriv ation of a few cuts only w ould give integral solutions. Ho wev er, we did not test this so far b ecause of the running times of this approach: in cases where we get fractional solutions the dual simplex metho d often needs as long as tw o w eeks and up to 12 GB memory! F rom exp erience with other linear programming problems we consider it likely that the interior p oint metho ds implemented in commercial solvers will b e m uch faster here (we exp ect less than a da y). A t the same time, we exp ect the memory consumption to b e considerably muc h higher, so the metho d would most probably b e unusable in practice. W e strongly b eliev e that a sp ecific integer linear solver should be dev elop ed rather than using general implemen tations. It is w ell kno wn that, for a few problems lik e the knapsac k problem [28][c hapter 24.6], their spe cific structure giv es rise to ad-hoc efficient approaches. Recalling that our incidence matrix is v ery sparse and well structured (the nonzero en tries of each column are either exactly tw o ( − 1), or exactly three 1) we strongly b eliev e that an efficient in teger solv er can b e developed and our approach can b e amenable to higher-resolution results in the near future. 5 Conclusion W e hav e shown that the minimization under b oundary constraints of mean cur- v ature based energies o ver surfaces, and in particular the Willmore energy , can b e cast as an in teger linear program. Unfortunately , this in teger program is not equiv alen t to its relaxation so the classical LP algorithms offer no warran ty 16 Thomas Schoenemann, Simon Masnou, and Daniel Cremers Fig. 5. A series of exp eriments (the result and the mesh edges) with increasing resolu- tion of the triangle space (and v arious b oundary constraints). An integral solution of the relaxed problem is obtained by a standard LP-solv er in b oth top cases. As for the last case, the triangle space resolution is now large enough for having configurations similar to the counterexample of figure 4. And indeed, an optimal solution is found for the relaxed problem that is not integral. The mesh on the b ottom-righ t sho ws actually t wo nested semi-spheres whose triangles hav e, at least for a few of them, non binary lab els. A linear programming approach to the discrete Willmore b oundary problem 17 that the integer optimal solution will b e found. This implies that pure integer linear algorithms must b e used, whic h are in general muc h more inv olved. W e b eliev e how ever that the particular structure of the problem pav es the wa y to a dedicated algorithm that would provide high-resolution glob al minimizers of the Willmore b oundary problem and generalizations. This is the purp ose of future researc h. References 1. T. Ach terb erg. Constr aint Inte ger Pr o gr amming . PhD thesis, T ec hnische Universit¨ at Berlin, 2007. 2. F. Almgren, J.E. T aylor, and L.-H. W ang. Curv ature-driven flows: a v ariational approac h. SIAM Journal on Contr ol and Optimization , 3:387–438, 1993. 3. A.A. Amini, T.E. W eymouth, and R.C. Jain. Using dynamic programming for solving v ariational problems in vision. IEEE T r ans. on Patt. Anal. and Mach. Intel l. , 12(9):855 – 867, September 1990. 4. A.I. Bobenko and P . Sc hr¨ oder. Discrete Willmore Flo w. In Eur o gr aphics Symp osium on Ge ometry Pro c essing , 2005. 5. Y. Boyk ov and V. Kolmogorov. An exp erimen tal comparison of min-cut/max-flo w algorithms for energy minimization in computer vision. In A.K. Jain M. Figueiredo, J. Zerubia, editor, Int. Workshop on Ener gy Minimization Metho ds in Computer Vision and Pattern R e c o gnition (EMMCVPR) , volume 2134 of LNCS , pages 359– 374. Springer V erlag, 2001. 6. K.A. Brakke. The surface evolv er. Exp erimental Mathematics , 1(2):141–165, 1992. 7. P . Camion. Characterization of totally unimo dular matrices. Pr o c. Am. Math. So c. , 16(5):1068–1073, 1965. 8. U. Clarenz, U. Diew ald, G. Dziuk, M. Rumpf, and R. Rusu. A finite element metho d for surface restoration with smo oth b oundary conditions. Computer Aide d Ge ometric Design , 21(5):427–445, 2004. 9. M. Droske and M. Rumpf. A level set form ulation for Willmore flow. Interfac es and F r e e Boundaries , 6(3), 2004. 10. G. Dziuk. Computational parametric Willmore flow. Numer. Math. , 111:55–80, 2008. 11. S. Esedoglu, S. Ruuth, and R.Y. Tsai. Threshold dynamics for high order geometric motions. T ec hnical rep ort, UCLA CAM rep ort, 2006. 12. L.R. F ord and D. F ulkerson. Flows in Networks . Princeton Univ ersity Press, Princeton, New Jersey , 1962. 13. L. Grady . Minimal surfaces extend shortest path segmentation metho ds to 3d. IEEE T r ans. on Patt. Anal. and Mach. Intel l. , 2009. T o app ear. http://cns- web.bu.e du/ lgr ady/ . 14. R. Grzibovskis and A. Hein tz. A con volution-thresholding scheme for the Willmore flo w. T echnical Rep ort Preprin t 34, Chalmers Universit y of T echnology , G¨ oteb org, Sw eden, 2003. 15. K. Hildebrandt, K. P olthier, and M. W ardetzky . On the conv ergence of metric and geometric prop erties of p olyhedral surfaces. Ge ometriae De dic ata , 123:89–112, 2005. 16. L. Hsu, R. Kusner, and J. Sulliv an. Minimizing the squared mean curv ature integral for surfaces in space forms. Exp erimental Mathematics , 1(3):191–207, 1992. 18 Thomas Schoenemann, Simon Masnou, and Daniel Cremers 17. R. Kusner. Comparison surfaces for the Willmore problem. Pacific J. Math. , 138(2), 1989. 18. S. Luc khaus and T. Sturzenhec k er. Implicit time discretization for the mean curv ature flow equation. Calculus of variations and p artial differ ential e quations , 3(2):253–271, 1995. 19. U.F. May er and G. Simonett. A n umerical scheme for axisymmetric solutions of curv ature-driv en free b oundary problems, with applications to the willmore flow. Interfac es and F r e e Boundaries , 4:89–109, 2002. 20. J.-M. Morv an. Gener alize d Curvatur es . Springer Publishing Company , Incorp o- rated, 1 edition, 2008. 21. N. Olischl¨ ager and M. Rumpf. Tw o step time discretization of willmore flow. In Pr o c. 13th IMA Int. Conf. on Math. of Surfac es , pages 278–292, Berlin, Heidelb erg, 2009. Springer-V erlag. 22. U. Pink all. Hopf tori in S 3 . Invent. Math. , 81:379–386, 1985. 23. U. Pink all and I. Sterling. Willmore surfaces. The Mathematic al Intel ligenc er , 9(2), 1987. 24. K. Polthier. Polyhe dr al surfac es of c onstant me an curvature . Habilitation thesis, TU Berlin, 2002. 25. K. Polthier. Computational asp ects of discrete minimal surfaces. In Glob al The ory of Minimal Surfac es , Proc. of the Cla y Mathematics Institute Summer Sc ho ol, 2005. 26. R. Rusu. An algorithm for the elastic flo w of surfaces. Interfac es and F r e e Bound- aries , 7:229–239, 2005. 27. R. Sc h¨ atzle. The Willmore boundary problem. Calc. V ar. Part. Diff. Equ. , 37:275– 302, 2010. 28. A. Schrijv er. The ory of line ar and inte ger pr o gr amming . Wiley-Interscience series in discrete mathematics. John Wiley and Sons, July 1994. 29. L. Simon. L e ctur es on Ge ometric Me asur e The ory , volume 3 of Pr o c. of the Center for Mathematic al Analysis . Australian National Univ ersity , 1983. 30. J.M. Sulliv an. A Crystal line Appr oximation The or em for Hyp ersurfac es . PhD thesis, Princeton Universit y , Princeton, New Jersey , 1992. 31. J.M. Sulliv an. Computing hypersurfaces whic h minimize surface energy plus bulk energy . Motion by Mean Curvatur e and R elate d T opics , pages 186–197, 1994. 32. J. V erdera, V. Caselles, M. Bertalmio, and G. Sapiro. Inpainting surface holes. In In Int. Confer enc e on Image Pr o c essing , pages 903–906, 2003. 33. M. W ardetzky , M. Bergou, D. Harmon, D. Zorin, and E. Grinspun. Discrete quadratic curv ature energies. Computer A ide d Ge ometric Design , 24(8–9):499–518, 2007.