Camille Couprie
email: c.couprie@esiee.fr

Leo Grady
email: leo.grady@siemens.com

Hugues Talbot
email: h.talbot@esiee.fr

Laurent Najman
email: l.najman@esiee.fr

Combinatorial Continuous Maximal Flows

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1. Introduction. Optimization methods have been used to address a wide variety of problems in computer vision. The early optimization approaches were formulated in terms of active contours and surfaces [START_REF] Kass | Snakes: Active contour models[END_REF] and then later level sets [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF]. These formulations were used to optimize energies of varied sophistication (e.g., using regional, texture, motion or contour terms [START_REF] Paragios | Geodesic Active Regions and Level Set Methods : Contributions and Applications in Artificial Vision[END_REF]) but generally converged to a local minimum, generating results that were sensitive to initial conditions and noise levels. Consequently, more recent focus has been on energy formulations (and optimization algorithms) for which a global optimum can be found.

Such energy formulations typically include a term which minimizes the boundary length (or surface area in 3D) of a region or the total variation of a scalar field in addition to a data term and/or hard constraints. In this paper, we focus on image segmentation as the example application on which to base our exposition. Indeed, segmentation has played a prominent (and early) role in the development of the use of global optimization methods in computer vision, and often forms the basis of many other applications [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF][START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF][START_REF] Darbon | Image restoration with discrete constrained total variation part ii: Levelable functions, convex priors and non-convex cases[END_REF].

Graph-based approaches. The max-flow/min-cut problem on a graph is a classical problem in graph theory, for which the earliest solution algorithm goes back to Ford and Fulkerson [START_REF] Ford | Maximal flow through a network[END_REF]. Initial methods for global optimization of the boundary length of a region formulated the energy on a graph and relied on max-flow/min-cut methods for solution [START_REF] Boykov | Interactive graph cuts for optimal boundary and region segmentation of objects in n-d images[END_REF]9]. It was soon realized that these methods introduced a metrication error. Metrication errors are clearly visible when gradient information is weak or missing, for example in the case of medical image segmentation or in some materials science applications like electron tomography, where weak edges are unavoidable features of the imaging modality. Metrication errors are also far more obvious in 3D than in 2D and increasingly so, as the dimensionality increases (see for example the 3D-lungs example in [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF]). Furthermore, metrication artifacts are even more present in other applications such as surface reconstruction and rendering, where the artifacts are a lot worse than in image segmentation. Various solutions for metrication errors were proposed. One solution involved the use of a highly connected lattice with a specialty edge weighting [START_REF] Boykov | Computing geodesics and minimal surfaces via graph cuts[END_REF], but the large number of graph edges required to implement this solution could cause memory concerns when implemented for large 2D or 3D images [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF].

Continuous approaches. To avoid the metrication artefacts without increasing memory usage, one trend in recent years has been to pursue spatially continuous formulations of the boundary length and the related problem of total variation [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF][START_REF] Nikolova | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF]. Historically, a continuous max-flow (and dual min-cut problem) was formulated by Iri [START_REF] Iri | Theory of flows in continua as approximation to flows in networks[END_REF] and then Strang [START_REF] Strang | Maximum flows through a domain[END_REF]. Strang's continuous formulation provided an example of a continuization (as opposed to discretization) of a classically discrete problem, but was not associated to any algorithm. Work by Appleton and Talbot [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF] provided the first PDE-based algorithm to find Strang's continuous max-flows and therefore optimal min-cuts. This same algorithm was also derived by Unger et al. from the standpoint of minimizing continuous total variation [START_REF] Unger | TVSeg -Interactive total variation based image segmentation[END_REF]. Adapted to image segmentation, this algorithm is shown to be equivalent [START_REF] Unger | Interactive globally optimal image segmentation[END_REF] to the Appleton-Talbot algorithm and has been demonstrated to be fast when implemented on massively parallel architectures. Unfortunately, this algorithm has no stopping criteria and has not been proved to converge. Works by Pock et al [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF], Zach et al [19] and Chambolle et al [START_REF] Chambolle | A convex approach for computing minimal partitions[END_REF] present different algorithms for optimizing comparable energies for solving multilabel problems, but again, those algorithms are not proved to converge. Some other works have been presenting provably converging algorithms, but with relatively slow convergence speed. For example, Pock and Chambolle introduce in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] a general saddle point algorithm that may be used in various applications, but needs half an hour to segment a 350 × 450 image on a CPU and still more than a dozen seconds on a GPU (Details are given in the experiments section). Links and differences with total variation minimization. G. Strang [START_REF] Strang | Maximum flows and minimum cuts in the plane[END_REF] has shown the continuous max-flow problem for the l 2 norm to be the dual of the total variation (TV) minimization problem under some assumptions. The problem of total variations, introduced successively in computer vision by Shulman and Hervé [START_REF] Shulman | Regularization of discontinuous flow fields[END_REF] and later Rudin, Osher and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] as a regularizing criterion for optical flow computation and then image denoising, has been shown to be quite efficient for smoothing images without affecting contours. Moreover, a major advantage of TV is that it is a convex problem. Thus, a straightforward algorithm such as gradient descent may be applied to find a minimum solution. However, there is a need for faster methods, and significant progress have been achieved recently with primaldual approaches [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF], Nesterov's algorithm [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] and Split-Bregman/Douglas-Rachford methods [START_REF] Goldstein | The split bregman method for l1-regularized problems[END_REF][START_REF] Zhang | Bregmanized nonlocal regularization for deconvolution and sparse reconstruction[END_REF]. Most methods minimizing TV focus on image filtering as an application, and even if those methods are remarkably fast in denoising applications, segmentation problems require a lot more iterations for those algorithms to converge. As stated previously, continuous max-flow is dual with total variation in the continuous setting under restrictive regularity conditions [START_REF] Nozawa | Max-flow min-cut theorem in an anisotropic network[END_REF]. In fact, Nozawa [START_REF] Nozawa | Examples of max-flow and min-cut problems with duality gaps in continuous networks[END_REF] showed that there is a duality gap between weighted TV and weighted max-flow under some conditions in the continuous domain. Thus, it is important to note that weighted TV problems are not equivalent to the weighted maximum flow problem studied in this paper. Furthermore, many works assume that the continuous-domain duality holds algorithmically, but we show later in this paper that at least in the combinatorial case this is not true.

The previous works for solving the max-flow problem illustrate two difficulties with continuous-based formulation, that are (1) the discretization step, which is necessary for deriving algorithms, but may break continuous-domain properties; and (2) the convergence, both of the underlying continuous formulation, and the associated algorithm, which itself depend on the discretization. It is very well known that even moderately complex systems of PDEs may not converge, and even existence of solutions is sometimes not obvious [START_REF] Fefferman | Existence and smoothness of the Navier-Stokes equation[END_REF]. Even when existence and convergence proofs both exist, sometimes algorithmic convergence may be slow in practice. For these reasons, combinatorial approaches to the maximal flow and related problems are beginning to emerge.

Combinatorial approaches. Discrete calculus [START_REF] Hirani | Discrete exterior calculus[END_REF][START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF] has been used in recent years to produce a combinatorial reformulation of continuous problems onto a graph in such a manner that the solution behaves analogously to the continuous formulation (e.g., [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF][START_REF] Grady | The piecewise smooth Mumford-Shah functional on an arbitrary graph[END_REF]).

In particular, Lippert presented in [START_REF] Lippert | Discrete approximations to continuum optimal flow problems[END_REF] a combinatorial formulation of an isotropic continuous max-flow problem on a planar lattice, making it possible to obtain a provably optimal solution. However, in Lippert's work, parameterization of the capacity constraint is tightly coupled to the 4-connected squared grid, and the generalization to higher dimensions seems non-trivial. Furthermore it involves the multiplication of the number of capacity constraints by the degree of each node, thus increasing the dimension of the problem. This formulation did not lead to a fast algorithm. The author compared several general solvers, quoting an hour as their best time for computing a solution on a 300 × 300 lattice.

Motivation and contributions.

In this paper, we pursue a combinatorial reformulation of the continuous max-flow problem initially formulated by Strang, which we term combinatorial continuous maximum flow (CCMF). Viewing our contribution differently, we adopt a discretization of continuous max-flow as the primary problem of interest, and then we apply fast combinatorial optimization techniques to solve the discretized version.

Our reformulation of the continuous max-flow problem produces a (divergencefree) flow problem defined on an arbitrary graph. Strikingly, CCMF are not equivalent to the discretization of continuous max-flows produced by Appleton and Talbot or that of Lippert (if for no other reason than the fact that CCMF is defined on an arbitrary graph). Moreover, CCMF is not equivalent to the classical max-flow on a graph. In particular, we will see that the difference lies in the fact that capacity constraints in classical max-flow restrict flows along graph edges while the CCMF capacity constraints restrict the total flow passing through the nodes.

The CCMF problem is convex. We deduce an expression of the dual problem, which allows us to employ a primal-dual interior point method to solve it, such as interior point methods have been used for graph cuts [START_REF] Bhusnurmath | Graph cuts via ℓ 1 norm minimization[END_REF] and second order cone problems in general [START_REF] Goldfarb | Second-order cone programming based methods for total variation image restoration[END_REF]. The CCMF problem has several desirable properties, including:

1. the solution to CCMF on a 4-connected lattice avoids the metrication errors. Therefore, the gridding error problems may be solved without the additional memory required to process classical max-flow on a highly-connected graph; 2. in contrast to continuous max-flow algorithms of Appleton-Talbot (AT-CMF) and equivalent, the solution to the CCMF problem can be computed with guaranteed convergence in practical time; 3. the CCMF problem is formulated on an arbitrary graph, which enables the algorithm to incorporate nonlocal edges [START_REF] Buades | A non-local algorithm for image denoising[END_REF][START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF][START_REF] Grady | Faster graph-theoretic image processing via small-world and quadtree topologies[END_REF] or to apply it to arbitrary clustering problems defined on a graph; and 4. the algorithm for solving the CCMF problem is fast, easy to implement, compatible with multi-resolution grids and is straightforward to parallelize. Our computation of the CCMF dual further reveals that duality between total variation minimization and maximum flow does not hold for CCMF and combinatorial total variation (CTV). The comparison between those two combinatorial problems is motivated by several interests:

1. for clarification: in the continuous domain, the duality between TV and MF holds under some regularity constraints. In the discrete anisotropic domain (i.e. Graph Cuts), this duality always holds. However, In the isotropic weighted discrete case, i.e. in the CTV, AT-CMF or CCMF cases, we are not aware of any discretization such that the duality holds. It is not obvious to realize that the duality does not hold in the discrete setting even though it may in the continuous case. We present clearly the differences between the two problems, theoretically, and in term of results; 2. to expose links between CCMF and CTV: We prove that the weak duality holds between the two problems; and 3. for efficiency: in image segmentation, the fastest known algorithms to optimize CTV are significantly slower than CCMF. Therefore, there is reason to believe that CCMF may be used to efficiently optimize energies for which TV has been previously shown to be useful (e.g., image denoising). In the next section, we review the formulation of continuous max-flow, derive the CCMF formulation and its dual and then provide details of the fast and provably convergent algorithm used to solve the new CCMF problem.

2. Method. Our combinatorial reformulation of the continuous maximum flow problem leads to a formulation on a graph which is different from the classical maxflow algorithm. Before proceeding to our formulation, we review the continuous maxflow and the previous usage of this formulation in computer vision.

2.1. The continuous max-flow (CMF) problem. First introduced by Iri [START_REF] Iri | Theory of flows in continua as approximation to flows in networks[END_REF], Strang presents in [START_REF] Strang | Maximum flows through a domain[END_REF] an expression of the continuous maximum flow problem max F st ,

s.t. ∇ • -→ F = 0, || -→ F || ≤ g. (2.1)
Here we denote by F st the total amount of flow going from the source location s to the sink location t, -→ F is the flow, and g is a scalar field representing the local metric distortion. The solution to this problem is the exact solution of the geodesic active contour (GAC) (or surface) formulation [START_REF] Caselles | Geodesic active contours[END_REF]. In order to solve the problem (2.1), the Appleton-Talbot algorithm (AT-CMF) [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF] solves the following partial differential equation system:

∂P ∂τ = -∇ • -→ F , ∂ -→ F ∂τ = -∇P, (2.2)
s. t. || -→ F || ≤ g.
Here P is a potential field similar to the excess value in the Push-Relabel maximum flow algorithm [START_REF] Goldberg | A new approach to the maximum flow problem[END_REF]. AT-CMF is effectively a simple continuous computational fluid dynamics (CFD) simulation with non-linear constraints. It uses a forward finitedifference discretization of the above PDE system subject to a Courant-Friedrich-Levy (CFL) condition, also seen in early level-sets methods. At convergence, the potential function P approximates an indicator function, with 0 values for the background labels, and 1 for the foreground, and the flow -→ F has zero-divergence almost everywhere. However, there is no guarantee of convergence for this algorithm and, in practice, many thousands of iterations can be necessary to achieve a binary P , which can be very slow.

Although Appleton-Talbot is a continuous approach, applying this algorithm to image processing involves a discretization step. The capacity constraint || -→

F || ≤ g is interpreted as max F x(i) 2 + max F y(i) 2 ≤ g 2 i , with F x(i)
the outgoing flow of edges linked to node i along the x axis, and F y(i) the outgoing flow linked to node i along the y axis. We notice that the weights are associated with point locations (pixels), which will correspond later to a node-weighted graph.

In the next section, we use the operators of discrete calculus to reformulate the continuous max-flow problem on a graph and show the surprising result that the graph formulation of the continuous max-flow leads to a different problem from the classical max-flow problem on a graph.

2.2.

A discrete calculus formulation. Before establishing the discrete calculus formulation of the continuous max-flow problem, we specify our notation. A graph consists of a pair G = (V, E) with vertices v ∈ V and edges e ∈ E ⊆ V × V . A transport graph G(V, E) comprises two additional nodes, named a source s and a sink t, and additional edges linking some nodes to the sink and some to the source. Including the source and the sink nodes, the cardinalities of G are given by n = |V | and m = |E|. An edge, e, spanning two vertices, v i and v j , is denoted by e ij . In this paper we deal with weighted graphs that include weights on both the edges and nodes. An edge weight is a value assigned to each edge which may be viewed as a capacity in the context of a maximum flow problem. The weight of an edge, e ij , is denoted by gij . In this work, we assume that gij ∈ R and gij > 0, and use g to denote the vector of R m containing the gij for every edge e ij of G. In addition to edge weights, we may also assign weights to nodes. The weight of node v i is denoted by g i . In this work, we also assume that g i ∈ R and g i > 0. We use g to denote the vector of R n containing the g ij for every edge e ij of G. We define a flow through edge e ij as F ij where F ij ∈ R and use the vector F ∈ R m to denote the flows through all edges in the graph . Each edge is assumed to be oriented, such that a positive flow on edge e ij indicates the direction of flow from v i to v j , while a negative flow indicates the direction of flow from v j to v i .

The incidence matrix is a key operator for defining a combinatorial formulation of the continuous max-flow problem. Specifically, the incidence matrix A ∈ R m×n is known to define the discrete calculus analogue of the gradient, while A T is known to define the discrete calculus analogue of the divergence (see [START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF] and the references therein). The incidence matrix maps functions on nodes (a scalar field) to functions on edges (a vector field) and may be defined as

A eij v k =      +1 if i = k, -1 if j = k, 0 otherwise, (2.3)
for every vertex v k and edge e ij . In our formulation of continuous max-flow, we use the expression |A| to denote the matrix formed by taking the absolute value of each entry individually. Given these definitions, we now produce a discrete (combinatorial) version of the continuous max-flow of (2.1) on a transport graph. As in [START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF][START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF][START_REF] Grady | The piecewise smooth Mumford-Shah functional on an arbitrary graph[END_REF], the continuous vector field indicating flows may be represented by a vector on the edge set, F . Additionally, the combinatorial divergence operator allows us to write the first constraint in (2.1) as A T F = 0. The second constraint in (2.1) involves the comparison of the point-wise norm of a vector field with a scalar field. Therefore, we can follow [START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF][START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF][START_REF] Grady | The piecewise smooth Mumford-Shah functional on an arbitrary graph[END_REF] to define the point-wise ℓ 2 norm of the flow field F as |A T |F 2 . In our notations here, as in the rest of the paper, F 2 = F • F denotes a element-wise product, "•" denoting the Hadamard (element-wise) product between the vectors, and the square root of a vector is also here and in the rest of the paper the vector composed of the square roots of every elements.

Putting these pieces together, we obtain max F st ,

s.t. A T F = 0, |A T |F 2 ≤ g 2 .
(2.4)

Compare this formulation to the classical max-flow problem on a graph, given in our notation as max F st ,

s.t. A T F = 0, |F | ≤ g. (2.5)
By comparing these formulations of the traditional max-flow with our combinatorial formulation of the continuous max-flow, it is apparent that the key difference between the classical formulation and our combinatorial continuous max-flow is in the capacity constraints. In both formulations, the flow is defined through edges, but in the classical case the capacity constraint restricts flow through an edge while the CCMF formulation restricts the amount of flow passing through a node by taking in account the neighboring flow values. This contrast-weighting applied to nodes (pixels) has been a feature of several algorithms with a continuous formulation, including geodesic active contours [START_REF] Caselles | Geodesic active contours[END_REF], CMF [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF], TVSeg [START_REF] Unger | TVSeg -Interactive total variation based image segmentation[END_REF], and [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. On the other hand, the problem defined in (2.4) is defined on an arbitrary graph in which contrast weights (capacities) are almost always defined on edges. The node-weighted capacities fit Strang's formulation of a scalar field of constraints in the continuous max-flow formulation and therefore our graph formulation of Strang's formulation carries over these same capacities. The most important point here is not having expressed the capacity on the nodes of the graph but rather how the flow is enforced to be bounded by the metric in each node. Bounding the norm of neighboring flow values in each node simulates a closer behavior of the continuous setting than if no contribution of neighboring flow was made as in the standard max-flow problem. Figure 2.1 illustrates the relationship of the edge capacity constraints in the classical max-flow problem to the node capacity constraints in the CCMF formulation. The null-divergence constraint is also essential in our formulation since it permits us to obtain constant partitions almost everywhere in the dual solution introduced in the next section, in particular binary ones. The difference between classical max-flow on a graph with the combinatorial continuous max-flow (CCMF) on a graph is that classical max-flow uses edge-weighted capacities while CCMF uses node-weighted capacities. This difference is manifest in the different solutions obtained for both algorithms and the algorithms required to find a solution. Specifically, the solution to the CCMF problem on a lattice does not exhibit metrication bias.

In the context of image segmentation, the vector g varies inversely to the image gradient. We propose to use, as in [START_REF] Unger | TVSeg -Interactive total variation based image segmentation[END_REF],

g = exp(-β||∇I|| 2 ), (2.6)
where I indicates the image intensity. For simplicity, this weighting function is defined for greyscale images, but g could be used to penalize changes in other relevant image quantities, such as color or texture. Before addressing the solution of the CCMF problem, we consider the dual of the CCMF problem and, in particular, the sense in which it represents a minimum cut. Since the cut weights are present on the nodes rather than the edges, we must expect the minimum cut formulation to be different from the classical minimum cut on a graph.

2.3.

The CCMF dual problem. Classically, the max-flow problem is dual to the min-cut problem, allowing a natural geometric interpretation of the objective function. In order to provide the same interpretation, we now consider the dual problem to the CCMF and show that we optimize a node-weighted minimum cut. In the following proposition, we denote the element-wise quotient of two vectors u

= [u 1 , . . . , u k ] and v = [v 1 , . . . , v k ] by u • /v = [u 1 /v 1 , . . . , u k /v k ].
We also denote 1 n a unit vector [1, . . . , 1] of size n, and recall that the square exponent v 2 of a vector v represents the resulting vector of the element-wise mulitiplication v • v.

Proposition 2.1. In a transport graph G with m edges, n nodes, we define a vector c of R m , composed of zeros except for the element corresponding to the source/sink edge which is 1. Let λ and ν be two vectors of R n . The CCMF problem max

F ∈R m c T F, s. t. A T F = 0, (2.7) |A T |F 2 ≤ g 2 .
has for dual

min λ∈R n , ν∈R n λ T g 2 + 1 4 1 n • /(|A|λ) T (c + Aν) 2 , (2.8) s. t. λ ≥ 0,
equivalently written in (2.15), and the optimal solution (F * , λ * , ν * ) verifies

max F c T F = c T F * = 2λ * T g 2 , (2.9)
and the n following equalities

λ * • |A T | (c + Aν * ) • /(|A|λ * ) 2 = 4λ * • g 2 .
(2.10)

Proof. The Lagrangian of (2.7) is

L(F, λ, ν) = (c T + ν T A T )F + λ T |A T |F 2 -λ T g 2 ,
with λ ∈ R n and ν ∈ R n two Lagrange multipliers. We have to find F in the dual function is such as

∇ F L(F, λ, ν) = 0. ∇ F L(F, λ, ν) = c + Aν + 2|A|λ • F = 0 ⇔ F = (c + Aν) • /(-2|A|λ). (2.11)
Substituting F in the Lagrangian function , we obtain

L(λ, ν) = (c T +ν T A T ) (c+Aν)•/(-2|A|λ) +λ T |A T | (c + Aν) • /(-2|A|λ) 2 -λ T g 2 .
The Lagrangian may be simplified as:

L(λ, ν) = -1 n • /(4|A|λ) T (c + Aν) 2 -λ T g 2 .
The dual of (2.7) is also min

λ,ν λ T g 2 + 1 4 1 n • /(|A|λ) T (c + Aν) 2 , (2.12) s. t. λ ≥ 0.
Furthermore, the KKT optimality conditions give

λ * • |A T |F * 2 -g 2 = 0. (2.13)
Using the expression of F * from equation (2.11) , and substituting it in (2.13), we obtain

λ * • |A T | (c + Aν * ) • /(|A|λ * ) 2 = 4λ * • g 2 .
Taking the sum of right hand side and left hand side,

(c + Aν * ) 2 • /(|A|λ * ) = 4λ * T g 2 . (2.14)
Finally, if we substitute (2.14) in the dual expression (2.8), we have max

F c T F = c T F * = 2λ * T g 2 . (a) (b) λ (c) ν (d) Threshold ν at .5 Fig. 2.2.
The dual problem to CCMF is a node-weighted minimum cut in which the variable λ is a weighted indicator vector labeling boundary nodes and the variable ν is a nearly binary vector indicating the source/sink regions. As a result, the contours of ν are slightly blurry. This is due to the equilibrium effect between the two dual variables. In practice, as λ is nonzero only in presence of a contour, ν is binary almost everywhere, except on a very thin line.

The expression of the CCMF dual may be written in summation form as

min λ,ν v i ∈V weighted cut λig 2 i + smoothness term 1 4 e ij ∈E\{s,t} (νi -νj) 2 λi + λj + source/sink enforcement 1 4 (νs -νt -1) 2 λs + λt s. t. λi ≥ 0 ∀i ∈ V.
(2.15)

Interpretation: The optimal value λ * is a weighted indicator of the saturated vertices (a vertex v i is saturated if |A T | i F 2 = g 2 i where |A T | i indicates the ith row of |A T |): λ * (v i ) > 0 if |A T | i F 2 = g(v i ) 2 , = 0 otherwise. (2.16)
The variables ν s and ν t are not constrained to be set to 0 and 1, only their difference is constrained to be equal to one. Thus their value range between a constant and a constant plus one. Let us call this constant δ, and without loss of generality one can consider that δ = 0. The term ν is at optimality a weighted indicator of the source/sink/saturated vertices partition:

ν * (v i ) =    0 + δ if v i ∈ S, a number between (0 + δ) and (1 + δ) if |A T | i F 2 = g(v i ) 2 , 1 + δ if v i ∈ T.
The expression (2.9) of the CCMF dual shows that the problem is equivalent to finding a minimum weighted cut defined on the nodes.

Finally, the "weighted cut" is recovered in (2.15), and the "smoothness term" is compatible with large variations of ν at the boundary of objects because of a large denominator (λ) in the contour area. An illustration of optimal λ and ν on an image is shown on Fig. 2.2.

Solving the Combinatorial Continuous Max-Flow problem.

When considering how to optimize the CCMF problem (2.4), the first key observation is that the constraints bound a convex feasible region. The real valued functions f i : R m → R defined for every node i = 1, . . . , n as

f i (F ) = |A T | i F 2 -g 2
i are non-negative quadratic, so convex functions. We define the vector f

(F ) of R n as f (F ) = [f 1 (F ), . . . , f n (F )] T .
Since the constraints are convex, the CCMF problem may be solved with a fast primal dual interior point method (see [START_REF] Boyd | Convex Optimization[END_REF]), which we now review in the specific context of CCMF.

Algorithm 1: Primal dual interior point algorithm

Data: The primal dual interior point (PDIP) algorithm iteratively computes the primal F and dual λ, ν variables so that the Karush-Kuhn-Tucker (KKT) optimality conditions are satisfied. This algorithm solves the CCMF problem (2.7) by applying Newton's method to a sequence of a slightly modified version of the KKT conditions. The steps of the PDIP procedure are given in Algorithm 1. Specifically for CCMF, the system M ∆ y = r system may be written

F ∈ R m = 0, λ ∈ (R * + ) n , ν ∈ R n , µ > 1, ǫ > 0.
  n i=1 λ i ∇ 2 f i (F ) Df (F ) T A -diag(λ)Df (F ) -diag(f (F )) 0 A T 0 0     ∆ F ∆ λ ∆ν   = -   r d = c + Df (F ) T λ + Aν rc = -diag(λ)f (F ) -(1/t) rp = A T F   (2.17)
with r d , r c , and r p representing the dual, central, and primal residuals. Additionally, the derivatives are given by

∇f i (F ) = 2|A T | i • F, Df (F ) =    ∇f 1 (F ) T . . . ∇f n (F ) T    , ∇ 2 f i (F ) = diag    2|A T | i1 . . . 2|A T | im    ,
with "•" denoting the Hadamard (element-wise) product between the vectors. Consequently, the primary computational burden in the CCMF algorithm is the linear system resolution required by (2.17). In the result section we present execution times obtained in our experiments using a conventional CPU implementation.

Observe that, although this linear system is large, it is very sparse and is straightforward to parallelize on a GPU [START_REF] Bolz | Sparse matrix solvers on the GPU: Conjugate gradients and multigrid[END_REF][START_REF] Kruger | Linear algebra operators for GPU implementation of numerical algorithms[END_REF][START_REF] Grady | Random walks for interactive organ segmentation in two and three dimensions: Implementation and validation[END_REF], for instance using an iterative GPU solver. If it does not necessarily imply a faster solution, the asymptotic complexity of modern iterative and modern direct solvers is about the same and, in our experience, there has been a strong improvement in the performance of a linear system solve when going from a direct solver CPU solution to a conjugate gradient solver GPU solution, especially as the systems get larger.

Comparison between CCMF and existing related approaches.

3.1. Min cut. As detailed in Section 2.3, there is a difference between the CCMF and the classic max flow formulation in the capacity constraint. Therefore, the dual problems are different. As proved by Ford and Fulkerson [START_REF] Ford | Maximal flow through a network[END_REF], and formalized for example in [START_REF] Bhusnurmath | Graph cuts via ℓ 1 norm minimization[END_REF], the min cut problem writes with our notations of section 2.3 min u gT |Au + c|.

(3.1)

We note that the ℓ 1 norm of the gradient of the solution u of the above min-cut expression turns into an ℓ 2 norm of the gradient of the solution ν in the CCMF dual (2.8).

3.2. Primal Dual Total variations. Unger et al. [START_REF] Unger | Interactive globally optimal image segmentation[END_REF] proposed to optimize the following primal dual TV formulation min

u max F i,j (u i+1,j -u i,j )F 1 i,j + (u i,j+1 -u i,j )F 2 i,j + data fidelity, s.t. (F 1 i,j ) 2 + (F 2 i,j ) 2 ≤ g i,j . (3.2)
We can note that the CCMF flow capacity constraint in a 4-connected lattice is similar to the constraint in (3.2) with a slight modification, and would be with finite-element notations

(F 1 i-1,j ) 2 + (F 2 i,j-1 ) 2 + (F 1 i,j ) 2 + (F 2 i,j ) 2 ≤ g i,j . (3.3)
In contrast to this finite element discretization, the provided CCMF formulation of continuous max-flow is defined arbitrary graphs. Furthermore, we note that the optimization procedure employed to solve (3.2) generalizes Appleton-Talbot's algorithm.

3.3. Combinatorial total variations. We now compare the dual CCMF problem with the Combinatorial Total Variation (CTV) problem. We note that the CCMF dual is different than CTV and discuss the weak duality of the two problems.

We recall the continuous expression of total variation given by Strang where Ω represents the image domain and Γ the boundary of Ω. The source and sink membership is represented by f , such that f (x, y) > 0 if (x, y) belongs to the sink, f (x, y) < 0 if (x, y) belongs to the source, and f = 0 otherwise. Considering a transport graph G, and a vector u defined on the nodes, this continuous problem may be written with combinatorial operators

min u∈R n g T |A T |(Au) 2 , s. t. u s -u t = 1, (3.5)
where the square root operator of a vector

v = [v 1 , . . . , v k ] is an element-wise square root operator √ v = [ √ v 1 , . . . , √ v k ].
Another way to write the same problem is min

u vi∈V g i eij ∈E (u i -u j ) 2 , s. t. u s -u t = 1. (3.6)
We note that in these equations the capacity g must be defined on the vertices. Although we describe this energy as "combinatorial total variation", due to its derivation in discrete calculus terms, it is important to note that this formulation is very similar to the discretization which has appeared previously in the literature from Gilboa and Osher [START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF], and Chambolle [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] (if we allow g = 1 everywhere).

3.3.1. CCMF and CTV are not dual. Strang proved that the continuous max-flow problem is dual to the total variation problem under some assumptions on g. Remarkably, this duality is not observed in the discrete case in which the CCMF dual and CTV problems are given by

min λ,ν λ T g 2 + 1 4 1 n • /(|A|λ) T (c + Aν) 2 , ≡ min u g T |A T |(Au) 2 , s. t. λ ≥ 0, s. t. u s -u t = 1.
We note that the analytic expressions are different and also not equivalent. The duality of the classical max-flow problem with the minimum cut holds because in the expression of the Lagrangian function, is possible to deduce a value of λ by substituting it into the dual problem so that it only depends on ν. However, λ in the dual CCMF problem (2.9) depends on several values of neighboring λ and ν. Thus, the CCMF dual problem cannot be simplified by removing a variable, for example by identifying ν and u in the CTV problem. Said differently, combining a null-divergence constraint with an anisotropic capacity constraint in the classical max-flow formulation allows the duality to hold; in contrast, in the CCMF formulation, combining a null-divergence constraint with an isotropic capacity constraint does not allow such duality to hold.

Numerically we can also show on examples that the value max F st of the flow optimizing CCMF is not equal to the minimum value of CTV (See Table 3.1).

Theoretical links between CCMF and CTV.

Even if CCMF is not dual with CTV, the two problems are weakly dual.

In the combinatorial setting the weak duality property is given by

||F || ≤ g ⇒ F T (Au) ≤ g T ||Au||. (3.7)
The next property shows that the norm ||F || = |A T |F 2 verifies weak duality. Proposition 3.1. Let G be a transport graph, F a flow in G verifying the capacity constraint |A T |F 2 ≤ g. Let u be a vector of R n (defined on nodes of G). Then

F T Au ≤ g T |A T |(Au) 2 .
Proof. Since we know that |A T |F 2 ≤ g, the following statement is true

|A T |F 2 T |A T |(Au) 2 ≤ g T |A T |(Au 2 .
We can now show that

F T Au ≤ (|A T |F 2 ) T |A T |(Au) 2 .

Transport graph with weights g on the nodes

Combinatorial Continuous Max-Flow (CCMF)

max F Fst, s. t. A T F = 0, |A T |F 2 ≤ g 2 .
Combinatorial Total Variation (CTV)

min u g T |A T |(Au)
1 0 0 1 1 0 0 1 0 0 0.2 Fst = 1.63 CT V (u) = i gi j (ui -uj ) 2 = 2 + 1 × √ 2 + 1 × √ 2 +1 × 1 + 0.8 2 +2 × 0.8 2 + 0.2 2 +2 × 0.2 = 8.16.
Table 3.1 Example illustrating the difference between optimal solutions of combinatorial continuous maximum flow problem (CCMF), and combinatorial total variation (CTV).

Using summation notation, then by the Cauchy-Schwartz inequality,

i∈V eij ∈E F ij (u i -u j ) ≤ i∈V   eij ∈E F ij 2     eij ∈E (u i -u j ) 2   .
We conclude that

F T Au ≤ (|A T |F 2 ) T |A T |(Au) 2 ≤ g T |A T |(Au) 2 .
In terms of energy value, this proposition means that the CCMF energy, that is to say the flow, is always smaller than the combinatorial total variation. Proposition 3.2. Let F be a compatible flow verifying the constraints in (2.4), and u a vector of R n such as u s -u t = 1. Then,

F st ≤ g T A T (Au) 2 . (3.8)
Proof. In the combinatorial setting, the Green formula gives

F T Au = u T A T F.
Let a be a vector of R n defined for each node v i as

a i =    -1 if v i belongs to the sink, 1
if v i belongs to the source, 0 otherwise.

(3.9)

We can provide the following equivalent formulation of the CCMF problem (2.4) using an incidence matrix A of the transport graph deprived of the sink-source edge:

max F st , s.t. A T F = F st a, |A T |F 2 ≤ g. (3.10)
As F verifies the divergence free constraint

A T F = F st a, u T A T F = u T F st a,
and as the equation u s -u t = 1 may be written equivalently a T u = 1, we conclude that

u T F st a = F st ,
and by weak duality (Property 3.1),

F st ≤ g T A T (Au) 2 .
This property should be of interest for extending CCMF to applications other than segmentation, which is outside the scope of this paper. In the next section, we present the performance of our approach in the context of segmentation.

Results

. We now present applications of Combinatorial Continuous Maximum Flow in image segmentation. To be used as a segmentation algorithm, three solutions are possible: the dual variable ν may be used directly if the user needs a matted result, otherwise ν may be thresholded, or finally an isoline or isosurface may be extracted from ν.

In the introduction we discussed how several works are related to CMF. Some are equivalent or slight modifications of AT-CMF for non-segmentation applications [START_REF] Unger | TVSeg -Interactive total variation based image segmentation[END_REF][START_REF] Unger | Interactive globally optimal image segmentation[END_REF][START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF]19]. As in the original AT-CMF work, none of these come with a convergence proof. In contrast, the work of [START_REF] Pock | An algorithm for minimizing the Mumford-Shah functional[END_REF] and [START_REF] Lippert | Discrete approximations to continuum optimal flow problems[END_REF] are provably convergent but both are very slow, and do not generalize easily to 3D image segmentation. Consequently, it seems reasonable to compare our segmentation method to the original AT-CMF as representative of the continuous approach, as well as graph cuts, which represent the purely discrete case. Finally, even if the total variation problem is different than the CMF problem (existence of a duality gap in the continuous [START_REF] Nozawa | Examples of max-flow and min-cut problems with duality gaps in continuous networks[END_REF]), the two problems are related enough to be compared. Specifically, we will also compare with Chambolle and Pock's recent work [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] which presented an algorithm for optimizing a TV-based energy for image segmentation.

Our validation is intended to establish three properties of the CCMF algorithm. First, we establish that the CCMF does avoid the metrication artifacts exhibited by conventional graph cuts (on a 4-connected lattice). This property is established by examples on a natural image and the recovery of the classical catenoid structure as the minimal surface spanning two rings. Second, we compare the convergence of the CCMF algorithm to the AT-CMF algorithm to show that the CCMF algorithm converges more quickly and in a more stable fashion. Finally, we establish that our formulation of the CMF problem does not degrade segmentation performance on a standard database. In fact, because of the reduction in metrication error our algorithm gives improved numerical results.For this experiment, we use the GrabCut database to compare the quality of the segmentation algorithms. In addition to the above tests, we demonstrate through examples that the CCMF algorithm is also flexible enough to incorporate prior (unary) terms and to operate in 3D.

4.1. Metrication artifacts and minimal surfaces. We begin by comparing the CCMF segmentation result with the classical max-flow algorithm (graph cuts). Figure 4.1 shows the segmentation of a brain, in which the contours obtained by graph cuts are noticeably blocky in the areas of weak gradient, while the contours obtained by both AT-CMF and CCMF are smooth.

In the continuous setting, the maximum flow computed in a 3D volume produces a minimal surface. The CCMF formulation may be also recognized as a minimal surface problem. In the dual formulation, the objective function is equivalent to a weighted sum of surface nodes. In [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF], Appleton and Talbot compared the surfaces obtained from their algorithm with the analytic solution of the catenoid problem to demonstrate that their algorithm was a good approximation of the continuous minimal surface and was not creating discretization artifacts. The catenoid problem arises from consideration of two circles with equal radius whose centers lie along the z axis. The minimal surface which forms to connect the two circles is known as a catenoid. The catenoid appears in nature, for example by creating a soap bubble between two rings. In order to demonstrate that CCMF is also finding a minimal surface, we performed the same catenoid experiment as was in [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF]. The results are displayed in Figure 4.2, where we show that CCMF approximates the analytical solution of the catenoid with even greater fidelity than the AT-CMF example.

In order to verify that the solution obtained by CCMF approximates perimeters of planar objects in 2D, we segmented several shapes -squares, discs -with known perimeters, scaled at different size, and compared the analytic perimeters with the energies obtained by CCMF. The results are reported in Fig. 4.3, where we observe that the true boundary length is proportional to the energy obtained with CCMF.

4.2. Stability, convergence and speed. We may compare the segmentation results using ν to Appleton-Talbot's result using P . We recall that AT-CMF solves the partial differential equation system (2.3) in order to solve the continuous maximum flow problem (2.1), but no proof was given for convergence. The potential function P approximates an indicator function, with 0 values for the background labels, and 1 for the foreground. It can be difficult to know when to stop the AT-CMF algorithm, since the iterations used to solve the AT-CMF algorithm may oscillate, as displayed in Figure 4.4 on a synthetic image. In contrast, the CCMF algorithm is guaranteed to converge and smoothly approaches the optimum solution.

We also compare segmentation results with results obtained by the recent algorithm of Chambolle and Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] optimizing a total variation based energy using a dual variable. In comparison to our work, the dual variable is not a flow (no constraint on the divergence and no proof of convergence toward a minimal surface as defined in [START_REF] Strang | Maximum flows through a domain[END_REF] or [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF]). In the literature, the continuous-domain duality between TV and max-flow is assumed, but in computational practice they are really different. Nevertheless, even as the method of Chambolle-Pock is solving a different problem, we performed numerical comparative tests. The problem of Chambolle and Pock (now denoted CP-TV) for binary segmentation optimizes min where A ∈ R m×n is the incidence matrix of the graph of n nodes and m edges defined on the image, g ∈ R n is the metric on the nodes, defined in equation (2.6).

u∈[0,1] n max F ∈R m ,||F ||∞≤1 F T (Au) + g T u + hard constraints attachment, (4.1)
The CCMF algorithm is faster than both AT-CMF and CP-TV for 2D image segmentation. We have implemented CCMF in Matlab. We used an implementation of Appleton-Talbot's algorithm in C++ provided by the authors, and a Matlab software implementing CP-TV on CPU, also provided by the respective authors. The three implementations make use of multi-threaded parallelization, and the CPU times reported here were computed on a Intel Core 2 Duo (CPU 3.00GHz) processor, with 2 Gb of RAM. The CCMF average computation time on a 321 × 481 image is 181 seconds after 21 iterations. For AT-CMF, 80000 iterations require 547 seconds. For CP-TV, the average number of iterations needed for convergence is 36000, requiring an average time of 1961 seconds on CPU, and the median computation time is 1406 seconds. Note that a GPU implementation of CP-TV, also provided by the authors, achieves a speed gain factor of about 100 on a NVidia Quadro 3700. Although it was not pursued here, it should be noted that the CCMF optimization approach could also fully benefit from parallelization (e.g., on a GPU) because the core computation is to solve a linear system of equations. Realizing the benefit of a GPU solution would require the use of an iterative algorithm such as conjugate gradients or multigrid [START_REF] Bolz | Sparse matrix solvers on the GPU: Conjugate gradients and multigrid[END_REF][START_REF] Kruger | Linear algebra operators for GPU implementation of numerical algorithms[END_REF], which would make comparisons to the direct solver used here less immediate. However, previous examples in the literature have suggested that this change in solver has not created difficulties in order to realize the benefits of a GPU architecture [START_REF] Grady | Random walks for interactive organ segmentation in two and three dimensions: Implementation and validation[END_REF]. Sometimes, it may not always be necessary to wait so long for the complete convergence of CCMF (or AT-CMF) to obtain acceptable results. In cases where images exhibiting sufficiently strong gradients, one iteration may be enough to obtain a satisfying segmentation. On such an image, one iteration of CCMF, and 100 iterations of AT-CMF show acceptable approximate results reached only after about 2 seconds for either algorithm. However, one cannot always rely on strong edges, and the power of our CCMF formulation is to behave in a predictable fashion even when edges are weak.

We may also compare the computation time of CCMF to the optimization of total variation using the Split Bregman method [START_REF] Goldstein | The split bregman method for l1-regularized problems[END_REF]. Optimizing TV on a 100 × 100 image requires 5000 iterations and takes 23 seconds with a sequential implementation of Split Bregman in C. On the same image, CCMF required only 17 iterations to reach the convergence criteria ||r d || < 1 and || η|| < 2, taking 4.7 seconds with a Matlab implementation (using a 2-threaded solver). We conclude that TV-based methods appear to be quite slow in the context of image segmentation, although we employed the Split Bregman algorithm which is known as one of the fastest algorithms for TV optimization for image denoising, and the very recent CP-TV algorithm. This difference in speed between the denoising and segmentation applications is due to the very large number of iterations required to convergence to a binary segmentation. Table 4.1 Dice coefficient (percentage) computed between the segmentation mask and the ground truth image (provided by GrabCut database). The lines above the double bar show results for the first set of seeds, and the lines below the double bar show results obtained with the second set of seeds.

Segmentation quality.

In this experiment, we compared our CCMF formulation to the AT-CMF formulation and conventional graph cuts on the problem of image segmentation, to determine if there were strong performance differences between the formulations. We expect that there would not be, since in principle all three formulations are trying to "minimize the cut", but have slightly different definitions of the cut length. The primary advantage that we expect from our formulation is in the reduction of metrication artifacts (as compared to conventional graph cuts) and speed/convergence (as compared to AT-CMF).

Our experiment consists of testing the Graph Cut, AT-CMF, CP-TV, and Combinatorial Continuous Max-Flow algorithms on a database with the same seeds. We used the Microsoft 'GrabCut' database available online [START_REF] Rother | GrabCut" -Interactive foreground extraction using iterated graph cuts[END_REF], which is composed of fifty images provided with seeds. From the seed images provided, it is possible to extract two different possible markers for the background seed. Examples of such seeds are shown in Fig. 4.5(a). Different convergence criteria are available for CCMF, such as the duality gap and norms of the residuals. However, for the CMF algorithm, we do not have any satisfying criteria. Bounding the number of non binary occurrences of P does not mean that the convergence is reached, because of possible oscillations. An intermediate result after ten thousand iterations may be significantly different from the result reached after one hundred thousand. Consequently, we have run the AT-CMF algorithm until we were convinced to have reached convergence, i.e. when P was nearly binary and did not change significantly when we doubled the number of iterations. For half of the images in the GrabCut database, AT-CMF algorithm required more than 20000 iterations to reach convergence, for a third of the images, more than 80000 iterations, and for 1/4 of the images, more than 160000 iterations. Binary convergence was still not reached after even 500000 iterations for the rest of the images (1/4), but we stopped the computation anyway. The Chambolle-Pock TV based algorithm for binary segmentation is provably convergent and benefits also from several possible convergence criteria. The convergence criteria used is a ratio of changed pixels from two successive intervals of iterations being smaller than an epsilon that we fixed to 10 -7 in order to obtain satisfying results. In practice, 1/4 of the images converged in less than 60000 iterations and for the rest of them we increased the maximum number of iterations up to 200000. In contrast, CCMF only needed 21 iterations on average, and never more than 27, to reach the convergence criteria ||r d || < 1 and || η|| < 2. We noted that for all methods, segmentation quality degraded quickly if these conditions were not respected.

Figure 4.3 displays the performance results for these algorithms. The Dice coefficient is a similarity measure between sets (segmentation and ground truth), ranging from 0 to 1.0 for no match and a perfect match respectively. All the tested algorithms show very good results, with a Dice coefficient of 0.95-0.97 for the well positioned seeds, and 0.89-0.92 for the second set of seeds, further away from the objects. The CCMF and AT-CMF results are really close, and the mean is better than the GC and the CP-TV results.

4.4. Extensions. The primary focus of our experiments were to demonstrate that our CCMF algorithm achieves a solution which does not exhibit metrication artifacts, that it is fast with provable convergence and that the segmentation quality is high. In the remainder of this section, we show that the algorithm can also incorporate unary terms and be equally applied in 3D. In fact, since CCMF is defined for an arbitrary graph, it could even be used to perform clustering in this more generalized framework. However, the benefit of avoiding gridding artifacts is less meaningful when performing clustering on an arbitrary graph, and therefore we omit any examples of this nature.

Unary terms.

A simple modification of the transport graph G permits the use of unary terms to automatically specify objects to be segmented. This approach is similar to the use of unary terms in the graph cuts computation [START_REF] Boykov | Interactive graph cuts for optimal boundary and region segmentation of objects in n-d images[END_REF]. The placement of unary terms for adding data attachment constraints in the max-flow problem is performed by adding edges linking every node of the lattice to the source and to the sink. In the case of the CCMF problem, as weights are defined on the nodes, we need to add intermediary nodes between the grid and source on the one hand, and the grid and sink on the other. However, since these intermediary nodes have just two edges incident on them, these node weights are equivalent to edge weights, meaning that our construction is equivalent to the use of unary terms for AT-CMF that was pursued by Unger et al. [START_REF] Unger | Interactive globally optimal image segmentation[END_REF]. In our case, we used weighted intermediary nodes simply to keep the CCMF framework consistent. Considering that the original lattice is composed of n nodes, we add for each node an "upper" node linked to the source and a "lower" node linked to the sink. An example of this construction is shown in Example of graph construction to incorporate unary terms for unsupervised segmentation. Each node is connected to a source node and sink node through an intermediary node which is weighted to reflect the strength of the positive and negative unary term. This construction is then applied with a simple appearance prior to automatically segment two images.

The weights of the additional nodes may be set to reflect the node priors for a particular application. For image segmentation, given mean background and foreground color BC and F C, we can set the capacities of the background nodes to g i = exp(-β(BC i -I i ) 2 ) and foreground nodes to g i = exp(-β(F C i -I i ) 2 ). Examples of results using these appearance priors are shown on Fig. 4.6.

Classification.

As the general CCMF formulation is defined on arbitrary graphs, it is possible to employ it in tasks beyond image segmentation. For instance CCMF may be employed to solve general problems of graph clustering, even on networks with no meaningful embedding such as social networks. We show a possible application in a classification example, considering the social network studied by Zachary [START_REF] Zachary | An information flow model for conflict and fission in small groups[END_REF]. After the split of a university karate club due to a conflict between its two leaders, Zachary's goal was to study if it was possible to predict the side joined by the members, based only on the social structure of the club. Classically, the graph is built by associating each member to a node, and edges link two members when special affinities are known between them. As weights are associated here to the strength of coupling between members, different strategies are possible for assigning the weights onto the nodes, which is necessary for using CCMF. For example, a given node/member may be assigned by the mean of affinity with its neighboring members in the graph. This weighting strategy works well in this example (See Fig. 4.7). Another strategy would be to compute the solution in the dual graph, that is in the graph where nodes are replaced by edges and (weighted) edges are replaced by (weighted) nodes. We have shown in an example the ability of CCMF for classification, a task that can not be treated with finite element-based methods such as AT-CMF. Evaluating the performances of CCMF for classification may be a topic for future research. 4.4.3. 3D segmentation. For 3D image segmentation, the minimal surface properties of CCMF generate good quality results, as shown in Fig 4 .8. The CCMF formulation applies equally well in 2D or 3D, since CCMF is formulated on an arbitrary graph (which may be a 2D lattice, a 3D lattice or an even more abstract graph). In 3D, our CCMF implementation is suffering from memory limitations in the direct Just one node is misclassified with CCMF. However, this node is known classically to be an unusual situation within the social network (see [START_REF] Zachary | An information flow model for conflict and fission in small groups[END_REF]) which is not captured by any known algorithm. solver we used, limiting its performances. Future work will address this issue using a dedicated solver.

Conclusion.

In this paper, we have presented a new combinatorial formulation of the continuous maximum flow problem and a solution using an interior point convex primal-dual algorithm. This formulation allows us to optimize the maximum flow problem as well as its dual for which we provide an interpretation as a minimal surface problem. This new combinatorial expression of continuous max-flow avoids blockiness artifacts compared to graph cuts. Furthermore, the formulation of CMF on a graph reveals that it is actually the fact that the capacity constraints are applied to point-wise norms of flow through nodes that allows us to avoid metrication errors, as opposed to the conventional graph cut capacity constraint through individual edges. Additionally, it was shown that our CCMF formulation provides a better approximation to the analytic catenoid than the conventional AT-CMF discretization of the continuous max-flow problem. Contrary to graph cuts, we showed that CCMF estimates the true boundary length of objects. Finally, unlike finite-element based methods such as AT-CMF, the CCMF formulation is expressed for arbitrary graphs, and thus can be employed in a large variety of tasks such as classification.

We provide in this paper an exact analytic expression of the dual problem, convergence of the algorithm being guaranteed by the convexity of the problem. In terms of speed, when an approximate solution is sufficient, our implementation of CCMF in Matlab is competitive to the Appleton-Talbot approach, which uses a system of PDEs, and a C++ implementation. The Appleton-Talbot algorithm has the significant drawback of not providing a criterion for convergence. In practice, this translates to long computation times when convergence is difficult for the AT-CMF algorithm. In contrast, our algorithm in this case is much faster, as we observe that the number of iterations required for convergence does not vary much for CCMF (less than a factor of two). The CCMF algorithm is simple to implement, and may be applied to arbitrary graphs. Furthermore, it is straightforward to add unary terms to perform unsupervised segmentation.

We have also compared our algorithm to known weighted combinatorial TV minimization methods (CTV optimized with split-Bregman, CP-TV). We have shown that our results are generally better for segmentation than combinatorial TV-based methods, and that our implementation on CPU is much faster and more predictable.

Specifically, the number of iterations required for convergence does not vary much for CCMF, whereas it can vary of more than a factor of ten for CP-TV. The deep study of the relationships between CCMF and combinatorial TV reveals that, in contrast with expectations from their duality under restrictive assumptions in the continuous domain, their duality relationship is only weak in the combinatorial setting. In fact max-flow and total variation are different problems with different constraints, yielding different algorithms and different results. One key difference between max-flow and total variation is that max-flow algorithms were developed as segmentation algorithms, following the graph cuts framework. One consequence is that max-flow formulations have a null divergence objective, which is not present in total variation formulations. Null divergence can be viewed as a consequence or a necessity in order to obtain constant partitions almost everywhere, in particular binary ones. This difference is important because in the proposed CCMF framework we impose a tight null divergence constraint, which to our knowledge is novel for an isotropic formulation (graph-cuts have always had this constraint). AT-CMF and similar frameworks achieve null divergence if and when convergence is achieved. There is no null divergence constraint at all in total variation frameworks. As a consequence, while in practice it is possible to compare total variation and max flow formulations, they should be treated differently. However, the strong computational performance and segmentation quality results of CCMF as compared to TV (using the strongest known optimizations methods) suggests that it may be advantageous to apply our CCMF formulation to problems for which TV has proven effective (such as filtering).

Several further optimizations of CCMFs are possible, for instance using multigrid implementations, the possibility to use GPU to solve the iterative linear system, the use of a dedicated solver for the particular sparse linear system involved in the computation. We also plan to compare the efficiency obtained with the interior point method to a first order algorithm for solving the CCMF problem.

Ultimately, we hope to employ combinatorial continuous maximum flow as a powerful segmentation algorithm which avoids metrication artifacts and provides a fast solution with provable convergence. Furthermore, we intend to explore the potential of CCMF to optimize other energy functions for which graph cuts or total variation have proved useful, such as surface reconstruction or efficient convex filtering.

  Fig.2.1. The difference between classical max-flow on a graph with the combinatorial continuous max-flow (CCMF) on a graph is that classical max-flow uses edge-weighted capacities while CCMF uses node-weighted capacities. This difference is manifest in the different solutions obtained for both algorithms and the algorithms required to find a solution. Specifically, the solution to the CCMF problem on a lattice does not exhibit metrication bias.

  Result: F solution to the CCMF problem (2.7) such as F st is maximized under the divergence free and capacity constraints, and ν, λ solution to the CCMF dual problem (2.9). repeat 1. Compute the surrogate duality gap η = -f (F ) T λ and set t = µn/ η. 2. Compute the primal-dual search direction ∆ y such as M ∆ y = r. 3. Determine a step length s > 0 and set y = y + s∆ y . (y = [F, λ, ν] T ) until ||r p || 2 ≤ ǫ, ||r d || 2 ≤ ǫ, and η ≤ ǫ.

Fig. 4 . 1 .

 41 Fig. 4.1. Brain segmentation. (a) Original image with foreground and background seeds. (b,c,d) Segmentation obtained with (b) graph cuts (GC), (c) AT-CMF, threshold of P (obtained after 10000 iterations), (d) CCMF, threshold of ν (15 iterations).

Fig. 4 . 2 .

 42 Fig. 4.2. The catenoid test problem: The source is constituted by two full circles and sink by the remaining boundary of the image. (a) Surface computed analytically, (b) isosurface of P obtained by CMF, (c) isosurface of ν obtained by CCMF. The root mean square error (RSME) has been computed to evaluate the precision of the results to the surface computed analytically. The RSME for CMF is 1.98 and for CCMF 0.75. The difference between those results is due to the fact that the CMF algorithm enforces exactly the source and sink points, leading to discretization around the disks. In contrast, the boundary localized around the seeds of ν is smooth, composed of grey levels. Thus the resulting isosurface computed by CCMF is more precise.

Fig. 4 . 3 .

 43 Fig. 4.3. Comparison of perimeters computed analytically and with Graph cuts and CCMF. To estimate the boundary length with CCMF, segmentations have been produced using the input images in (b) with the foreground and background seeds that appears in (a). (c) : plot showing the linear relation between the true perimeters and the energies obtained with CCMF. The symbols used to represent the dots correspond to the objects segmented in the images. A similar relation of proportionality is obtained with CP-TV. (d) : value of the cut obtained for a diagonal and a vertical line of same length using graph cuts. The difference observed explains (e) : plot showing the nonlinear relation between the true perimeters and the energies obtained with Graph cuts.

Fig. 4 . 4 .

 44 Fig. 4.4. Segmentation of an artificial image with AT-CMF (top row) and CCMF (bottom row). Top row (AT-CMF): (a) Image where the black and white discs are seeds. AT-CMF result stopped in (b) after 100 iterations, (c) 1000 iterations, (d) 10000 iterations. Bottom row (CCMF): (e) Image where the black and white discs are seeds, CCMF result ν after (f ) 1 iterations, (g) ν after 15 iterations iterations and (h) threshold of the final ν.

Fig. 4 . 5 .

 45 Fig. 4.5. (a) : Seed images belonging to the first (top) and second (bottom) set of seeds provided by the GrabCut database. (b) Segmentations by Graph Cuts computed with seeds from the two sets, (c) Segmentations by CCMF.

Fig. 4

 4 Fig. 4.6.

  Fig. 4.7. Zachary's karate club network. The two leaders in conflict are represented by the left hand side and right hand side nodes.Just one node is misclassified with CCMF. However, this node is known classically to be an unusual situation within the social network (see[START_REF] Zachary | An information flow model for conflict and fission in small groups[END_REF]) which is not captured by any known algorithm.

Fig. 4 . 8 .

 48 Fig. 4.8. Liver segmentation in 3D by CCMF.