Antonin Chambolle
email: antonin.chambolle@polytechnique.fr

Daniel Cremers
email: cremers@tum.de

Thomas Pock
email: pock@icg.tugraz.at

A convex approach to minimal partitions

We describe a convex relaxation for a family of problems of minimal perimeter partitions. The minimization of the relaxed problem can be tackled numerically, we describe an algorithm and show some results. In most cases, our relaxed problem finds a correct numerical approximation of the optimal solution: we give some arguments to explain why it should be so, and also discuss some situation where it fails.

Introduction 1.Contribution

We present an approach which allows to numerically compute solutions to the minimal partition problem and to a few related problems. The considered optimization problems arise in many fields of science. Motivated by image analysis applications like image segmentation, they have been the subject of extensive study in the beginning of the 90s [START_REF] Morel | Segmentation of images by variational methods: a constructive approach[END_REF][START_REF] Bellettini | An almost everywhere regularity result for minimal partitions[END_REF][START_REF] Massari | Regularity properties of optimal segmentations[END_REF][START_REF] Tamanini | Optimal segmentation of unbounded functions[END_REF], in connection to the celebrated Mumford-Shah [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] segmentation problem.

After a first (unpublished) version of this paper [START_REF] Chambolle | A convex approach for computing minimal partitions[END_REF] had been written and a conference version [START_REF] Pock | A convex relaxation approach for computing minimal partitions[END_REF] had been published, several approaches were proposed in the image processing literature for tackling this problem. See, for instance, [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF][START_REF] Lellmann | Convex multi-class image labeling by simplex-constrained total variation[END_REF][START_REF] Brown | A convex relaxation method for a class of vector-valued minimization problems with applications to Mumford-Shah segmentation[END_REF][START_REF] Bae | Global minimization for continuous multiphase partitioning problems using a dual approach[END_REF]. However, it seems that most of them solve a problem which is quite "far" (at least, further than the approaches we will discuss here) from the original problem so that we still believe that our contribution is useful and, to some extent, helps clarifying what seems possible and what is not.

The problem of (numerically) finding a partition of a set Ω ⊂ R d , which minimizes the (d -1)-dimensional measure of the total interface, plus either boundary conditions, or some external field particular to each set, is a challenging task. Its discrete version, known as the "Potts' model" [START_REF] Potts | Some generalized order-disorder transformations[END_REF] (an extension of the ferromagnetic Ising model [START_REF] Ising | Beitrag zur Theorie des Ferromagnetismus[END_REF]), is described by an energy whose minimization is NP-hard. What we propose here is to derive some convexification of the problem, which is close enough to the convex envelope (in fact, it is the closest in a particular, reasonably large class), but, also, seems to be numerically tractable, at least when the number of labels is not too high. Then, minimizing numerically this envelope, we experience that in general, at least when d = 2, the minimizer is actually a minimal partition. Figure 1 shows an example where the proposed method is applied to solve a piecewise constant Mumford-Shah model with 10 regions. We point out that related experiments first appeared in a paper by K. Brakke of 1995 [START_REF] Brakke | Soap films and covering spaces[END_REF], and that it was known since that time that in some cases, the problem was indeed convex, see also [START_REF] Lawlor | Paired calibrations applied to soap films, immiscible fluids, and surfaces or networks minimizing other norms[END_REF]. We will describe both representations of [START_REF] Brakke | Soap films and covering spaces[END_REF] and [START_REF] Chambolle | A convex approach for computing minimal partitions[END_REF]. Although both are equivalent up to a change of variable, strangely enough, results which are not obvious in one representation are sometimes much simpler in the other, and vice versa.

It came as a surprise to us that the method works in many cases. Nevertheless, we also provide as a numerical counterexample a situation where it does not work. The goal of this paper is to describe our setting, and explain how we perform the minimization. Then, we provide arguments to explain why it works in many cases, and why it shouldn't in other cases.

In some sense, our approach is related to similar relaxations in discrete optimization (and in particular optimization of Markov Random Fields or MRFs), such as LP-relaxation or roof duality relaxation. The point of view closest to ours in the discrete literature seems to be roof duality [START_REF] Hammer | Roof duality, complementation and persistency in quadratic 0-1 optimization[END_REF], indeed, we also look for a convex relaxation of our problem which is local in some sense, and is obtained as the supremum of affine functions which satisfy local constraints.

Let us mention however that we also use representations of the problem in which our relaxation is, in a precise sense, optimal. Then, we base our analysis on the concept of "calibrations", which may be seen as the continuous counterpart of roof duality, and was developed independently for the study of minimal surfaces in the 70's and 80's [START_REF] Brakke | Soap films and covering spaces[END_REF][START_REF] Lawlor | Paired calibrations applied to soap films, immiscible fluids, and surfaces or networks minimizing other norms[END_REF][START_REF] Alberti | The calibration method for the Mumford-Shah functional and free-discontinuity problems[END_REF].

Our extension to the continuous setting and the representation follow the spirit of [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF] where two of the authors of this paper were extending recent approaches in MRFs optimization to the continuous setting. This is crucial, in particular, if one wants to reproduce precisely specific surface energies, and in particular isotropic interfacial energies.

Another point of view similar to ours is found in a recent paper of C. Zach et al. [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF]. It is also written in a continuous setting although, there again, the interaction potentials are eventually chosen anisotropic. See also [START_REF] Lellmann | Convex multi-class image labeling by simplex-constrained total variation[END_REF] for a variant. In fact, we show that with the anisotropy chosen in [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF], both approaches boil down to the same representation, while in general, we claim that ours is better (meaning, a tighter relaxation, closer to the original problem).

Eventually, we must point out that what we propose is quite different from approaches based on phase-field approximations or level-sets (see for instance [START_REF] Chan | Active contour and segmentation models using geometric PDE's for medical imaging[END_REF]). These do not aim at computing a global minimizer of the problem. Furthermore, the latter mentioned cannot approximate very general interfacial energies, or even uniform ones.

Two representations for multi-label problems

The generic problem we will address here is of the form min

{E i } k i=1 1 2 k i=1 Per(E i ; Ω) + k i=1 E i g i (x) dx (1)
where andk i=1 E i = Ω (up to Lebesgue-negligible sets). Here, g i are "potentials" (which in practice will take various forms). Since the optimization problem is clearly invariant to replacing all functions g i by g i + g with an arbitrary integrable function g, we will assume without loss of generality that all g i are nonnegative, more specifically that g i ∈ L 1 (Ω; R + ).

{E i } k i=1 is a partition of an open set Ω ⊂ R d , d ≥ 2, into k sets: E i ∩ E j = 0 of i = j,
The first term in [START_REF] Ahuja | Network flows[END_REF] is half the sum of the perimeters of the sets E i , i = 1, . . . , k, which will be defined precisely in Section 2.2. It is the same as the surface of the total interface i<j ∂E i ∩ ∂E j (as the common surface between E i and E j is counted twice in the sum of the perimeters) and defines what is called the total perimeter of the partition {E i } k i=1 . The most classical way to represent -and relax -this problem is to introduce the characteristic functions v i = χ E i associated with each set E i . These satisfy v i (x) ≥ 0, i = 1, . . . , k, for a.e. x ∈ Ω (2)

k i=1 v i (x) = 1, a.e. x ∈ Ω (3)
v i (x) ∈ {0, 1}, i = 1, . . . , k, a.e. x ∈ Ω (4)

Then, since by definition (see Section 2.2), the perimeter of E i in Ω is the total variation Ω |Dv i | of its characteristic function, the energy in (1) can be rewritten

E(v) = 1 2 k i=1 Ω |Dv i | + k i=1 Ω v i (x)g i (x) dx (5)
Then, the most straightfoward relaxation of the problem, considered for instance by Zach et al. in [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF] (but this is also the standard approach in discrete optimization, see the pages on multiway cut problems in [START_REF] Ahuja | Network flows[END_REF]), consists in minimizing E on the convex set of BV functions with values in the simplex S = {z ∈ R k + : k i=1 z i = 1}:

S = v = (v 1 , . . . , v k ) ∈ BV (Ω; [0, 1] k ) : v i ≥ 0 , k i=1 v i = 1 a.e. in Ω . (6)
That is, we have kept the constraints (2-3), and relaxed (4) into v i ∈ [0, 1] for each i. This is essentially what is also done in [START_REF] Brown | A convex relaxation method for a class of vector-valued minimization problems with applications to Mumford-Shah segmentation[END_REF][START_REF] Bae | Global minimization for continuous multiphase partitioning problems using a dual approach[END_REF]. It turns out that this approach can give wrong results. In Section 5, we will show both theoretically and experimentally that this happens more often than for the method presented in this paper -see Fig. 13 and Proposition 5.1.

Another, different representation of the problem is proposed in [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF] and was adopted in the previous version of this paper [START_REF] Chambolle | A convex approach for computing minimal partitions[END_REF]. It consists in using, instead of v = (v 1 , . . . , v k ), the variables u = (u 1 , . . . , u k-1 ) defined by

u i = 1 - i j=1 v j , i = 1, . . . , k -1 (7)
which satisfy the monotonicity constraint

1 ≥ u 1 (x) ≥ • • • ≥ u k-1 (x) ≥ 0 , for a.e. x ∈ Ω . ( 8
)
Letting by convention (which is coherent with [START_REF] Bae | Global minimization for continuous multiphase partitioning problems using a dual approach[END_REF] and will be implicit in the rest of the paper) u 0 ≡ 1 and u k ≡ 0, we recover v from u by the change of variable

v i = u i-1 -u i , i = 1, . . . , k. (9)
Linear (TV) The function u i is also the characteristic function of the upper level set {ι(x) > i} of the label ι(x) = {i : v i (x) = 1} ∈ {1, . . . , k}, defined for a.e. x ∈ Ω.

The u representation makes more sense than the previous one whenever the labels form an ordered set (for instance, it they correspond to disparity measurements in stereo reconstruction [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF]). In particular, we will also consider interfacial energies σ i,j (which penalize the surface of the interfaces ∂E i ∩ ∂E j ) which are different from 1 and can depend on (i, j), and this representation is particularly well suited for penalizations of the form σ(ij), such as the total variation |i -j| (linear interaction potential, see Fig. 2, left), or the "truncated total variation" min{|i -j|, T } (truncated linear potential, Fig. 2, right).

On the other hand, the v representation is more justified when the labels have an arbitrary meaning (such as potatoes, tomatoes, other...) and the interaction energy is unrelated to their values, such as in the minimal partition problem (1) (or the "Pott's" model), see Fig. 2, middle.

In the discrete setting, the u representation has been first introduced in the seminal works of Ishikawa and Geiger [START_REF] Ishikawa | Segmentation by grouping junctions[END_REF][START_REF] Ishikawa | Exact optimization for Markov random fields with convex prio rs[END_REF]. It was then adapted to continuous problems of image reconstruction in [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF].

Outline

The paper is organized as follows: a preliminary section gives the definition and some basic facts about functions with bounded variation and finite-perimeter partitions. In Section 3 we formalize the minimal partition problem and introduce first its convex, lower-semicontinuous envelope by means of Legendre-Fenchel conjugation. This abstract approach, however, does not in general lead to a problem which is computationally tractable: for this reason, we introduce a notion of "local" convex envelope, which is less tight but easier to implement. We show (Prop. 3.3) that the tightest "local" convex envelope reduces to the duality approach through "paired calibrations," which were introduced in [START_REF] Lawlor | Paired calibrations applied to soap films, immiscible fluids, and surfaces or networks minimizing other norms[END_REF][START_REF] Brakke | Soap films and covering spaces[END_REF] in order to provide optimality conditions for the partition problem (or more general minimal surfaces problems). From the theoretical point of view, one of the main contributions of this paper are the bounds [START_REF] Ishikawa | Segmentation by grouping junctions[END_REF], which, in regards of the bounds [START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF] for the global convex envelope, show that the local convex envelope is close to optimality. This can be related to similar estimates for the discrete Potts model (see for instance [START_REF] Ahuja | Network flows[END_REF]).

In Section 4 we present a simple numerical approximation of the minimal partition problem. We use a simple yet provably convergent primal-dual algorithm to solve the arising convex problem (which we first described in [START_REF] Pock | A convex relaxation approach for computing minimal partitions[END_REF]). In Section 4.3, we present a number of promising experimental results on a variety of segmentation and geometric inpainting problems which show the potential and limitations of our approach.

In Section 5, we compare our approach to the convexifications proposed in the recent papers [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF] and [START_REF] Lellmann | Convex multi-class image labeling by simplex-constrained total variation[END_REF]. Another important contribution of this paper is the proof, not only experimental, that the envelope studied in this paper is, in general, strictly higher (Prop. 5.1). Nevertheless, the latter are simpler (and [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF] is, in fact, equivalent to our approach in 1D as well as some anisotropic cases). In Section 6, we consider generalizations of the partition problem by assigning weights to respective label transitions. This allows to tackle, for instance, the case of truncated linear potentials. We show experimental comparisons of different models. An important result in that section is Prop. 6.1, where we show that in the representation [START_REF] Bae | Global minimization for continuous multiphase partitioning problems using a dual approach[END_REF], one can get rid of some constraints and solve a simpler problem, and yet recover the correct solution. This is what was experimentally observed in [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF].

A few technical appendices follow, where we explain how to perform practically the pro-jections on some of the convex sets introduced throughout the paper, and how to compute the value of our convexification in the 3-labels case -while there is no closed form in general. Of particular interest is Appendix A.3, where we show that in many cases the problem, locally, behaves like a problem with fewer phases (in practice, the saturated phases, which are 2 or 3 in dimension 2).

Preliminaries

Functions with bounded variation

We recall here some results on functions with bounded variations, see [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF][START_REF] Giusti | Minimal surfaces and functions of bounded variation[END_REF][START_REF] Evans | Measure theory and fine properties of functions[END_REF][START_REF] Ziemer | Weakly differentiable functions[END_REF] for more details. Let d, k ≥ 1, and Ω be an open subset of R d . Given a function v = (v 1 , . . . , v k ) ∈ L 1 (Ω; R k ), we define its total variation as sup -

Ω k i=1 v i div ξ i : ξ ∈ C 1 c (Ω; R d ) k , |ξ(x)| 2 = k i=1 |ξ i | 2 ≤ 1 ∀x ∈ Ω . ( 10
)
We say that v has bounded variation whenever the value of ( 10) is finite. In this case, one shows (using Riesz' representation theorem) that the distributional derivative Dv is a bounded (R k×d -valued) Radon measure, whose total mass Ω |Dv| is precisely the value of [START_REF] Birchfield | A pixel dissimilarity measure that is insensitive to image sampling[END_REF]. We denote by BV (Ω; R k ) the space of (k-dimensional vector valued) functions with bounded variation. Endowed with the norm v BV = v L 1 (Ω;R k ) + Ω |Dv|, it is a Banach space. We use the notation BV (Ω) whenever k = 1, and BV (Ω; K) for functions v such that v(x) ∈ K a.e., where K is a closed, convex subset of R k .

The distributional derivative Dv has the decomposition [START_REF] Federer | Geometric measure theory[END_REF][START_REF]′ pert. Spaces BV and quasilinear equations[END_REF]:

Dv = ∇v dx + Cv + (v + -v -) ⊗ ν v H d-1 J v . (11)
Here, ∇v dx is the part of Dv which is absolutely continuous with respect to the Lebesgue measure. The approximate gradient ∇v(x) is a L 1 vector field, which corresponds to the weak gradient whenever v ∈ W 1,1 . The jump set J v is the sets of points x ∈ Ω where there exist

v -(x) = v + (x) ∈ R k and ν v (x) ∈ S d-1 such that (y → v(x + εy)) L 1 (B(0,1);R k ) -→ y → χ {y•νv(x)>0} v + (x) + χ {y•νv(x)>0} v -(x)
as ε → 0. Of course, here, we could replace the triplet (v

+ , v -, ν v ) with (v -, v + , -ν v ): when k = 1, the convention is to choose u + > u -, if k > 1,
we can for instance choose that the component with the lowest index i such that (v

+ ) i = (v -) i satisfies (v + ) i > (v -) i . Then, the tensor product (v + -v -) ⊗ ν v is the k × d matrix ((v + -v -) i (ν v ) j ) i=1,...,
Ψ(x, p) = lim t→0 tΨ(x, p/t) (12)
which is a convex, one-homogeneous function (possibly taking the value +∞). Then, [START_REF] Bouchitté | Integral representation of convex functionals on a space of measures[END_REF][START_REF] Rockafellar | Integrals which are convex functionals[END_REF]]

Ω Ψ(x, Dv) = Ω Ψ(x, ∇v(x)) dx + Ω Ψ x, Cv |Cv| d|Cu| + Jv Ψ(x, (v + (x) -v -(x)) ⊗ ν v (x)) dH d-1 (x). ( 13
)
Moreover, if Ψ * (x, •) is the Legendre-Fenchel conjugate [START_REF] Ekeland | Convex analysis and variational problems[END_REF][START_REF] Rockafellar | Convex analysis[END_REF] of Ψ(x, •) (w.r. the "p" variable), one also have the dual representation

Ω Ψ(x, Dv) = sup Ω v(x) • div ξ(x) -Ψ * (x, ξ(x)) dx : ξ ∈ C ∞ c (Ω; R k×d ) (14)
(here div ξ is the vector ( d j=1 (∂(ξ i ) j /∂x j )) k i=1 ). We have the following approximation result, whose proof follows Meyers-Serrin's classical proof for Sobolev functions and is found for instance in [START_REF] Giusti | Minimal surfaces and functions of bounded variation[END_REF]:

Theorem 2.1. For any v ∈ BV (Ω; R k ), there exists a sequence (v n ) n≥1 of functions in C ∞ (Ω; R k ) such that v n → v in L 1 (Ω; R k ) and lim n→∞ Ω |∇v n (x)| dx = |Dv|(Ω). (15)
Combining this result with a celebrated theorem of Reshetnyak [3, Theorem 2.39], we get in addition that, if Ψ is a convex, continuous and one-homogeneous function over R k×d ,

lim n→∞ Ω Ψ(∇v n (x)) dx = Ω Ψ(Dv) . ( 16
)

Caccioppoli sets

If E ⊂ Ω is a measurable set, then its perimeter in Ω is defined as the total variation [START_REF] Birchfield | A pixel dissimilarity measure that is insensitive to image sampling[END_REF] of χ E . A Caccioppoli set, or set with finite perimeter, is a set such that χ E ∈ BV (Ω): then we let Per(E, Ω) := Ω |Dχ E |. In this case, the (rectifiable) jump set J χ E is also called the "reduced boundary" of E, denoted by ∂ * E, and is equal, up to a set of H d-1 measure zero, to Ω \ (E 1 ∪ E 0 ), where E 1 , resp., E 0 , are the points where E has Lebesgue density 1, resp., 0. The normal vector ν χ E (as defined in the previous section), which we also denote by ν E , is the inner normal vector to ∂ * E, defined H d-1 -a.e. on the boundary. A Caccioppoli partition of Ω is a (finite or countable) sequence of subsets of Ω,

(E i ) i≥1 such that |E i ∩ E j | = 0 if i = j, |Ω \ i E i | = 0 (equivalently, i χ E i = 1 a.e.

in Ω), and

i Per(E i , Ω) < +∞. The total perimeter of the partition is half the sum of the perimeters of the sets, since in the latter sum each interface between two sets E i and E j is counted twice.

We recall eventually that if u ∈ BV (Ω) is a scalar-valued BV function, it enjoys the co-area formula:

Ω |Du| = +∞ -∞ Per({u > s}; Ω) ds . (17)
That is, the total variation of a function is recovered as the total sum of the surfaces of its level sets.

In the next section, we introduce a class of variational problem whose unknown is a Caccioppoli partition, with a fixed, maximal number of sets k.

Minimal partitions

The convex envelope of the minimal partition problem

Let us first, to simplify, focus on the problem of computing a minimal partition of Ω, bounded open subset of R d . We will generalize the problem afterwards. We assume we want to find a partition of a set Ω into k (at most) sets E 1 , . . . , E k , which solves problem (1), for given functions g 1 , . . . , g k ∈ L 1 (Ω; R + ).

As already mentioned, the first term in energy (1) is the total length of the boundaries of the partition: the weight 1/2 is there to take into account that each interface between to sets E i and E j is contained twice in the sum, as the boundary of each set.

The existence of a solution to (1) is a straightforward result of calculus of variations (which follows from Rellich's theorem and the lower-semicontinuity of the total variation, see for instance [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF]), and we will only focus on the problem of actually finding it.

If we introduce the functions v i = χ E i , then a partition is a vector v = (v 1 , . . . , v k ) ∈ BV (Ω; {0, 1} k ) which satisfies conditions [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF]. Then, the interfacial energy of the partition is nothing else as H d-1 (J v ), and, as mentioned in the introduction, it coincides with F(v) where the convex, one-homogeneous function F : BV (Ω; R k ) → [0, +∞) is defined as

F(v) = 1 2 k i=1 Ω |Dv i | and extended to L 2 (Ω; R k ) by letting F(v) = +∞ if v ∈ BV (Ω; R k ).
We therefore define J : L 2 (Ω; R k ) → [0, +∞] by

J (v) =    F(v) if v ∈ BV (Ω; {0, 1}) , k i=1 v i = 1 a.e. , +∞ else.
Then, problem (1) can be written equivalently min

v∈L 2 (Ω;R k ) J (v) + Ω v(x) • g(x) dx (18)
where g = (g 1 , . . . , g k ) ∈ L 1 (Ω; R k + ) (in fact, we should require g ∈ L 2 in order for the last integral to make sense, but observe that even for g ∈ L 1 is this integral a continuous functional of v on the domain of J , which is the set where it is finite. We will not bother about these details, since in practice these functions are bounded...) Problem ( 18) is nonconvex since, although F is convex, the domain of J , defined as the set dom J = {v ∈ L 2 (Ω; R k ) : J (v) < +∞}, which clearly is the set

S 0 = v = (v 1 , . . . , v k ) ∈ BV (Ω; {0, 1} k ) : v i ≥ 0 , k i=1 v i = 1 a.e. in Ω , (19)
is not convex. However, as we already observed, it is standard that (18) admits a solution.

How could we compute this solution? A first, natural idea would be to identify the convex envelope of J (and hence of [START_REF] Chambolle | Theoretical Foundations and Numerical Methods for Sparse Recov ery[END_REF], since the other term in the problem is linear). Indeed, if we let, for w ∈ L 2 (Ω; R k )

J * (w) = sup v∈L 2 (Ω;R k ) Ω v(x) • w(x) dx -J (v)
be the Legendre-Fenchel conjugate of J (see [START_REF] Ekeland | Convex analysis and variational problems[END_REF][START_REF] Rockafellar | Convex analysis[END_REF]) and then, again, for v ∈ L 2 (Ω; R k )

J * * (v) = sup w∈L 2 (Ω;R k ) Ω v(x) • w(x) dx -J * (w)
be the Legendre-Fenchel conjugate of J * , it is well known [START_REF] Ekeland | Convex analysis and variational problems[END_REF][START_REF] Rockafellar | Convex analysis[END_REF] that J * * is the convex, lower-semicontinuous envelope of J , and, as well, that the minimizers of ( 18) are also minimizers of min

v∈L 2 (Ω;R k ) J * * (v) + Ω v(x) • g(x) dx . (20)
Conversely, minimizers of (20) can only be convex combinations of minimizers of the original problem [START_REF] Chambolle | Theoretical Foundations and Numerical Methods for Sparse Recov ery[END_REF] (in case of non uniqueness, of course, it might be impossible to identify a minimizer of the latter from a minimizer of the former).

However, we do not know any explicit form of J * * , and believe there is none in general. It is standard that another way to compute J * * is through the formula

J * * (v) = inf lim inf n→∞ Nn α=1 θ n α J (v n α ) : N n ≥ 1 , θ n α ≥ 0 , Nn α=1 θ n α = 1 , Nn α=1 θ n α v n α → v as n → ∞
which proves useful in some special cases (for instance, if v i = c i = constant for each i, with c i ≥ 0 and i c i = 1, we easily deduce J * * (v) = 0), but does not help in general.

Still, we can compute the domain of J * * and provide some useful estimates:

Proposition 3.1. The domain of J * * is dom J * * = S, given by (6). We have J * * (v) = J (v) whenever v ∈ S 0 = dom J , and for each v ∈ S,

F(v) ≤ J * * (v) ≤ (k -1)F(v) . (21)
If, moreover, (k -2)/(k(k -1)) ≤ v i ≤ 2/k a.e.
in Ω, for each i = 1, . . . , k, then the estimate is improved:

F(v) ≤ J * * (v) ≤ 2 k -1 k F(v) . (22)
The last condition expresses the fact that v lies close to the "mixture" (1/k, . . . , 1/k) where all k phases are uniformly spread in the domain: more precisely, v = θw+(1-θ)(1/k, . . . , 1/k) for some w ∈ S and θ ∈ [0, 1/(k -1)].

Observe that for k = 2, the estimates are optimal and just show that J * * = F: it is well-known in that case that, thanks to the co-area formula [START_REF] Brown | A convex relaxation method for a class of vector-valued minimization problems with applications to Mumford-Shah segmentation[END_REF], the problem of finding a partition into two sets is essentially convex, and is relaxed by minimizing (5) over S. Then, a partition is found by letting E 1 = χ {v 1 >s} , for any s ∈ [0, 1), and E 2 = Ω \ E 1 , see for instance [START_REF] Bellettini | Convex approximations of functionals with curvature[END_REF][START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF][START_REF] Chambolle | Theoretical Foundations and Numerical Methods for Sparse Recov ery[END_REF].

Proof. First, since F is itself convex, lower-semicontinuous, and F ≤ J (J = F on S 0 = dom J and +∞ elsewhere), the inequality

F ≤ J * * is straightforward. In particular, if v ∈ S 0 , F(v) ≤ J * * (v) ≤ J (v) = F(v), and we deduce J * * (v) = J (v).
If v ∈ S, it must be that J * * (v) = +∞: indeed, it is classical that dom J * * ⊂ convdom J , which shows that in particular, if J * * (v) < +∞, then v i ≥ 0 and k i=1 v i = 1 a.e. in Ω. On the other hand, v must have bounded variation, otherwise J (v) * * ≥ F(v) = +∞, a contradiction. We deduce that dom J * * ⊂ S.

Let now v ∈ S and u ∈ BV (Ω; [0, 1] k-1 ) be defined by ( 9) (we also let u 0 = 1, u k = 0). For s ∈ [0, 1], we let u s = (u s 1 , . . . , u s k-1 ) where u s i = χ {u i >s} for i = 0, . . . , k: it is for each x a binary vector with 1 = u s 0

(x) ≥ u s 1 (x) ≥ • • • ≥ u s k-1 (x) ≥ u s k (x) = 0.
The co-area formula [START_REF] Brown | A convex relaxation method for a class of vector-valued minimization problems with applications to Mumford-Shah segmentation[END_REF] yields that for each i and a.e. s ∈ [0, 1], u s i ∈ BV (Ω; {0, 1}), moreover, we have u i (x) = 1 0 u s i (x) ds for a.e. x and each i.

Letting therefore, following [START_REF] Bae | Global minimization for continuous multiphase partitioning problems using a dual approach[END_REF], v s i = u s i-1u s i , we find that v s ∈ S 0 for a.e. s ∈ [0, 1], and

v i = 1 0 v s i ds.
By convexity, we deduce that (it requires, in fact, an approximation argument which is standard and which we skip)

J * * (v) ≤ 1 0 J (v s ) ds. (23)
Now, for a.e. s,

J (v s ) = H d-1 (J v s ) = H d-1 (J u s ) is the total surface of the jump set of u s . Since J u s ⊂ k-1 i=1 J u s i , it follows that J (v s ) ≤ k-1 i=1 H d-1 (J u s i ) = k-1 i=1 Ω |Du s i | . (24)
From ( 23), ( 24), ( 17) and ( 7), we find

J * * (v) ≤ k-1 i=1 Ω |Du i | ≤ k-1 i=1 i j=1 Ω |Dv j | = k-1 j=1 (k -j) Ω |Dv j |. (25)
The last estimate (25) cannot, in fact, depend on the order of the labels, and we deduce that for any permutation σ of {1, . . . , k},

J * * (v) ≤ k-1 j=1 (k -j) Ω |Dv σ(j) |.
In particular, we may assume (after a suitable rearrangement) without loss of generality that

Ω |Dv 1 | ≤ Ω |Dv 2 | ≤ • • • ≤ Ω |Dv k |, in which case one can show that k-1 j=1 (k -j) Ω |Dv j | ≤ k -1 2 k j=1 Ω |Dv j | ,
which together with (25) yields the right-hand side of [START_REF] Chan | Active contour and segmentation models using geometric PDE's for medical imaging[END_REF]. We conclude also that dom J * * = S. This proof is strongly inspired from similar estimates for the discrete multiway cut problem, see [START_REF] Ahuja | Network flows[END_REF]. Eventually, we prove [START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF]. Consider first w of the following form: for some i 1 , i 2 , one has

w i 1 = w, w i 2 = 1 -w, w i = 0 if i ∈ {1, . . . , k} \ {i 1 , i 2 }, where w ∈ BV (Ω; [0, 1]). Then, w ∈ S, moreover, w = 1 0 w s ds where for s ∈ (0, 1), w s i 1 = χ { w>s} , w s i 2 = χ { w≤s} , w s i = 0 for i = i 1 , i 2 .
We deduce from the convexity of J * * and the co-area formula (17) that

J * * (w) ≤ 1 0 J (w s ) ds = Ω |D w| = F(w) ,
which together with [START_REF] Chan | Active contour and segmentation models using geometric PDE's for medical imaging[END_REF] shows that for two-phases states such as w, one has, in fact,

J * * (w) = F(w) . ( 26
) Let now w ∈ S and v = (1/(k -1))w + (1 -1/(k -1))(1/k, . . . , 1/k) (which spans all v ∈ S with (k -2)/(k(k -1)) ≤ v i ≤ 2/k a.e.
). We will write v as a convex combination of states with only two phases for which the energy is estimated with [START_REF] Evans | Measure theory and fine properties of functions[END_REF]. Let ŵi , i = 1, . . . , k, be as follows:

ŵi j =    w i if j = i , 1 -w i k -1 else.
Then, ŵi is the average of the (k -1) states for which the ith component is w i and one of the other components is 1w i , so that by [START_REF] Evans | Measure theory and fine properties of functions[END_REF],

J * * (w i ) ≤ Ω |Dw i | = (k -1) Ω |Dv i | .
where we have used

v = (1/(k -1))w + constant. Now, direct calculation shows that v = (1/k) k i=1
ŵi , so that

J * * (v) ≤ k -1 k k i=1 Ω |Dv i | = 2 k -1 k F(v) ,
which shows ( 22), as expected.

Remark 3.2. Observe that in the proof of ( 21), it is useful to consider not only the "v" representation of the phases, but also the "u" representation obtained through the change of variable ( 9).

The "local" convex envelope

Since it does not seems tractable to design an algorithm minimizing the unknown functional J * * , we will now look for a "convex envelope" which shares, in addition, the property of being local, in the sense that it can be written roughly in the form [START_REF] Bouchitté | Integral representation of convex functionals on a space of measures[END_REF], where we will assume in addition that Ψ(x, •) is even. In fact, a convex local functional may have a more general form, as detailed in [START_REF] Bouchitté | Integral representation and relaxation of convex local functionals on BV(Ω)[END_REF], however, we do not believe that considering such form would improve the readability of this paper, nor give essentially different results, especially in view of Proposition 3.3. We consider therefore a nonnegative function Ψ : Ω × R k×d → [0, +∞), continuous, and convex in the last variable, such that

Ω Ψ(x, Dv) ≤ J (v) , v ∈ L 2 (Ω; R k ). ( 27
)
We will moreover assume (as this is true, for instance, if Ψ(x, p) = 1 2 k i=1 |p i |, and we hope to do better than this integrand), that

Ω Ψ(x, Dv) = J (v) , v ∈ dom J = S 0 . (28)
A first observation is that we must have Ψ(x, 0) = 0, and this implies that for any p, t → Ψ(x, tp)/t is nondecreasing, so that Ψ ≤ Ψ (the recession function defined in [START_REF] Bouchitté | Integral representation and relaxation of convex local functionals on BV(Ω)[END_REF]).

It follows that for any v,

Ω Ψ(x, Dv) ≤ Ω Ψ(x, Dv)
and since [START_REF] Federer | Geometric measure theory[END_REF] only involves Ψ (this is because, in fact, J (v) is finite only on binary vectors, for which Dv is purely singular, so that

Ω Ψ(x, Dv) = Jv Ψ(x, ν v (x)) dH d-1
), clearly the largest possible choice is to take Ψ = Ψ, that is, to restrict the choice to convex, 1-homogeneous integrands. Therefore, from now on, we also assume that Ψ(x, •) is 1-homogeneous. We will show that, in this class, we can exhibit a maximal integrand Ψ such that ( 27) holds for all v.

Let now x 0 ∈ Ω and ρ > 0 such that B(x 0 , ρ) ⊂ Ω, choose ν ∈ S d-1 a direction and let

Γ + = ∂B(x 0 , ρ) ∩ {x : (x -x 0 ) • ν > 0}, Γ -= ∂B(x 0 , ρ) ∩ {x : (x -x 0 ) • ν < 0}, ∆ = B(x 0 , ρ) ∩ {x : (x -x 0 ) • ν = 0} Denote (f i ) k i=1 the canonical basis of R k . For i = j, we denote λ ± i,j = Γ ± Ψ(x, (f i -f j ) ⊗ ➉ ♥❻(x)) dH d-1 (x) , λ i,j = ∆ Ψ(x, (f i -f j ) ⊗ ν) dH d-1 (x)
where ➉ ♥❻(x) = (xx 0 )/|xx 0 | denotes the outer normal to B(x 0 , ρ).

Then, considering all possible v piecewise constant with J v ⊆ Γ + ∪ Γ -∪ ∆, and using ( 27) and ( 28), we find the following relationships:

λ + i,j + λ - i,j = H d-1 (Γ + ) + H d-1 (Γ -) λ + i,j + λ j,i = H d-1 (Γ + ) + H d-1 (∆) λ - i,j + λ i,j = H d-1 (Γ -) + H d-1 (∆)
for all i, j ∈ {1, . . . , k}. Combining these three equations, we deduce:

λ i,j + λ j,i = 2H d-1 (∆) = 4ρ
and dividing by 2ρ and letting then ρ → 0, it implies

Ψ(x 0 , (f i -f j ) ⊗ ν)) + Ψ(x 0 , (f j -f i ) ⊗ ν)) = 2.
Since we have assumed that Ψ(x, •) was even, we find that Ψ(x 0 , (f if j ) ⊗ ν) = 1 for all i, j. It follows, using also the homogeneity, that Ψ(x, p) = |p| for any p of the form

(f i -f j ) ⊗ p, p ∈ R d .
Hence the largest possible Ψ (in the considered class) is independent on x, and given by the largest convex function below

Ψ 0 (p) =    |p| if p = (f i -f j ) ⊗ p, 1 ≤ i, j ≤ k, p ∈ R d +∞ else.
It remains to compute (Ψ 0 ) * * : first, for q ∈ R k×d ,

(Ψ 0 ) * (q) = sup 1≤i<j≤k p∈R k (q i -q j ) • p -|p| = max 1≤i<j≤k    0 if |q i -q j | ≤ 1 , +∞ else.
Letting therefore

K = q = (q 1 , . . . , q k ) T ∈ R k×d : |q i -q j | ≤ 1 ∀i < j (29)
we find that (Ψ 0 ) * * is the support function of K, that is, (Ψ 0 ) * * = sup q∈K q • p. We have shown the following:

Proposition 3.3. The largest convex, local functional of the form (13), with Ψ(x, •) a nonnegative, even convex function which satisfies both ( 27) and (28), is

J(v) = Ω Ψ(Dv)
where Ψ(p) = sup q∈K q • p, K given by [START_REF] Giusti | Minimal surfaces and functions of bounded variation[END_REF].

Denote now, for all i, j:

K i,j = {q : |q i -q j | ≤ 1} so that K = i<j K i,j . Then, defining Ψ i,j (q) = sup q∈K i,j q • p =    |p| if p = (f i -f j ) ⊗ p , p ∈ R d +∞ else,
one can express Ψ as the inf-convolution of the Ψ i,j 's, that is

Ψ(p) = min i,j p i,j =p i<j Ψ i,j (p i,j )
which is also

Ψ(p) = min p i,j ∈R k p i = j>i p i,j -j<i p j,i i<j |p i,j | . ( 30
)

The convex partition problem

Our approach for computing minimizers of (1) (equivalently, ( 18)) will hence be the following: we will look for a minimizer v of

min v∈S J(v) + Ω v(x) • g(x) dx . (31)
Then, two situations may occur:

1. v ∈ S 0 , that is, v(x) ∈ {0, 1} k a.e. with k i=1 vi (x) = 1 a.e. in Ω: in this case, we have J(v) = J (v), and since J ≤ J we have solved [START_REF] Chambolle | Theoretical Foundations and Numerical Methods for Sparse Recov ery[END_REF] and solved the partition problem (with field g).

2. v ∈ S 0 : this covers, in fact, two situations quite different in nature:

(a) Either J(v) + Ω g • v dx = inf v J (v) + Ω g • v dx: in this case, one can show that
v is a convex combination of minimizers of [START_REF] Chambolle | Theoretical Foundations and Numerical Methods for Sparse Recov ery[END_REF], which are non-unique. This is tricky, since except for k = 2 (where a thresholding will do) there is no clear way to find back a partition from the values of v. The best way might be to perturb a little the field g, hoping to fall back in a case of uniqueness. This is illustrated in Fig. 7.

(b) Or J(v) + Ω g • v dx < inf v J (v) + Ω g • v dx:
in which case there seems to be no hope to find the optimal partition from the result. This is the case in the example of Fig. 6.

What is a maybe surprizing is that we found many more examples which fall in the first of these two categories than in the second. For relaxations which are less tight such as ( 48) and ( 49) -see Sec. 5 -this happens less frequently. In the second case, the relaxed solution needs to be binarized -see [START_REF] Lellmann | Continuous multiclass labeling approaches and algorithms[END_REF] for a discussion of different binarization strategies, and [START_REF] Lellmann | Optimality bounds for a variational relaxation of the image partitioning problem[END_REF] for interesting bounds on the energy of the binarized solution.

Before explaining how ( 31) is practically solved, let us analyse a little further the properties of the local convex envelope J.

Estimates for Ψ

We check here that Ψ satisfies basic estimates which show that, in particular, ( 22) is nearly optimal.

Lemma 3.4. For all p = (p 1 , . . . , p k ) T ∈ R k×d : if k i=1 p i = 0, then Ψ(p) = +∞. Otherwise, one has:

1 2 k i=1 |p i | ≤ Ψ(p) ≤ k -1 k k i=1 |p i | . (32)
In particular, (22) holds with J * * replaced with J and for any v

∈ L 2 (Ω; R k ) with k i=1 v i = 1 a.e. in Ω.
Proof. The first inequality in [START_REF] Ishikawa | Segmentation by grouping junctions[END_REF] clearly follows from our construction. It is also a straightforward consequence of the inclusion B(0, 1/2) k ⊂ K. On the other hand, if q ∈ K, then q + (q, . . . , q) T ∈ K for any q ∈ R d , and we deduce that

Ψ(p) = sup q∈K q • p = sup q∈K,q∈R d q • p + q • k i=1 p i .
Hence, if k i=1 p i = 0, this is +∞, as claimed. On the other hand, if k i=1 p i = 0, we see that the sup can be taken only on the vectors q ∈ K with i q i = 0 (simply choose q = -(1/k) k i=1 q i ). But in this case, for any i = 1, . . . , k,

q i = - j =i q j = 1 - 1 k q i - 1 k j =i q j = 1 k j =i (q i -q j )
and we deduce |q i | ≤ (k -1)/k. The right-hand side of (32) follows.

Duality. Paired calibrations

The functional J(v) = Ω Ψ(Dv) has the dual representation [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF]:

J(v) = sup - Ω k i=1 v i div ξ i : ξ ∈ C ∞ c (Ω; R k×d ) , ξ(x) ∈ K ∀ x ∈ Ω .
Consider a Dirichlet minimimal partition problem: we assume that Ω is bounded with Lipschitz boundary, and we choose a partition (E 0 i ) k i=1 of ∂Ω. We then look for a minimal partition of Ω, with this given trace, that is, which minimizes min

(E i ) k i=1 1 2 k i=1 Per(E i , Ω) (where (E i ) k
i=1 is a partition) under the constraint that E i meets ∂Ω on E 0 i . In fact, it is standard that this constraint should be relaxed in the following way: min

(E i ) k i=1 1 2 k i=1 Per(E i , Ω) + ∂Ω |χ E 0 i -χ E i | . We let v 0 = (χ E 0 1 , . . . , χ E 0 i ): then the local convexification of this problems becomes min v∈S Ω Ψ(Dv) + ∂Ω Ψ((v 0 -v) ⊗ ➉ ♥❻ Ω (x)) dH d-1 (33)
where ➉ ♥❻ Ω is the outer normal to ∂Ω. It follows then from Corollary 6.7 in Section 6.6 that the constraint v ∈ S can be removed in [START_REF] Ising | Beitrag zur Theorie des Ferromagnetismus[END_REF]: a minimizer over v ∈ L 2 (Ω; R k ) will in fact belong to S. We deduce the following result:

Proposition 3.5. The "generalized partition" v minimizes (33) if and only if there exists a vector field ξ ∈ L ∞ (Ω; R k×d ), with div ξ = 0 in Ω, |ξ iξ j | ≤ 1 a.e., and such that

k i=1 ξ i • Dv i = Ψ(Dv) |Dv| -a.e. in Ω , (34)
k i=1 (v 0 i -v i )ξ i • ➉ ♥❻ Ω = Ψ((v 0 -v) ⊗ ➉ ♥❻ Ω (x))
a.e. on ∂Ω

Proof. An essential point is Corollary 6.7 in Section 6.6, which shows that the minimum of [START_REF] Ising | Beitrag zur Theorie des Ferromagnetismus[END_REF] The functional in [START_REF] Ising | Beitrag zur Theorie des Ferromagnetismus[END_REF], which we denote by J, is a convex, l.s.c. function from L 2 (Ω; R k ) to [0, +∞]. If v is a minimizer (which exists thanks to Rellich's compactness theorem), over L 2 (Ω; R k ), we must have 0 ∈ ∂J(v), the subgradient of J at v defined in a standard way by

among all v ∈ L 2 (Ω; R k ) is, in fact, reached for a field v ∈ BV (Ω; R k ) with v ∈ S, that is, v(x) ∈ S (
w ∈ L 2 (Ω; R k×d ) : J(v) + Ω w • (v ′ -v) dx ≤ J(v ′ ) ∀ v ′ ∈ L 2 (Ω; R k×d ) .
This subgradient is described in [START_REF] Moll | The anisotropic total variation flow[END_REF][START_REF] Andreu | The Dirichlet problem for the total variation flow[END_REF] as follows: it consists in the functions w ∈ L 2 (Ω; R k ) such that there exists a vector fields ξ ∈ L ∞ (Ω; R k×d ) with ξ ∈ K and w = -div ξ a.e., and which satisfy the conditions [START_REF] Lai | The convergence of a central-difference discretization of Rudin-Osher-Fatemi model for image denoising[END_REF] and [START_REF] Lawlor | Paired calibrations applied to soap films, immiscible fluids, and surfaces or networks minimizing other norms[END_REF]. If v is itself a partition (i.e., v ∈ S 0 ), then |Dv| is purely singular and carried by the jump set J v which is the boundary of the partition. Then the first condition expresses the fact that the flux of ξ through this boundary is equal to the energy. In general, it is equivalent to requiring that

Ω k i=1 ξ i (x) • Dv i (x) = J(v) ,
while "≤" is always true because ξ ∈ K a.e. in Ω. The second condition [START_REF] Lawlor | Paired calibrations applied to soap films, immiscible fluids, and surfaces or networks minimizing other norms[END_REF] does not enforce any constraint as long as v = v 0 on ∂Ω. If not, then the field ξ is determined by the difference v 0 -v and the normal ➉ ♥❻ Ω to ∂Ω. The fact that these conditions are sufficient for w = -div ξ for being a subgradient at v is easily checked by integration by parts, the necessity for w to be of this form is more difficult.

The proposition follows from the two previous observations.

A vector field which satisfies the conditions in Proposition 3.5 is called a "paired calibration", following a terminology of Lawlor and Morgan [START_REF] Lawlor | Paired calibrations applied to soap films, immiscible fluids, and surfaces or networks minimizing other norms[END_REF]. They and Brakke [START_REF] Brakke | Soap films and covering spaces[END_REF] observed first that the existence of a paired calibration was a sufficient condition for the minimality of a partition. This also follows from Proposition 3.5, since if we have a function v ∈ S 0 (i.e., a partition) for which a paired calibration, satisfying all the conditions in the Proposition. exists, then it is a minimizer of (33), hence a minimal partition.

Let us consider now very simple example in dimension two. It is known that the boundary of a minimal partition, in this case, is made of straight lines which meet at triple points (a "Taylor junction"), with an angle of 120 • [START_REF] Bellettini | An almost everywhere regularity result for minimal partitions[END_REF][START_REF] Bonnet | On the regularity of edges in image segmentation[END_REF]. The following very easy proposition is standard and its proof is found for instance in [START_REF] Brakke | Soap films and covering spaces[END_REF] Proposition 3.6. Let Ω ⊂ R 2 be a Lipschitz, connected set, with 0 ∈ Ω.

Let Ē2 = {x = (x 1 , x 2 ) ∈ Ω : x 2 > 0 , x 1 > - √ 3x 2 }, Ē3 = {x = (x 1 , x 2 ) ∈ Ω : x 2 < 0 , x 1 < √ 3x 2 }, Ē1 = Ω \ ( Ē1 ∪ Ē2 ). Then v = (χ E 0 , χ E 1 , χ E 2 ) is a solution of min v∈BV (Ω;R k ) {J(v) : v = v on ∂Ω} .
This means that the minimum of the local convex envelope J is, in fact, the minimal partition (which solves the partition functional J ). The proof is straightforward, as there exists a constant vector field ξ ∈ K (hence with div ξ = 0), which satisfies both [START_REF] Lai | The convergence of a central-difference discretization of Rudin-Osher-Fatemi model for image denoising[END_REF] and [START_REF] Lawlor | Paired calibrations applied to soap films, immiscible fluids, and surfaces or networks minimizing other norms[END_REF].

(1, 0, 0) (0, 0, 1) Just pick ξ 1 = (0, 0), ξ 2 = ( √ 3/2, 1/2), and

ξ 3 -ξ 2 (0, 1, 0) ξ 3 -ξ 1 ξ 2 -ξ 1
ξ 3 = ( √ 3/2, -1/2). Observe that ξ 2 -ξ 3 = (0, 1), see Fig. 3. Then, ξ ∈ K since |ξ 1 -ξ 2 | = |ξ 1 -ξ 3 | = |ξ 2 -ξ 3 | = 1, (ξ 2 -ξ 1 ) = ν E 2 on ∂E 2 ∩ ∂E 1 , (ξ 2 -ξ 3 ) = ν E 2 on ∂E 2 ∩ ∂E 3 , and (ξ 3 -ξ 1 ) = ν E 3 on ∂E 3 ∩ ∂E 1 : it follows that (34) holds.
The second condition is automatically satisfied since v is of course equal to its own trace on ∂Ω. We will check that "looser" convex approximations of the partition problem fail at reconstructing properly the triple point, see Section 5.

Remark 3.7. Although we have performed very few simulations in three dimensions, the structure of minimal surfaces in 3D and their singularities, as described by J. Taylor [START_REF] Taylor | The structure of singularities in soap-bubble-like and soap-film-like minimal surfaces[END_REF], is such that for all three types of singular minimal cones which she describes, there also will exist a calibration made of constant vector fields as in the 2D situation. We mention that two of these minimal cones only can appear in the three phases problem (a plane and three planes meeting on one line, that is the 2D triple point extended by translation in the orthogonal direction), while the third one (built upon a regular tetraedron) involves at least four phases. The first two cones inherit their calibration from the 2D situation, while it is observed in [START_REF] Lawlor | Paired calibrations applied to soap films, immiscible fluids, and surfaces or networks minimizing other norms[END_REF] that in any dimension, there is an obvious constant calibration for the cone built upon the (d -2)-dimensional edges of the regular simplex, which is made of the (appropriately resized) normal vectors to each facet of the simplex. For more details see Section 4.3.3.

Numerical approximation 4.1 Discretization

There are many ways to numerically tackle problems [START_REF] Ishikawa | Exact optimization for Markov random fields with convex prio rs[END_REF] and [START_REF] Ising | Beitrag zur Theorie des Ferromagnetismus[END_REF]. We present a simple approach which is quite efficient, in the framework of finite differences.

Approximation To simplify we consider only the 2D case (d = 2), and we assume that Ω = (0, 1) 2 . Let N > 1 and h = 1/N > 0 be a discretization step.

To allievate the notation, let us first consider the discretization of the simplified, scalar problem min

v∈BV (Ω;[0,1] Ω |Dv| + Ω g(x)v(x) dx (36)
If

v h = (v h i,j ) 0≤i,j≤N -1 ∈ X h ∼ R N ×N is a discrete function, which is identified with the function in L 2 (Ω) v h (x) = 0≤i,j<N v h i,j χ (ih,jh)+[0,h) 2 (x) , (37)
we can define its discrete total variation by

T V h (v h ) = h 2 0≤i,j<N (∇ h v h ) i,j (38)
where | • | is the standard Euclidean norm and ∇ h :

X h → X h × X h is defined by (∇ h v h ) i,j = 1 h              (v h i+1,j -v h i,j , v h i,j+1 -v h i,j ) T if 0 ≤ i, j < N -1 , (v h i+1,j -v h i,j , 0) T if 0 ≤ i < N -1, j = N -1 , (0, v h i,j+1 -v h i,j ) T if i = N -1, 0 ≤ j < N -1 , (0, 0) T if (i, j) = (N -1, N -1)
Then, assuming that T V h is extended to L 2 (Ω) by letting T V h (v) = +∞ whenever v is not of the form (37), we have the following elementary result. For the definition and properties of Γ-convergence, see [START_REF] Braides | Γ-convergence for beginners[END_REF][START_REF] Maso | An introduction to Γ-convergence[END_REF], and Remark 4.2.

Proposition 4.1. T V h Γ-converges to the total variation Ω |Dv| as h → 0.

Proof. Let just quickly sketch a proof of this result. By definition of the Γ-convergence we need to show that [START_REF] Braides | Γ-convergence for beginners[END_REF] for any v ∈ BV (Ω):

(i) If v h → v, then Ω |Dv| ≤ lim inf h→0 T V h (v h ); (ii) There exists v h → v with Ω |Dv| ≥ lim sup h→0 T V h (v h ).
The proof of (ii) is standard and based on the fact that any v ∈ BV (Ω) can be approximated with smooth functions

v n ∈ C ∞ (Ω) with Ω |∇v n | dx → Ω |Dv|, see Theorem 2.1. Then, it is easy to show that if v n ∈ C ∞ (Ω), by simply defining v n,h i,j = v n (ih, jh) we have lim h→0 T V h (v n,h ) = Ω |∇v n | dx (
with, in this case, an error bound depending on D 2 v n ). Point (ii) follows from a standard diagonal argument.

The proof of (i) is based on the dual formulation of the total variation. We consider (v h ) h , a sequence of discrete functions at scale h, h = 1/N → 0, which converges to v in L 2 (Ω). To show (i), it is enough to show that for any ξ ∈ C ∞ c (Ω; R 2 ) with |ξ| ≤ 1 in Ω, we have

- Ω v(x)div ξ(x) dx ≤ lim inf h→0 T V h (v h ). (39)
Then, taking the supremum over all such ξ's will yield (i). We have

- Ω v(x)div ξ(x) dx = lim h→0 - Ω v h (x)div ξ(x) dx = lim h→0 i,j h 0 ξ 1 ((i + 1)h, jh + s)(v h i+1,j -v h i,j ) ds + h 0 ξ 2 (ih + s, (j + 1)h)(v h i,j+1 -v h i,j ) ds . (40) Now, if s ∈ (0, h), |((i + 1)h, jh + s) -(ih + s, (j + 1)h)| = |(h -s, s -h)| ≤ √ 2h, so that ξ 1 ((i + 1)h, jh + s) 2 + ξ 2 (ih + s, (j + 1)h) 2 ≤ |ξ((i + 1)h, jh + s)| + √ 2 ∇ξ L ∞ (Ω) h ≤ 1 + Ch
and it follows that

ξ 1 ((i + 1)h, jh + s)(v h i+1,j -v h i,j ) + ξ 2 (ih + s, (j + 1)h)(v h i,j+1 -v h i,j ) ≤ (1 + Ch)h (∇ h v h ) i,j
We deduce from ( 40) that (39) holds, hence (i).

Remark 4.2. By standard properties of the Γ-convergence [START_REF] Braides | Γ-convergence for beginners[END_REF], it follows that if v h is a minimizer of

J h (v h ) + Ω g(x)v h (x) dx , (41)
then as h → 0, v h converges to a solution of (36) (up to subsequences if this solution is non unique). Hence, to solve our problem, we will minimize [START_REF] Moll | The anisotropic total variation flow[END_REF]. Observe that Propostion 4.1 does not provide any error estimate. A recent approach to error estimates for this type of problems is found in [START_REF] Lai | The convergence of a central-difference discretization of Rudin-Osher-Fatemi model for image denoising[END_REF].

The vectorial case is identical: now we consider at scale h = 1/N > 0, the discrete version of J, defined for

v h = ((v h 1 ) i,j , . . . , (v h k ) i,j ) 0≤i,j≤N -1 ∈ (X h ) k : J h (v h ) = h 2 0≤i,j<N Ψ (∇ h v h l ) i,j ) k l=1 ( 42
)
where Ψ is still the integrand in Prop. 3.3. Then, J h is considered as a functional in L 2 (Ω; R k ), v h being identified with the vectorial function

v h (x) = 0≤i,j<N v h i,j χ (ih,jh)+[0,h) 2 (x) , (43)
as previously in [START_REF] Lellmann | Optimality bounds for a variational relaxation of the image partitioning problem[END_REF]. The proof of the following result is then identical to the proof Proposition 4.1 in the scalar setting, provided we use in addition to the approximation Theorem 2.1, the convergence ( 16):

Proposition 4.3. As h → 0, J h Γ-converges to J(v) = Ω Ψ(Dv) (and J(v) = +∞ whenever v ∈ L 2 (Ω; R k ) \ BV (Ω; R k ))
This means that, as in Remark (4.2), we can approximate the minimization of (31) with the discrete minimization of min

v h l ≥0 , k l=1 v l h =1 J h (u h ) + Ω g(x) • v h (x) dx. (44)
The approximation of Dirichlet problems such as [START_REF] Ising | Beitrag zur Theorie des Ferromagnetismus[END_REF] follows the same lines.

A simple, converging algorithm

Now, how do we effectively solve (44)? There are essentially two difficulties. One comes from the fact that there is no tractable, simple expression of Ψ (except for (30) which could actually be used in an iterative way), so that it is more reasonable to use its dual definition (as the support function of K) in the optimisation. The second comes from the fact that (at least if k ≥ 3) the problem involves quite a lot of variables and it seems better to design a simple algorithm, especially if it is so simple that it can be implemented on basic parallel processing units such as GPUs (graphic cards). For these reasons, we propose here to use a primal-dual Arrow-Hurwicz [START_REF] Arrow | Studies in linear and non-linear programming[END_REF][START_REF] Ekeland | Convex analysis and variational problems[END_REF] type algorithm for minimizing [START_REF] Moll | The anisotropic total variation flow[END_REF]. The advantage of this algorithm is that it extends easily to the problems we will address in the next sections. It has been recently suggested in this framework (more precisely, total variation minimization for image denoisi¡ng and reconstruction), first by H. Talbot and B. Appleton [START_REF] Appleton | Globally optimal geodesic active contours[END_REF] (who, in fact, suggested a primal-dual flow in a continuous setting, which was then discretized using classical methods for hyperbolic schemes) and more recently, in a setting closer to ours, by M. Zhu and T. Chan [START_REF] Zhu | An efficient primal-dual hybrid gradient algorithm for total variation image restoration[END_REF]. The version we propose (first in [START_REF] Pock | A convex relaxation approach for computing minimal partitions[END_REF], with a first proof of convergence), which is slightly different, is improved in the sense that we can provide an estimate for the error in the objective [START_REF] Chambolle | Theoretical Foundations and Numerical Methods for Sparse Recov ery[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems withapplications to imaging[END_REF]. It is in fact a variant of the Douglas-Rachford splitting method [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF], and is inspired from similar algorithms in [START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF][START_REF] Popov | A modification of the arrow-hurwicz method for search of saddl e points[END_REF]. (We refer to [START_REF] Pock | A convex relaxation approach for computing minimal partitions[END_REF][START_REF] Chambolle | Theoretical Foundations and Numerical Methods for Sparse Recov ery[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems withapplications to imaging[END_REF] for details.)

We fix a scale h = 1/N > 0, and first discretize the external field by letting for l = 1, . . . , k

(G h l ) i,j = 1 h 2 (ih,jh)+[0,h) 2 g l (x) dx
(if not already given in a discrete form). We denote G h i,j = ((G h 1 ) i,j , . . . , (G h k ) i,j ). Then, we introduce a primal variable V = (V i,j ) ∈ (X h ) k and a dual Ξ = (Ξ i,j ) ∈ (X h × X h ) k . The problem we want to optimize may be written as min

V i,j ∈S max ((Ξ l ) i,j ) k l=1 ∈K h 2 i,j k l=1 (Ξ l ) i,j • (∇ h V l ) i,j + h 2 i,j V i,j • G h i,j , (45)
where we recall that the simplex S = {z ∈ R k + :

k i=1 z i = 1}.
The idea of a Arrow-Hurwicz algorithm is to follow simultaneously a gradient descent in V and ascent in Ξ, until convergence to the saddle-point.

First we initialize with Ξ 0 = 0, V 0 = V 0 = 0. Then, we fix two "time steps" τ, τ ′ > 0 and we update Ξ n , V n , V n by letting

(Ξ n+1 ) i,j = Π K Ξ n i,j + τ ′ (∇ h V n ) i,j
for all i, j

V n+1 i,j = Π S V n i,j + τ (div h Ξ n+1 ) i,j -G h i,j
for all i, j

V n+1 = 2V n+1 -V n . (46)
Here, the discrete divergence div h = -(∇ h ) * is the opposite of the adjoint of the discrete gradient. The projections Π K and Π S are respectively the projections on the convex sets K and S. The first is quite complicated to perform, and is done using Dykstra's iterative algorithm, see Appendix A. It motivates the use of a GPU since it can be easily parallelized. The projection onto the simplex S can be done in linear time. See for example [START_REF] Duchi | Efficient projections onto the l1-ball for learning in high dimensions[END_REF].

It is shown that the scheme converges as long as τ τ ′ ≤ 1/8 (for h = 1), and in some sense the speed of convergence is optimal [START_REF] Chambolle | Theoretical Foundations and Numerical Methods for Sparse Recov ery[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems withapplications to imaging[END_REF]. We found that it is a good compromise between efficiency and memory storage (most variables need be stored just once, the projection over K is done only once per iteration).

A theoretical advantage of this approach is that the primal-dual gap, given by max

Ξ i,j ∈S h 2 i,j Ξ i,j • (∇ h V n ) i,j + h 2 i,j V n i,j • G h i,j -min V i,j ∈S h 2 i,j (-div h Ξ n i,j ) • V i,j + h 2 i,j V i,j • G h i,j = J h (V n ) + h 2 i,j V n i,j • G h i,j + h 2 i,j max l=1,...,k [(div h Ξ n l ) i,j -(G h l ) i,j ] ≥ 0 (47)
may be computed at each step and goes to zero at convergence, and could be taken as a criterion for convergence. However, for k ≥ 4, the calculation of the primal energy J h (that is, of Ψ), is nearly intractable. For k = 3, it can be done quite efficiently, see Appendix C.

Examples

Problems with 3 labels

To test this approach, we have implemented a very simple program where the input is a color image and the weight on each set E 1 , E 2 , E 3 is equal to minus the level of red, green and blue respectively. The result is a kind of projection to the closest of these three colors (of course, this gives terrible output for most images). The example on Figure 4, right, is an output for an input image with areas dominantly red, green or blue. The example corresponds to a value of h = 1 (the discretization step in ( 44)), and the weight is (-2)× the level of each channel, normalized between 0 and 1. The image has almost 120000 points (412 × 291). The

Input image

Output Figure 4: An example of a three-region segmentation.

gap [START_REF] Popov | A modification of the arrow-hurwicz method for search of saddl e points[END_REF] at convergence is ≤ 1 (for an energy of -122355). The total length energy (J h (v)) of the output is 5897. (Hence, "most" of the energy is the external field term.) The output is "almost" binary (i.e., the field v is most of the time an extreme point of the simples S).

In fact, we can compute a "width" of the interface, defined for instance as the ratio of the number of points where u 1 = 1v 1 and u 2 = u 1v 2 are between 0.1 and 0.9, over the total "length" J h (u). In this example, we found .75 pixel units, which is quite narrow. This is not expected to be zero, because the discretization (42) of J requires a fuzzy interface to approximate precisely the length (just as [START_REF] Lellmann | Continuous multiclass labeling approaches and algorithms[END_REF], this is clear if one thinks of how the Γ-limit superior is established, point (ii) in the proof of Proposition 4.1).

We have also tried to solve the Dirichlet problem [START_REF] Ising | Beitrag zur Theorie des Ferromagnetismus[END_REF]. A simple way to do (approximately) so is to minimize [START_REF] Nemirovski | Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems[END_REF] with a very strong external field (forcing pure red, green or blue) except in a grey area where no particular color (i.e., set E i , i = 1, 2, 3) is favored. The result is impressive: what is expected, that is, a sharp discontinuity set with a triple point where all three interfaces meet with and angle of 120 • , is actually computed by the program: see Figure 5. In this case, the weights g i are 1/20 in the red, green, blue area and equal in the Input with missing region Inpainted image grey area. This makes the length term quite important in the total energy, and actually the total length J h (v) which is computed is about 603 which is not much more of the "true" (isotropic) length of the total interface, estimated between 599 and 600 pixels. The size of the image is 376 × 357. This is of course an illustration of Proposition 3.6.

In fact, in trying to understand why it so often works so well, we could also find a "counterexample," that is, a situation where we do not expect a binary solution (v i ∈ {0, 1} a.e.), but rather a mixed one.

Why, sometimes, it does not work. Now, let us consider the same geometry, but instead of a boundary datum, a weight g 1 , g 2 , g 3 such that on each Ēi , g i = 1 while g j = 0 for j = i, except on a disk centered at the triple point where we choose g 1 = g 2 = g 3 . Equivalently, we run our program with as input, Ē1 colored pure cyan ((R, G, B) = (0, 1, 1)), Ē2 pure yellow (1, 1, 0), Ē3 pure magenta (1, 0, 1), except a grey area in the middle (see Figure 6, left). Then, we find that the grey circle is completed with a perfect triple point, but this time with mixed colors, that is, cyan ((R, G, B) = (0, .5, .5)) in Ē1 , yellow (.5, .5, 0) in Ē2 , magenta (.5, 0, .5) in Ē3 . In other words, the optimal solution is a mixture of half E j and half E k in Ēi , for {i, j, k} = {1, 2, 3}. The issue is that, now, the opposite of the "calibration" of Proposition 3.6 can show that this is actually a local minimizer inside the grey circle, as before (with expected energy of 3/2 times the radius of the circile). In particular, the length energy J h (v) which is computed now is 306, about half the true length, as expected, thus strictly below the length of any binary solution with the same pattern of discontinuity (and, actually, any reasonable binary solution for this input data). We should add that convergence of our scheme is quite slow (and poor) for this particular example. We really believe that the example in Figure 6 corresponds to a situation where J(v) < J * * (v), that is, the local convex envelope is strictly below the true (nonlocal) convex envelope.

Problems with more labels

We have observed that our approach works quite well with a quite high number of labels. Figure 7 shows the completion of four regions with equidistant labels Here the values of (a i ) 4 i=1 are kept fixed, and span the four color intensity values (1, 0, 0), (0, 1, 0), (0, 0, 1), (1, 1, 1). The result is what is theoretically expected (two triple junctions, with angles of 120 • ). However, we must point out that the solution which is found here is only one out of two. It seems that the program selects this particular solution because the discretization makes it minimal. If one rotates the original image by 45 • , then both solutions have the same energy even for the discrete problem, and our program produces an output which is not binary, but a convex combination of the two minimal binary solutions, as shown in Figure 8.

We insist that it does not correspond to a situation where the convexification is too low (as in Fig. 6), but just a case of non uniqueness, where any convex combination of the binary solutions is also minimizing, as mentioned in point 2a, p. 15.

We also have implemented a "basic" image segmentation model (following a piecewise constant Mumford-Shah model, see [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF][START_REF] Chan | Active contour and segmentation models using geometric PDE's for medical imaging[END_REF]). The idea is to solve alternatively (for fixed k)

min (E i ,a i ) k i=1 min λ 2 k i=1 Per(E i , Ω) + k i=1 E i (I(x) -a i ) 2 dx
with respect to (E i ) k i=1 and then to (a i ) k i=1 . Here, I : Ω → [0, 1] 3 is the color information (intensity of red, green, blue channels) of the original image. The initial values (a i ) k i=1 (also vector-valued) are initialized using the k-means algorithm, and then, once a new partition (E i ) is found, each a i is updated by computing the average value of I in E i . Figure 1 was computed in this way, while Figure 9 shows another example, once with 4 and then with 10 labels. Again, as for three labels, in almost all our experiment the solution was nearly binary (up to a smoothing of the discontinuities due to the discretization).

Input 4 regions 10 regions

Figure 9: Basic piecewise constant segmentation.

3D Problems

We tested our approach on a few standard test cases in 3D. The generalization of the numerical scheme to 3D is straight forward. For the minimal surface problems presented below we enforced the boundary conditions by generating appropriate external fiels G h l at the faces of a cubic grid.

Figure 10 shows the result of a minimal surface problem with three phases. Similar to the 2D case, our method finds an almost binary solution. Figure 11 shows the result of a minimal surface problem with six phases. Clearly, the solution of the problem with six phases is not unique. In order to obtain a unique solution, we introduced a slight anisotropy in the size of the grid cells. Again, our algorithm finds an almost binary solution.

Comparison with other relaxations

In this section we compare our approaches with the recent proposals of [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF] and [START_REF] Lellmann | Convex multi-class image labeling by simplex-constrained total variation[END_REF]. We show that our approach is tighter, and will solve problems which cannot be solved by looser relaxations. This is illustrated numerically in Fig. 13.

In [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF] the problem is simply relaxed as

min v 1 2 k i=1 Ω |Dv i | + k i=1 Ω g i (x)v i (x) dx , (48)
that is, replacing J * * with F. Clearly, if v i ∈ {0, 1} a.e., this is the same. Proposition 21 shows that this relaxation is, in general, below ours (except for k = 2 where they are the same). It should therefore fail more often.

A variant is proposed in [START_REF] Lellmann | Convex multi-class image labeling by simplex-constrained total variation[END_REF]: the mixed 1, 2-norm of (Dv 1 , . . . , Dv k ) used in ( 48) is replaced with the 2-norm of the matrix (that is, a "standard" classical TV, this time normalized with 1/ √ 2). Its gives the relaxation:

min v 1 √ 2 Ω |Dv| + k i=1 Ω g i (x)v i (x) dx , (49)
where Ω |Dv| is defined by [START_REF] Birchfield | A pixel dissimilarity measure that is insensitive to image sampling[END_REF]. Again, if v ∈ S 0 , then this boils down to energy [START_REF] Chambolle | Theoretical Foundations and Numerical Methods for Sparse Recov ery[END_REF]. However, in general, it is below. A first remark is that the integrand is of the form Ψ L (Dv), where Ψ L is the support function of the Euclidean ball K L of radius 1/ √ 2 in R k×d :

K L = q = (q 1 , . . . , q k ) T ∈ R k×d : k i=1 q 2 i ≤ 1 2 . Now, clearly, K L ⊂ K (defined in (29)), since if q ∈ K L , (q i -q j ) 2 = q 2 i + q 2 j -2q i • q j ≤ 2(q 2
i + q 2 j ) ≤ 1 for any i < j, so that q ∈ K. This shows that Ψ L ≤ Ψ. We can show the following, more precise result: Proposition 5.1. Let Ω, Ēi , i = 1, 2, 3 be as in Proposition 3.6. Let vi = χ Ēi for each i, and

v = (v 1 , v2 , v3 ). Then v is not a minimizer of F(v), nor of Ω Ψ L (Dv) = (1/ √ 2) Ω |Dv|, with prescribed boundary condition v = v on ∂Ω.
Remark 5.2. The interfacial energy which is actually approached in [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF] is anisotropic and, in fact, a sum of 1D energies in each direction. In this case, one can show again that their approach is equivalent to (the anisotropic version of) ours. We claim that the correct isotropic generalization of the approach in [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF] is what we presented in this paper.

Proof. To simplify the notation, in this proof, let the sets Ēi be defined in all R 2 (and vi = χ Ēi ∩Ω for each i). We first prove that v is not a minimizer of F with prescribed boundary condition. We consider afterwards relaxation [START_REF] Rockafellar | Integrals which are convex functionals[END_REF].

Step 1 First, let us assume that Ω is a convex set and let us not worry about the boundary condition. We choose a rotationally symmetric, smooth mollifier ρ ∈ C ∞ c (B(0, 1); R + ), with B(0,1) ρ dx = 1. As usual, for ε > 0, we let ρ ε (x) = ρ(x/ε)/ε 2 . We then define, for i = 1, 2, 3, a smooth function v ε i as the restriction to Ω of ρ ε * χ Ēi . By linearity, clearly, we still have

3 i=1 v ε i = 1.
By the co-area formula [START_REF] Brown | A convex relaxation method for a class of vector-valued minimization problems with applications to Mumford-Shah segmentation[END_REF], the variation Ω |Dv ε i | is the average of the lengths of the level lines {v ε i = s} for s ∈ (0, 1). Let us choose a coordinate system (y 1 , y 2 ) such that Ē1 = {y 2 < -|y 1 |/ √ 3}, and the jump set ∂ Ē1 is the graph y 2 = w(y 1 ) := -|y 1 |/ √ 3. Assume for a while, to simplify, that Ω, near this graph, coincides with the set |y 1 | < 2. First, assume also ε = 1 and consider the function

v 1 1 = ρ * χ Ē1 .
It is nondecreasing with respect to y 2 , and even strictly nondecreasing at points (y 1 , y 2 ) such that the support of ρ(• -(y 1 , y 2 )) meets ∂ Ē1 , that is, as soon as 0 < v 1 1 < 1. Invoking the implicit functions theorem, we see that for any 0 < s < 1, the level line {v 1 1 = s} is the graph of a smooth, even function w s , defined by v 1 1 (y 1 , w s (y 1 )) = s. In particular,

w ′ s (y 1 ) = -(∂ 1 v 1 1 )/(∂ 2 v 1 1 )(y 1 , w s (y 1 )
). Outside of B(0, 1) (hence, in particular, for |y 1 | > 1), we clearly have w s = w + c s for a constant c s ∈ (-2/ √ 3, 2/ √ 3), and with c s = -c 1-s for each s ∈ (0, 1). For each s, hence,

w ′ s (y 1 ) = w ′ (y 1 ) = ±1/ √ 3 if |y 1 | > 1.
We check that 1/ √ 3 is an upper bound for the derivative of w ′ s : indeed, since any translate of Ē1 in a direction (t, -1) with |t| ≤ √ 3 is included in Ē1 , we also have

t∂ 1 v 1 1 -∂ 2 v 1 1 ≥ 0: hence tw ′ s + 1 ≥ 0 for all t with |t| ≤ √ 3, that is, |w ′ s | ≤ 1/ √ 3 = |w ′ |.
Since, moreover, w ′ s (0) = 0 (w s is even), we deduce

H 1 ({v 1 1 = s}) = 2 -2 1 + |w ′ s | 2 dy 1 < 4 × 2 √ 3 = H 1 (∂ Ē1 ∩ Ω) .
Hence, from the co-area formula we deduce

Ω |Dv 1 1 | < H 1 (∂ Ē1 ∩ Ω).
A simple scaling argument will then show that there exists c > 0 such that

Ω |Dv ε 1 | ≤ H 1 (∂ Ē1 ∩ Ω) -cε for ε ≤ 1. Now, if
the boundary of ∂Ω is made of two straight lines (not necessarily vertical) in the neighborhood of its intersection with ∂ Ē1 , the same is still true, because the possible increase or loss of length of {v ε 1 = s} at the boundary (with respect to the previous case) is compensated exactly by the loss or increase of the symmetric line {v ε 1 = 1 -s}, so that the average length remains the same as when ∂Ω is vertical. If Ω is a generic convex set, then the length of the lines {v ε 1 = s} in Ω are even shorter than if ∂Ω were replaced with two tangent straight lines at its intersection with ∂ Ē1 , so that the variation of v ε 1 is even lower.

Step 2 Now, we consider any open set Ω and η > 0 such that B(0, 2η) ⊂ Ω. We choose ε << η and let as above, for each i ∈ {1, 2, 3},

v ε i = ρ ε * χ Ēi in B(0, η), while v ε i = vi = χ Ēi in Ω \ B(0, 2η). In B(0, 2η) \ B(0, η), we build v ε
i by joining with straight lines the level lines of v ε i inside B(0, η) and outside B(0, 2η): one can check that the total variation of each v ε i in the crown is then of order at most 2η + ε 2 /η as ε → 0, while, using the first step, it is less than 2ηcε inside the ball B(0, η). Hence, if ε is small enough, we find that Ω |Dv ε i | < Ω |Dv i |, and this shows the first part of the proposition.

Next we show why v is not even a minimizer of Ω Ψ L (Dv) = (1/ √ 2) Ω |Dv|, with v = v on ∂Ω. This is a bit more tricky, since we did not succeed in building explicitly a competitor: it therefore relies on a (simple) calibration argument. Let R > 0 be a radius such that B(0, R) ⊂ Ω. We first observe that for any r ∈ (0, R), we can build a competitor v r with same energy as v: we simply introduce the equilateral triangle T r which has vertices at (r, 0), (-r/2, r √ 3/2), (-r/2, -r √ 3/2). Outside of this triangle, we let v r = v. Inside the triangle, we let v r = (1/3, 1/3, 1/3). We have

Ω |Dv r | = Ω\B R |Dv| + 3 √ 2(R -r) + 2 3 2 + 1 3 2 + 1 3 2 H 1 (∂T r ) = Ω\B R |Dv| + 3 √ 2(R -r) + 2 3 × 3 √ 3r = Ω\B R |Dv| + 3 √ 2R = Ω |Dv| (50)
so that the energy of v r is the same as the energy of the optimal partition v. If this were the minimal energy, then, reasoning as in the proof of Proposition 3.5 (and admitting that the (1-Lipschitz) projection onto the simplex S reduces the energy, so that the unconstrained minimizer of Ω |Dv| with boundary condition v = v on ∂Ω is in fact in S), we find that there should exist a vector field ξ ∈ L ∞ (R k×d ), with zero divergence, ξ ∈ K L a.e., and such that both [START_REF] Lai | The convergence of a central-difference discretization of Rudin-Osher-Fatemi model for image denoising[END_REF] and [START_REF] Lawlor | Paired calibrations applied to soap films, immiscible fluids, and surfaces or networks minimizing other norms[END_REF] hold with Ψ replaced with Ψ L , and v with any minimizer, such as v, v r , or any convex combination of these (for instance, (1/R) R 0 v r dr, which is Lipschitz with piecewise constant gradient in B(0, R)). This fully determines ξ in B(0, R) (see Fig. 12):

we must have ξ = ([v r ]/|[v r ]|) ⊗ ν Tr (the normal to ∂T r ). Hence, a.e. in E 1 ∩ B(0, R), ξ = (2/ √ 6, -1/ √ 6, -1/ √ 6) ⊗ (-1, 0), in E 2 ∩ B(0, R), ξ = (-1/ √ 6, 2/ √ 6, -1/ √ 6) ⊗ (1/2, 3/2), in E 3 ∩ B(0, R), ξ = (-1/ √ 6, -1/ √ 6, 2/ √ 6) ⊗ (1/2, - √ 3/2).
(1, 0, 0) (0, 0, 1) (0, 1, 0)

( 1 3 , 1 3 , 1 3 )
min {E i } k i=1 1≤i<j≤k σ i,j H d-1 (∂E i ∩ ∂E j ) + k i=1 E i g i (x) dx (51)
where the positive weights σ i,j = σ j,i satisfy

σ i,j ≤ σ i,l + σ l,j (52)
for any {i, j, l} ⊂ {1, . . . , k}. This condition is necessary for [START_REF] Strekalovskiy | Generalized ordering constraints for multilabel optimization[END_REF] to have a solution (otherwise, the surface energy is not lower-semicontinuous, since it will better to insert between E i and E j an infinitesimal layer of label l), however, if it is not satisfied our approach will automatically solve a relaxed (l.s.c. envelope) version of the problem.

To solve this, we introduce as before, for v ∈ S 0 (defined in [START_REF] Chambolle | A convex approach for computing minimal partitions[END_REF], v = (χ E 1 , . . . , χ E k )),

J (v) = Jv σ 0 (v + (x) -v -(x)) dH d-1 (x)
where σ 0 (q) = σ i,j if q i = -q j = ±11 , q l = 0 for l = i, j. We would like to compute the convex envelope J * * of J , just like in Section 3, however, it is not expected to have a tractable expression. So that we introduce again the "local" convex envelope J(v) = Ω Ψ(Dv), where now Ψ is the support function of the closed, convex set

K = q = (q 1 , . . . , q k ) T ∈ R k×d : |q i -q j | ≤ σ i,j ∀i < j . (53)
(Observe that if [START_REF] Strekalovskiy | Total variation for cyclic structures: Convex relaxation and efficient minimization[END_REF] is not met, then |q iq l | ≤ σ i,l and |q lq j | ≤ σ l,j will yield |q iq j | ≤ σ i,l + σ l,j < σ i,j , so K does not depend on the latter value.)

As before, we propose to solve instead of (51), its variant [START_REF] Ishikawa | Exact optimization for Markov random fields with convex prio rs[END_REF], with the new function J.

An equivalent representation

Now, in some cases, when the labels are in an ordered set, with a meaningful value (e.g., the disparity in a stereo reconstruction problem), the representation v might not be adapted.

In [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF], where the interaction potential which is considered is the total variation of the labels, the variable which is used is the variable u i = χ {ι>i} where ι(x) ∈ {1, . . . k} is the label (which of course might be a different ordered set, but for simplicity we will stick to this case). It is obtained from v i = χ {ι=i} by the change of variable [START_REF] Bae | Global minimization for continuous multiphase partitioning problems using a dual approach[END_REF], and satisfies [START_REF] Bellettini | An almost everywhere regularity result for minimal partitions[END_REF]. Observe that to represent k labels, one needs only the k -1 variables u 1 , . . . , u k-1 , since u 0 ≡ 1 and u k ≡ 0 are fixed. Problem [START_REF] Ishikawa | Exact optimization for Markov random fields with convex prio rs[END_REF] can be rewritten in term of u = (u 1 , . . . , u k-1 ) ∈ BV (Ω; R k-1 ). First, we use the change of variable [START_REF] Bellettini | Convex approximations of functionals with curvature[END_REF] to introduce (with u 0 ≡ 1, u 1 ≡ 0)

H(u) = J((u i-1 -u i ) k i=1 ) = Ω Ψ((D(u i-1 -u i )) k i=1 ) = Ω Φ(Du)
where, for p = (p 1 , . . . , p k-1 ) (and letting p 0 = p k = 0)

Φ(p) = Ψ((p i-1 -p i ) k i=1 ) = sup q∈K k i=1 q i • (p i-1 -p i ) = sup q∈K k-1 i=1 (q i+1 -q i ) • p i = sup r∈L r • p (54)
where L = r = (q i+1q i ) k-1 i=1 ∈ R (k-1)×d : q ∈ K .

If r i = q i+1q i , we have q i = i-1 l=1 r l + q 1 for i ≥ 2 and if i < j, q jq i = j-1 l=i r l so that it follows from (53) that L = r = (r 1 , . . . , r k-1 ) ∈ R (k-1)×d :

j-1 l=i r l ≤ σ i,j ∀ i < j (55)
On the other hand, the term

Ω v • g dx = Ω k i=1 (u i-1 -u i )g i dx = Ω k-1 i=1 u i (g i+1 -g i ) dx + Ω g 1 dx = Ω u • g dx + const.
where g ∈ L 1 (Ω; R k-1 ) is given by gi (x) = g i+1 (x)g i (x), and we have used the convention u 0 ≡ 1, u k ≡ 0. We see that [START_REF] Ishikawa | Exact optimization for Markov random fields with convex prio rs[END_REF] is equivalent to the minimization min u∈C 0

H(u) + Ω u(x) • g(x) dx (56)
where

C 0 = u ∈ BV (Ω; R k-1 ) : 1 ≥ u 1 ≥ • • • ≥ u k-1 ≥ 0 a.e. . (57)
In [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF], the authors use σ i,j = |i -j| (the total variation of the labels): in this case, one checks easily that the convex L is simply given by

L = r = (r 1 , . . . , r k-1 ) ∈ R (k-1)×d : |r i | ≤ 1 ∀ i (58)
on which it is straightforward to project. On the other hand, the representation with the variable v would require to project onto K = {q : |q i+1q i | ≤ 1}, which is quite tricky. This explains why this representation is also useful, depending on the particular form of the surface tension σ i,j between the phases i and j.

Another case which is better represented in terms of the variable u is the "truncated total variation", σ i,j = min{|i -j|, T } for some threshold T > 0, see Fig. 2, right. In that case, L = r = (r 1 , . . . , r k-1 ) ∈ R (k-1)×d :

|r i | ≤ 1 ∀ i , and j-1 l=i r l ≤ T ∀ i < j with j -i ≥ T . (59)
The projection onto the convex set L given by ( 59) is clearly more complicated than onto [START_REF] Ziemer | Weakly differentiable functions[END_REF], however, it is still simpler than on the corresponding convex set K.

A simpler variant

In this section and all that follows, we assume that

σ i,j = σ(|j -i|)
for a function σ : (0, +∞) → (0, +∞) which is concave, positive, nondecreasing. This clearly implies [START_REF] Strekalovskiy | Total variation for cyclic structures: Convex relaxation and efficient minimization[END_REF] , which is our claim. All the particular energies considered up to now meet these conditions. Solving (56) with the numerical approach described in Section 4.2 is quite simple, but still requires to perform two relatively difficult tasks, which are (i) to project a vector onto L (ii) to project a vector onto

C 0 = {u ∈ R k-1 : 1 ≥ u 1 ≥ • • • ≥ u k-1 ≥ 0} (
which is in general much easier than the first projection, except in case L is given by ( 58)). It is suggested in [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF] to avoid this last step by minimizing, instead of (56), the unconstrained variant min

u∈BV (Ω;R k-1 ) H(u) + Ω k i=1 g i (x)|u i-1 (x) -u i (x)| dx (60)
(where still, u 0 ≡ 1 and u k ≡ 0). We can show the following proposition: Proposition 6.1. Assume g i > 0 a.e, and for each i = 1, . . . , k. Then, any solution of (60)

satisfies 0 ≤ u k-1 ≤ u k-2 ≤ • • • ≤ u 1 ≤ 1 a.e.
in Ω, and hence is a solution of (56)

a.e. x ∈ Ω. Indeed, at a.e. x, and for any α ∈ {1, . . . , d},

|∂ α Π C (u)(x)| = lim ε→0 1 ε |Π C (u)(x + εe α ) -Π C (u)(x)| ≤ lim ε→0 1 ε |u(x + εe α ) -u(x)| = |∂ α u(x)|
(where e α = (δ α,β ) d β=1 is the canonical basis of R d ). From which our claims follows, as

|∇Π C (u)(x)| 2 = d α=1 |∂ α Π C (u)(x)| 2 ≤ |∇u(x)| 2 .
In particular, the projection decreases the Euclidean total variation.

As seen above, we can cover Ω by finitely many closed sets A associated each to a particular partition of {1, . . . , k -1} into l clusters {k n , . . . , k n+1 }, with the same notation as above. Letting û be defined by ûi

= 1 k n+1 -k n k n+1 -1 j=kn u j whenever k n ≤ i < k n+1 , the set A is simply {x ∈ Ω : Π C (u) = û}. This implies ∇Π C (u) = ∇û a.e. in A, hence A Φ(∇û(x)) dx = A Φ(∇Π C (u)(x)) dx . (63)
Hence to prove the lemma, it is enough to show that Φ(∇û(x)) ≤ Φ(∇u(x))

a.e. in A. But in Ω we have:

∇û i = 1 k n+1 -k n k n+1 -1 j=kn ∇u j for k n ≤ i < k n+1 . Now, if (ξ i ) k-1 i=1 ∈ L, we have for any x ∈ Ω k-1 i=1 ξ i • ∇û i (x) = l n=1 k n+1 -1 i=kn k n+1 -1 j=kn ξ j k n+1 -k n • ∇u i (x) = k-1 i=1 ξi • ∇u i (x) (65)
where for each i with

k n ≤ i < k n+1 , we have let ξi = 1 k n+1 -k n   k n+1 -1 j=kn ξ j   .
If we show that ( ξi ) k-1 i=1 ∈ L, it will follow from (65) that i ξ i ∇û i (x) ≤ Φ(∇u(x)), from which (64) will follow.

We let ξ 0 i = ξ i for each i, and for n = 1, . . . , l we let ξ n i = ξi if i < k n+1 , and

ξ n i = ξ n-1 i = ξ i if i ≥ k n+1 .
In other words, ξ n is obtained from ξ n-1 by averaging all components between k n and k n+1 -1, and leaving the other unchanged. Let us show by induction that (ξ

n i ) k-1 i=1 ∈ L, for each n ≤ l: since ( ξi ) k-1 i=1 = (ξ l i ) k-1 i=1
, the thesis will follow.

6.6 Proof of proposition 6.1

Proof. Let u solve (60). Let u ′ = Π C 0 u. By Lemma 6.4, we have J(u ′ ) ≤ J(u). Let us now consider the other terms in (60), that is,

Ω g 1 (x)|1 -u 1 (x)| + k-1 i=2 g i (x)|u i-1 (x) -u i (x)| + g k (x)|u k-1 (x)| dx
We see that when we replace u with, first, u ′′ = Π C u in this expression, some terms of the form |u i-1 (x)u i (x)| will vanish (hence decreasing the energy, since g i > 0 a.e.), while, thanks to (62), if

|u ′′ i-1 (x) -u ′′ i (x)| = 0, we have u i-1 (x) ≥ u ′′ i-1 (x) > u ′′ i (x) ≥ u i (x) so that |u ′′ i-1 (x) -u ′′ i (x)| ≤ |u i-1 (x) -u i (x)|.
On the other hand, the effect of projecting then u ′′ onto [0, 1] k-1 , to find u ′ , does not alter most terms in the sum from i = 1 to k -1, but strictly reduces the first and last term (|1u 1 (x)| and |u k-1 (x)|) if u ′′ was not already in [0, 1] k-1 , while it might also reduce some other terms of the form |u i-1 (x)u i (x)| whenever u i-1 (x) > u i (x). Hence, the energy is strictly reduced by replacing u with u ′ , unless u = u ′ , that is, u ∈ C 0 . Remark 6.5. The same proof would show that an unconstrained minimizer of H with prescribed Dirichlet conditions in L 1 (∂Ω; C 0 ) will also be in C 0 a.e., hence be a minimizer of the constrained problem of minimizing over C 0 . Corollary 6.6. Assume g i > 0 a.e., for i = 1, . . . , k, and let v solve

min i v i =1 a.e. J(v) + k i=1 Ω g i (x)|v i (x)| dx .
Then v solves [START_REF] Ishikawa | Exact optimization for Markov random fields with convex prio rs[END_REF].

It simply follows from Proposition 6.1 after the changes of variable [START_REF] Bae | Global minimization for continuous multiphase partitioning problems using a dual approach[END_REF] and [START_REF] Bellettini | Convex approximations of functionals with curvature[END_REF]. We also have the following variant which regards the Dirichlet problem (33): Corollary 6.7. Let v ∈ L 2 (Ω; R k ) minimize the energy in (33) without the constraint v ∈ S. Then v ∈ S, and is a solution of (33).

Proof. The fact that Ψ(p) = +∞ if

i p i = 0 yields that any v of finite energy is in BV (Ω; R k ), with i v i = 1 a.e. in Ω (which is inherited from the boundary condition v 0 and the fact i v 0 i = 1 on ∂Ω). Then, using the changes of variable ( 7) and ( 9), together with Remark 6.5, we find that v ∈ S.

Numerical analysis and results

We do not detail our implementation here: indeed, it follows the same lines as in Section 4. The discrete approximation is of the same type and the Γ-convergence result in Proposition 4.3 will still holds in case Ψ is replaced with Φ and J with H. Then, we use the same type of Arrow-Hurwicz' iteration, as described in [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF]. To avoid the reprojection onto C 0 (wich is not so complicated, see Appendix B, we introduce also a dual variable for the terms |u i-1u i | which is updated exactly like the other dual variable.

The slowest part in this method is obviously the computation of the projection onto L of the dual variable (Ξ 1 , . . . , Ξ k-1 ) ∈ R (k-1)×d . This motivates again the use of a GPU. This step is done using Dykstra's algorithm as described in Appendix A. The set L is described as the intersection of k(k + 1)/2 "simple" convex sets and an iterate of Dykstra's algorithm must therefore involve as many projections and dual variables. Also, the primal-dual gap, which has an expression similar to [START_REF] Popov | A modification of the arrow-hurwicz method for search of saddl e points[END_REF], can be computed as before only in the 2 or 3-labels cases (since evaluating Φ, like Ψ, requires an iterative algorithm if k ≥ 4), so that we usually choose an arbitrary stopping criterion for our iterations (like the variation between two successive iterates). We show an example of stereo reconstruction, as described in [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF]. We have implemented with k = 17 labels a disparity reconstruction algorithm, where the data term is made of absolute differences of the color channels between the left and right image, with a linear interpolation as described in [START_REF] Birchfield | A pixel dissimilarity measure that is insensitive to image sampling[END_REF]. We show here the reconstruction with three different interaction potentials: the total variation as in [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF], a "Potts" energy (the length of the discontinuity), and the "truncated TV". Results are shown on Figure 15. The best, and more precise results are obtained with the truncated TV. More precisely, we can count the number of "wrong" pixels where the error in disparity is more than one: in this case, we have found 2.55% of wrong pixels in the result with the TV interaction, 4.22% with the Potts model, and 1.88% with the truncated TV.

A Algorithm for projecting onto K Let us explain here how we implemented, in our practical computation, the projection of a vector (ξ i ) k i=1 in R k×d onto the convex K (or with the same technique, of (ξ i ) k-1 i=1 in R (k-1)×d onto the convex L). These two sets are described as a finite intersection of simple sets, on which a projection is easily computed. In this case, a converging algorithm was first proposed by Dykstra in the 1980's [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF].

A.1 Dykstra's algorithm

The idea is essentially as follows: we are given k convex sets in R N , K 1 , . . . , K k , each on which has a simple algorithm for projecting onto it, and you want to compute the projection of a vector x onto k i=1 K i . Letting ψ i (z) = 0 if z ∈ K i , and +∞ else, it means you want to solve

Then, a good idea is to consider the dual problem, which in this case can be written [START_REF] Ekeland | Convex analysis and variational problems[END_REF][START_REF] Rockafellar | Convex analysis[END_REF] min

y∈R N |x -y| 2 2 + k i=1 ψ i * (y) , (67)
moreover, z solves (66) if and only if ȳ = xz solves (67). Here, ( k i=1 ψ i ) * denotes the Legendre-Fenchel conjugate [START_REF] Ekeland | Convex analysis and variational problems[END_REF][START_REF] Rockafellar | Convex analysis[END_REF] of k i=1 ψ i , that is,

k i=1 ψ i * (y) = sup z∈R N y • z - k i=1 ψ i (z)
Now, it is known (and easy to check) that the Legendre-Fenchel conjugate of a sum of convex, lower semicontinuous functions is given by the inf-convolution of the Legendre-Fenchel conjugate of each function. In this case, that means:

k i=1 ψ i * (y) = inf k i=1 ψ * i (y i ) : k i=1 y i = y , (68)
so that another way to write (67) is: inf

(y i ) k i=1 ∈(R N ) k 1 2 x -k i=1 y i 2 + k i=1 ψ * i (y i ) , (69)
The idea behind Dykstra's algorithm is to minimize now (69) alternatively in each variable (y i ), for i = 1, . . . , k, until convergence (this point of view was is found for instance in Gaffke and Mathar [START_REF] Gaffke | A cyclic projection algorithm via duality[END_REF]).

Hence, we start with y 0 i = 0 for each i, and (y n+1 i ) k i=1 is found from (y n i ) k i=1 by letting, for i = 1, . . . , k:

y n+1 i = arg min y∈R N 1 2
x -j<i y n+1 j -j>i y n jy

2 + ψ * i (y) . (70)
Now, let us set for each n ≥ 1, y n = k i=1 y n i , x n = xy n (and x 0 = x). We also set x n 0 = x n-1 , x n i = x -( j≤i y n j + j>i y n-1 j ) (in particular x n k = x n = x n+1 0

). The primal version of (70) is, for n ≥ 1 and i = 1, . . . , k:

x n+1 i = arg min z∈R N 1 2 x n i-1 + y n i -z 2 + ψ i (z) (71)
while y n+1 i is then given by y n+1 i = x n i-1 + y n ix n+1 i . But, as ψ i is defined, it means nothing else than

x n+1 i = Π K i (x n i-1 + y n i ) y n+1 i = x n i-1 + y n i -x n+1 i (72)
This is the main iteration of Dykstra's algorithm, as described in the original work. Moreover, it is shown there that as n → ∞, x n → Π K (x) (strongly, if we replace R N with a Hilbert space).

A.2 The projection onto K

In order to implement the numerical method described in Section 4, we need to project onto the convex set K ⊂ R k×d defined in [START_REF] Giusti | Minimal surfaces and functions of bounded variation[END_REF] or [START_REF] Tamanini | Optimal segmentation of unbounded functions[END_REF]. which may be written as the intersection of convex sets

1≤i 1 <i 2 ≤k K i 1 ,i 2 , with K i 1 ,i 2 = (q i ) k i=1 ∈ R k×d : |q i 2 -q i 1 | ≤ σ i 1 ,i 2 .
Given q = (q i ) k i=1 ∈ R k×d , the projection onto K i 1 ,i 2 is given by

(Π K i 1 ,i 2 (q)) i =        q i if i = i 1 , i 2 ;
q i 1 + 1 2 q i 1 ,i 2 if i = i 1 q i 2 -1 2 q i 1 ,i 2 if i = i 2 where q i 1 ,i 2 = (|q i 2q i 1 |σ i 1 ,i 2 ) + q i 2q i 1 |q i 2q i 1 | .

Dykstra's algorithm is now simple to implement: we are given a vector q we want to project onto K. We fist let q i 1 ,i 2 = 0, and recursively update as follows:

1. For 1 ≤ i 1 < i 2 ≤ k:

(a) Compute q = q i 2q i 1 + q i 1 ,i 2 ;

(b) Let then q = (|q|σ i 1 ,i 2 ) + (q/|q|);

(c) Assign to q i 1 the value of q i 1 + 1 2 ( qq i 1 ,i 2 ) and to q i 2 the value q i 2 -1 2 ( qq i 1 ,i 2 ); (d) Let q i 1 ,i 2 take the value q.

2. Test how much q = (q i ) k i=1 has changed during the loop: if change is less than some tolerance, end. Else: go back to step 1.

The projection onto L is slightly more complicated but the idea is the same, we skip the details.

A.3 Accelerating the projection onto K

The number of constraints in K grow quadratically in k, whereas for the weak relaxations, the constraints stay linear in k. For very specific functionals such as the total cyclic variation [START_REF] Strekalovskiy | Total variation for cyclic structures: Convex relaxation and efficient minimization[END_REF], the quadratic number of constraints can actually be expressed in an equivalent form which only requires linearly many constraints. Yet in the general case, Dykstra's algorithm becomes quite impractical for a larger number of labels (e.g. k > 5). However, we observed that only a very small fraction (usually below 1%) of the constraints are violated during the iterates of the primal-dual algorithm. We show that this can be used to significantly speeding up the projection onto K, in case the surface tensions satisfy a triangular condition a bit stronger than the natural condition [START_REF] Strekalovskiy | Total variation for cyclic structures: Convex relaxation and efficient minimization[END_REF], see (73) below.

In this Section, we show that indeed only those constraints which are actually violated need to be considered for computing the projection onto K. This leads to a dramatic reduction of the computing time.

Given the surface tensions (σ i,j ) (with σ i,j = σ j,i ), the sets K i,j (1 ≤ i < j ≤ k) are defined as in the previous Section by K i,j = p ∈ R d×k : |p ip j | ≤ σ i,j and we have K = i<j K i,j .

Let us analyse the way the boundaries of the sets K i,j intersect each other. First, the normal at p ∈ ∂K i,j (so that |p ip j | = σ i,j ) is ν i,j (p) = 1 √ 2σ i,j (p ip j ) ⊗ (e ie j )

where here, for p ∈ R d and x ∈ R k , p ⊗ x ∈ R d×k is (p n x i ) i=1,...,k n=1,...,d , and (e i ) i=1,...,k is the canonical basis of R k .

In particular, we see that if i < j, l < m are all different indices, then ∂K i,j and ∂K l,m intersect perpendicularly, as ν i,j (p), ν l,m (p) = 0. (Throughout this section, we use •, • to denote the full scalar product in R d×k , in particular, p ⊗ x, q ⊗ y = (p • q)(x • y).)

Figure 1 :

 1 Figure 1: Color input image and segmentation obtained by minimizing the piecewise constant Mumford-Shah functional with 10 regions.

Figure 2 :

 2 Figure 2: Various popular interaction potentials.

  the simplex) a.e. in Ω. Then, the result follows from quite standard convex analysis, in the Hilbert space L 2 (Ω; R k ), and the characterization of the subgradient of the (anisotropic) total variation with Dirichlet boundary conditions, see[START_REF] Moll | The anisotropic total variation flow[END_REF] Prop. 3] and [4, Thm. 2].

Figure 3 :

 3 Figure 3: The "calibration" for the triple point.

Figure 5 :

 5 Figure 5: Reconstruction of a triple point.

Figure 6 :

 6 Figure 6: Example of a nonbinary solution.

Figure 7 :

 7 Figure 7: Completion of four regions.

Figure 8 :

 8 Figure 8: Completion of four regions: in case of non uniqueness, the method may find a combination of the solutions.

Figure 10 :

 10 Figure 10: Minimal surface problem with three phases in 3D. The red and the gree phases are represented as solid volumes and the blue phase is represented by a slice cutting through the volume. Note that the phases meet at non-trival (curved) surfaces.

Figure 11 :

 11 Figure 11: Minimal surface problem with six phases in 3D. The left image shows a volume rendering of the result computed by our approach. The right image shows a natural minimal surface generated by a soap film in a cubic frame. Note the little square in the middle of the cube which is spanned by all six phases.

Figure 12 :Figure 13 :

 1213 Figure 12: A calibration for the triple point and problem (49) does not exist

Figure 13

 13 Figure13illustrates the negative result of Proposition 5.1: the relaxation J is the only one which reconstructs correctly the triple point.

6

  Extension to more general partition problems 6.1 The convex, local envelope of the general partition problem Now, we explain how our setting is generalized to more general interaction energies. As before, d ≥ 1 is the space dimension, k ≥ 2 the number of labels, Ω a bounded open subset of R d . By a general partition problem, we mean the problem of finding a partition (E i ) k i=1 which minimizes:

  , as in general it follows that σ(a) + σ(b) ≥ σ(a + b) for any a, b > 0. Indeed, a = (a/(a + b))(a + b) + (1a/(a + b)) × 0, b = (b/(a + b))(a + b) + (1b/(a + b)) × 0, so that (letting σ(0) = inf t>0 σ(t) = lim t→0 σ(t) ≥ 0) σ(a) ≥ a a + b σ(a + b) + b a + b σ(0) and σ(b) ≥ b a + b σ(a + b) + a a + b σ(0) , from which follows σ(a) + σ(b) ≥ σ(a + b) + σ(0) ≥ σ(a + b)

Figure 14 :Figure 15 :

 1415 Figure 14: A stereo pair

2 Figure 16 :

 216 Figure 16: Construction of the optimal point x

  k,j=1,...,d . The set J v is shown to be a (d -1)-dimensional set which is rectifiable in the sense of Federer (that is, can be covered by countably many C 1 hypersurfaces, up to a set which is negligible for the (d -1)-dimensional Hausdorff measure, H d-1 ). In particular, ν v coincides H d-1 -a.e. with a normal to J v . Eventually, the cantor part Cv is the part of the measure which, essentially, has dimension between d and d -1: it is singular with respect to the Lebesgue measure, but also satisfies |Cu|(A) = 0 for any set A with H d-1 (A) < +∞.Given a convex, continuous function Ψ : Ω × R k×d → [0, +∞), one can define the integral

Ω Ψ(x, Dv) as follows: one introduces the recession function

It is straightforward to extend the problem to nonsymmetric weight σi,j = σj,i, and more general weights depending also on the direction of the jump, but to simplify we stick to the isotropic case. For a recent generalization allowing to impose ordering constraints and direction-dependent penalties we refer to[START_REF] Strekalovskiy | Generalized ordering constraints for multilabel optimization[END_REF].

Acknowledgements

We thank François Alouges, Bastian Goldlücke, Robert V. Kohn, Olivier Pantz and Margherita Solci for fruitful discussions, as well as the Hausdorff Institute for Mathematics in Bonn for financial support. We also thank Bernhard Kainz for helping with the 3D renderings, and Frank Pacard for providing the cubic frame and soap for Fig. 11. The first author acknowledges the hospitality of the University of Bonn, and in particular the Computer Science Department and the Hausdorff Institute for Mathematics. He was partially supported by the Agence Nationale de la Recherche, project ANR-06-BLAN0082. The second author acknowledges support through an ERC starting grant. The third author acknowledges support from the Austrian Science Fund (FWF) under the grant P22492-N23.

Remark 6.2. As already observed, there is no loss of generality in assuming that the g i are positive, as soon as they are bounded from below (even, in fact, by a given L 1 function which does not depend on the label i), as the problem remains unchanged if the same integrable function is added to all the functions g i .

The proof of Proposition 6.1 relies on some truncation properties of the functional H and will be a consequence of the Lemmas 6.3 and 6.4 which are proved in the following sections 6.4 and 6.5: we postpone it to section 6.6.

Ordering of the vector u

Let us denote by

we denote by Π C (u) the orthogonal projection of u onto C. The projection u ′ = Π C (u) is a vector of the following form: there exist 1 = k 1 ≤ k 2 ≤ • • • ≤ k l+1 = k such that for each n = 1, . . . , l and each i with k n ≤ i < k n+1 ,

Indeed, the sets {k n , . . . , k n+1 -1}, when containing more than one index, are just the maximal clusters of indices whose associated coefficients become equal in the projection, and u ′ is then simply the projection of u onto {z ∈ R k-1 :

In addition, we claim that the projection u ′ must satisfy, for each n = 1, . . . , l,

Indeed, if for instance u kn > u ′ kn for some n, we define a vector u ′′ ∈ C by

Observe in particular that (62) yields that if

An recursive algorithm for computing the projection Π C (u) can be deduced from this remark, see Appendix B. Now, let us show that the projection onto C decreases the energy H:

Here, Π C (u) denotes the function x → Π C (u(x)).

Proof. First, using Theorem 2.1 and the convergence ( 16), together with the lower-semicontinuity of H, it is enough to show the lemma for a smooth u ∈ C ∞ (Ω; R k-1 ). In particular, Π C (u) is (Lipschitz) continuous: we even have (since

This is true for n = 0, by assumption. Assume n ≥ 1 and (ξ n-1 i ) k-1 i=1 ∈ L. The only non-obvious conditions to check to prove that (ξ n i ) k-1 i=1 ∈ L are of the kind

In the first case, for instance, simple algebra shows that this sum is the convex combination

hence its norm is less than, using the concavity of σ:

and the thesis follows. The second case is treated in the same way.

This achieves the proof of the lemma.

Truncation of the coordinates

Now, we define

Then, we have in addition:

Here again,

e, moreover, as before, we may assume that u is smooth. Then, the same arguments as in the proof of the previous lemma show that the inequality is true provided we can show that for any i 1 , i 2 with 1 ≤ i 1 ≤ i 2 ≤ k -1, and any p ∈ R (k-1)×d , if p is defined by pi = p i when i 1 ≤ i ≤ i 2 and pi = 0 else, Φ(p) ≤ Φ(p). But this is an obvious consequence of the fact that if (ξ i ) k-1 i=1 ∈ L, also ( ξi ) k-1 i=1 , defined by ξi = ξ i whenever i 1 ≤ i ≤ i 2 and ξi = 0 else, belongs to L, which follows from the assumption that σ is nondecreasing.

If now p ∈ K, p ∈ ∂K i,j ∩ ∂K l,j , we find in the same way

In all cases, we see that if the surface tensions satisfy

for all triple of distinct indices i, j, l, then we have that the intersection of any two surfaces ∂K i,j and ∂K l,m satisfies ν i,j (p), ν l,m (p) ≥ 0.

(74)

Now, we recall that for any convex set C, its normal cone at x ∈ C is the set N C (x) of vectors ν such that yx, ν ≤ 0 ∀y ∈ C .

Moreover, it coincides with the subgradient ∂δ C (x) of the characteristic function

In case C, C ′ are two convex set and C ∩ C ′ has nonempty interior, we deduce that [50, Cor. 23.8.1]

In the present case, we deduce that for any I ⊂ {(i, j) : i < j}, the normal cone to K I = (i,j)∈I K i,j at p is made of all the vectors λ i,j ν i,j (p), where (λ i,j ) are non-negative numbers and (i, j) runs on all the pairs (i, j) ∈ I such that |p ip j | = σ i,j . In particular, we deduce that if p ∈ ∂K I ∩ ∂K I ′ , then

The following result follows:

Proposition A.1. Assume (73) holds for all i, j, l = 1, . . . , k, and let p ∈ R d×k . Let

Proof. By continuity it is enough to prove the proposition when p is in the interior of K p , which we will assume from now on. Let q = Π K (p) = Π L p ∩K p (p). Then, q ∈ ∂(L p ∩ K p ) which can be split into three parts, ∂L p ∩ int(K p ), int(L p ) ∩ ∂K p , and ∂L p ∩ ∂K p . If q ∈ ∂L p ∩ int(K p ), then the first-order condition of optimality for q will easily show that q = Π L p (p), which is our thesis.

Hence it is enough to rule out the two other possibilities, which are q ∈ int(L p ) ∩ ∂K p and q ∈ ∂L p ∩ ∂K p . The first one is obviously impossible, as in that case, for ε > 0 small, q + ε(pq) ∈ K = L p ∩ K p is closer to p than q, a contradiction. Hence it remains to show that q ∈ ∂L p ∩ ∂K p .

Assume q ∈ ∂L p ∩ ∂K p . First of all, one can find a unit vector ν K p ∈ N K p (q) such that for t > 0 small, qtν K p is in the interior of a cone of vertex q contained in K p (as it will lie in the interior of the dual cone to N K p (q). A possible choice, here, is to take ν K p equal to the renormalized sum of all the normals ν i,j (p) ((i, j) ∈ I p , |q iq j | = σ i,j ) extremal to the cone N K p (q), which, thanks to (75), will satisfy ν K p , ν i,j (p) > 0 for all i, j.

If for t small, we also have qtν K p ∈ L p then we're done, indeed, we deduce that pq, -tν K p ≤ 0, but since p ∈ K p and ν K p ∈ N K p (q), one also has pq, ν K p ≤ 0. Hence pq, ν K p = 0 which is possible only if p ∈ ∂K p , a contradiction (we have assumed p in the interior).

If qtν K p ∈ L p for any t > 0 small, let us show that it is not too far anyway, that is:

Let z t ∈ L p be such that qtν K pz t = dist(qtν K p , L p ). Clearly, this distance is less than t (since q ∈ L p ). Assume by contradiction that for a sequence (t n ) n ↓ 0, it is larger than δt n for some δ > 0. Possibly extracting a further subsequence we may assume that (qz tn )/t nν K p converges to some vector ν of norm between δ and 1, moreover, clearly, ν ∈ N L p (q) and in particular ν, ν K p ≥ 0. As z tn is the projection onto L p of qt n ν K p and q ∈ L p we also have

hence in the limit ν, ν K p ≤ -δ 2 , a contradiction.

Hence, (76) holds and qtν K pz t = o(t), in particular, (qz t )/t → ν K p as t → 0. We also deduce that z t ∈ K p if t is small enough, hence z t ∈ K p ∩ L p = K, so that pq, z tq ≤ 0 Dividing by t and sending t → 0, we deduce that pq, ν K p ≥ 0, but since p ∈ int(K p ), again, and ν K p ∈ N K p (q), it follows that pq, ν K p = 0 and p ∈ ∂K p , a contradiction.

B Algorithm for projecting onto C

Let u ∈ R k-1 : we give here a simple algorithm for computing the projection Π C u of u onto

Of course, as in the previous section, we could still use Dykstra's algorithm: however, this will usually converge in infinitely many iterations (although the error decreases very fast), whereas the alternative solution we propose here needs at most k -2 iterates to produce the exact solution.

Let us first modify slightly the problem, and consider the minimization:

where ω i are positive weight. Then, the same analysis as in section 6.4 holds true, and in particular (61): we hence know that if for some index i ∈ {1, . . . , k -2}, u i ≤ u i+1 , the solution z of (77) satisfies z i = z i+1 (while if there is no such index, u ∈ C and z = u is the solution).

In this case, (77) is obviously equivalent to the minimization of

over all z ∈ C with z i = z i+1 , which (since the last line in (78) does not depend on z) is a new problem of the form (77), but now in dimension (k -2), and with the coordinates u i , u i+1 replaced with their average with respective weights ω i and ω i+1 , and their common weight in the new distance replaced with the sum of the weights ω i + ω i+1 . Hence, a straightforward recursive algorithm for computing Π C u is as follows:

1. Let first z = u, and define k -1 "clusters" C i = {i}, i = 1, . . . , k -1, of only one element.

2. Identify two coordinates z i , z i+1 with z i < z i+1 . If there are no such coordinates, z ∈ C: the procedure ends and z = Π C u.

3. Replace the clusters C j , j ∈ C i ∪ C i+1 , with the union C i ∪ C i+1 . Replace then z i and z i+1 with the average j∈C i z j /(♯ C i ).

4. go back to 2.

C Evaluation of Ψ for three labels

Let us quickly mention how, in case k = 3, the value of Ψ(q) defined in [START_REF] Hammer | Roof duality, complementation and persistency in quadratic 0-1 optimization[END_REF] can be computed in order to estimate the primal-dual gap [START_REF] Popov | A modification of the arrow-hurwicz method for search of saddl e points[END_REF].

We use the formula (54) to write, given q = (q 1 , q 2 , q 3 ) with q 1 + q 2 + q 3 = 0 (otherwise Ψ(q) = +∞), Ψ(q 1 , q 2 , q 3 ) = Φ(p 1 , p 2 )

with p 1 = -q 1 , p 2 = -q 1q 2 = q 3 . Using (55) (for weights σ i,j = 1), we have that (here p = (p 1 , p 2 ))

Standard duality shows that this can be rewritten

that is, q is the point whose total distance to 0, p 1 and p 2 in R 2 is minimal. Hence, the problem boils down to finding a simple way to compute, given a, b, c ∈ R 2 , the point x ∈ R 2 which minimizes |x -a| + |x -b| + |x -c|.

If a, b, c are aligned (in our case, if det(p 1 , p 2 ) = 0), then x is optimal when it is equal to a point which lies in between the two other (or is equal to one of the other). For instance, if b ∈ [a, c], x = b is optimal and the solution is |a -b| + |b -c|.

If the points are not on the same line, then they form a non-degenerate triangle. It is shown, then, that if x is optimal, it must lie inside the triangle formed by the three points. Then, from the optimality conditions, one checks that either x is one of the points a, b, c (say, for instance, a), in which case the angle bac must be greater than, or equal to 120 • , or x ∈ {a, b, c}, in which case the angles abc, bac, acb are all three smaller than 120 • . In this last non-degenerate case, the three angles axb, bxc, cxa must be exactly 120 • , so that x lies at the intersection of three circles, which are the circles made of the points which "see" two vertices of the triangle under an angle of 120 • . This is very easy to compute. The centers of these circles are A, B, C: A is at distance |bc|/(2 √ 3) of the segment [b, c] and its projection on this segment is the middle point (b + c)/2 (and a and A must be separated by (bc)). B is the same with respect to a and c, and C with respect to a and b. Then, the circles of center A and radius |bc|/ √ 3 and of center B and radius |ac|/ √ 3 have the two intersection points c and x: x is then the symmetric of c with respect to the line (AB). See Fig. 16.