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 [1] and then later level sets [2]. These formulations were used to optimize energies of varied sophistication (e.g., using regional, texture, motion or contour terms [3]) 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 [4, 5, 6].
Graphbased approaches
The maxflow/mincut problem on a graph is a classical problem in graph theory, for which the earliest solution algorithm goes back to Ford and Fulkerson [7]. Initial methods for global optimization of the boundary length of a region formulated the energy on a graph and relied on maxflow/mincut methods for solution [8, 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 3Dlungs example in [10]). 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 [11], but the large number of graph edges required to implement this solution could cause memory concerns when implemented for large 2D or 3D images [10].
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 [12, 10, 13]. Historically, a continuous maxflow (and dual mincut problem) was formulated by Iri [14] and then Strang [15]. 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 [10] provided the first PDEbased algorithm to find Strang’s continuous maxflows and therefore optimal mincuts. This same algorithm was also derived by Unger et al. from the standpoint of minimizing continuous total variation [16]. Adapted to image segmentation, this algorithm is shown to be equivalent [17] to the AppletonTalbot 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 [18], Zach et al [19] and Chambolle et al [20] 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 [21] a general saddle point algorithm that may be used in various applications, but needs half an hour to segment a 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 [22] has shown the continuous maxflow problem for the 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é [23] and later Rudin, Osher and Fatemi [24] 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 [12], Nesterov’s algorithm [25] and SplitBregman/DouglasRachford methods [26, 27]. 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 maxflow is dual with total variation in the continuous setting under restrictive regularity conditions [28]. In fact, Nozawa [29] showed that there is a duality gap between weighted TV and weighted maxflow 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 continuousdomain 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 maxflow problem illustrate two difficulties with continuousbased formulation, that are (1) the discretization step, which is necessary for deriving algorithms, but may break continuousdomain 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 [30]. 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 [31, 32] 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., [33, 34]).
In particular, Lippert presented in [35] a combinatorial formulation of an isotropic continuous maxflow 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 4connected squared grid, and the generalization to higher dimensions seems nontrivial. 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 lattice.
Motivation and contributions
In this paper, we pursue a combinatorial reformulation of the continuous maxflow problem initially formulated by Strang, which we term combinatorial continuous maximum flow (CCMF). Viewing our contribution differently, we adopt a discretization of continuous maxflow as the primary problem of interest, and then we apply fast combinatorial optimization techniques to solve the discretized version.
Our reformulation of the continuous maxflow problem produces a (divergencefree) flow problem defined on an arbitrary graph. Strikingly, CCMF are not equivalent to the discretization of continuous maxflows 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 maxflow on a graph. In particular, we will see that the difference lies in the fact that capacity constraints in classical maxflow 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 primaldual interior point method to solve it, such as interior point methods have been used for graph cuts [36] and second order cone problems in general [37]. The CCMF problem has several desirable properties, including:

the solution to CCMF on a 4connected lattice avoids the metrication errors. Therefore, the gridding error problems may be solved without the additional memory required to process classical maxflow on a highlyconnected graph;

in contrast to continuous maxflow algorithms of AppletonTalbot (ATCMF) and equivalent, the solution to the CCMF problem can be computed with guaranteed convergence in practical time;

the algorithm for solving the CCMF problem is fast, easy to implement, compatible with multiresolution 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:

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, ATCMF 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;

to expose links between CCMF and CTV: We prove that the weak duality holds between the two problems; and

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 maxflow, 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 maxflow (CMF) problem
First introduced by Iri [14], Strang presents in [15] an expression of the continuous maximum flow problem
(1) 
Here we denote by the total amount of flow going from the source location to the sink location , is the flow, and 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 [40]. In order to solve the problem (1), the AppletonTalbot algorithm (ATCMF) [10]
solves the following partial differential equation system:
(2)  
Here is a potential field similar to the excess value in the PushRelabel maximum flow algorithm [41]. ATCMF is effectively a simple continuous computational fluid dynamics (CFD) simulation with nonlinear constraints. It uses a forward finitedifference discretization of the above PDE system subject to a CourantFriedrichLevy (CFL) condition, also seen in early levelsets methods. At convergence, the potential function approximates an indicator function, with values for the background labels, and for the foreground, and the flow has zerodivergence 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 , which can be very slow.
Although AppletonTalbot is a continuous approach, applying this algorithm to image processing involves a discretization step. The capacity constraint is interpreted as , with the outgoing flow of edges linked to node along the axis, and the outgoing flow linked to node along the axis. We notice that the weights are associated with point locations (pixels), which will correspond later to a nodeweighted graph.
In the next section, we use the operators of discrete calculus to reformulate the continuous maxflow problem on a graph and show the surprising result that the graph formulation of the continuous maxflow leads to a different problem from the classical maxflow problem on a graph.
2.2 A discrete calculus formulation
Before establishing the discrete calculus formulation of the continuous maxflow problem, we specify our notation. A graph consists of a pair with vertices and edges . A transport graph comprises two additional nodes, named a source and a sink , and additional edges linking some nodes to the sink and some to the source. Including the source and the sink nodes, the cardinalities of are given by and . An edge, , spanning two vertices, and , is denoted by . 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, , is denoted by . In this work, we assume that and , and use
to denote the vector of
containing the for every edge of . In addition to edge weights, we may also assign weights to nodes. The weight of node is denoted by . In this work, we also assume that and . We use to denote the vector of containing the for every edge of . We define a flow through edge as where and use the vector to denote the flows through all edges in the graph . Each edge is assumed to be oriented, such that a positive flow on edge indicates the direction of flow from to , while a negative flow indicates the direction of flow from to .The incidence matrix is a key operator for defining a combinatorial formulation of the continuous maxflow problem. Specifically, the incidence matrix is known to define the discrete calculus analogue of the gradient, while is known to define the discrete calculus analogue of the divergence (see [32] 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
(3) 
for every vertex and edge . In our formulation of continuous maxflow, we use the expression 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 maxflow of (1) on a transport graph. As in [32, 33, 34], the continuous vector field indicating flows may be represented by a vector on the edge set, . Additionally, the combinatorial divergence operator allows us to write the first constraint in (1) as . The second constraint in (1) involves the comparison of the pointwise norm of a vector field with a scalar field. Therefore, we can follow [32, 33, 34] to define the pointwise norm of the flow field as . In our notations here, as in the rest of the paper, denotes a elementwise product, “” denoting the Hadamard (elementwise) 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
(4) 
Compare this formulation to the classical maxflow problem on a graph, given in our notation as
(5) 
By comparing these formulations of the traditional maxflow with our combinatorial formulation of the continuous maxflow, it is apparent that the key difference between the classical formulation and our combinatorial continuous maxflow 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 contrastweighting applied to nodes (pixels) has been a feature of several algorithms with a continuous formulation, including geodesic active contours [40], CMF [10], TVSeg [16], and [21]. On the other hand, the problem defined in (4) is defined on an arbitrary graph in which contrast weights (capacities) are almost always defined on edges. The nodeweighted capacities fit Strang’s formulation of a scalar field of constraints in the continuous maxflow 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 maxflow problem. Figure 1 illustrates the relationship of the edge capacity constraints in the classical maxflow problem to the node capacity constraints in the CCMF formulation. The nulldivergence 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.
In the context of image segmentation, the vector varies inversely to the image gradient. We propose to use, as in [16],
(6) 
where indicates the image intensity. For simplicity, this weighting function is defined for greyscale images, but 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 maxflow problem is dual to the mincut 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 nodeweighted minimum cut. In the following proposition, we denote the elementwise quotient of two vectors and by . We also denote a unit vector [1, …, 1] of size , and recall that the square exponent of a vector represents the resulting vector of the elementwise mulitiplication .
In a transport graph with edges, nodes, we define a vector of , composed of zeros except for the element corresponding to the source/sink edge which is . Let and be two vectors of . The CCMF problem
s. t.  
has for dual
(8)  
s. t. 
equivalently written in (15), and the optimal solution verifies
(9) 
and the following equalities
(10) 
The Lagrangian of (2.3) is
with and two Lagrange multipliers. We have to find in the dual function is such as .
(11) 
Substituting in the Lagrangian function , we obtain
The Lagrangian may be simplified as:
The dual of (2.3) is also
(12)  
s. t. 
Furthermore, the KKT optimality conditions give
(13) 
Using the expression of from equation (11) , and substituting it in (13), we obtain
Taking the sum of right hand side and left hand side,
(14) 
Finally, if we substitute (14) in the dual expression (8), we have
The expression of the CCMF dual may be written in summation form as
(15) 
Interpretation: The optimal value is a weighted indicator of the saturated vertices (a vertex is saturated if where indicates the th row of ):
(16) 
The variables and 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 . The term is at optimality a weighted indicator of the source/sink/saturated vertices partition:
The expression (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 (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.4 Solving the Combinatorial Continuous MaxFlow problem
When considering how to optimize the CCMF problem (4), the first key observation is that the constraints bound a convex feasible region.
The real valued functions defined for every node as are nonnegative quadratic, so convex functions. We define the vector of as .
Since the constraints are convex, the CCMF problem may be solved with a fast primal dual interior point method (see [42]), which we now review in the specific context of CCMF.
The primal dual interior point (PDIP) algorithm iteratively computes the primal and dual , variables so that the KarushKuhnTucker (KKT) optimality conditions are satisfied. This algorithm solves the CCMF problem (2.3) 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 system may be written
(17) 
with , , and representing the dual, central, and primal residuals. Additionally, the derivatives are given by
with “” denoting the Hadamard (elementwise) product between the vectors.
Consequently, the primary computational burden in the CCMF algorithm is the linear system resolution required by (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 [43, 44, 45], 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.
3 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 [7], and formalized for example in [36], the min cut problem writes with our notations of section 2.3
(18) 
We note that the norm of the gradient of the solution of the above mincut expression turns into an norm of the gradient of the solution in the CCMF dual (8).
3.2 Primal Dual Total variations
Unger et al. [17] proposed to optimize the following primal dual TV formulation
(19) 
We can note that the CCMF flow capacity constraint in a 4connected lattice is similar to the constraint in (19) with a slight modification, and would be with finiteelement notations
(20) 
In contrast to this finite element discretization, the provided CCMF formulation of continuous maxflow is defined arbitrary graphs. Furthermore, we note that the optimization procedure employed to solve (19) generalizes AppletonTalbot’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
(21) 
where represents the image domain and the boundary of . The source and sink membership is represented by , such that if belongs to the sink, if belongs to the source, and otherwise. Considering a transport graph , and a vector defined on the nodes, this continuous problem may be written with combinatorial operators
(22) 
where the square root operator of a vector is an elementwise square root operator . Another way to write the same problem is
(23) 
We note that in these equations the capacity 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 [46], and Chambolle [12] (if we allow everywhere).
3.3.1 CCMF and CTV are not dual
Strang proved that the continuous maxflow problem is dual to the total variation problem under some assumptions on . Remarkably, this duality is not observed in the discrete case in which the CCMF dual and CTV problems are given by
,  
, 
We note that the analytic expressions are different and also not equivalent. The duality of the classical maxflow 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 (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 in the CTV problem. Said differently, combining a nulldivergence constraint with an anisotropic capacity constraint in the classical maxflow formulation allows the duality to hold; in contrast, in the CCMF formulation, combining a nulldivergence constraint with an isotropic capacity constraint does not allow such duality to hold.
Numerically we can also show on examples that the value of the flow optimizing CCMF is not equal to the minimum value of CTV (See Table 1).
Transport graph with weights on the nodes 
Combinatorial Continuous MaxFlow (CCMF)

Combinatorial Total Variation (CTV)





3.3.2 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
(24) 
The next property shows that the norm verifies weak duality.
Let be a transport graph, a flow in verifying the capacity constraint . Let be a vector of (defined on nodes of ). Then
Since we know that , the following statement is true
We can now show that
Using summation notation, then by the CauchySchwartz inequality,
We conclude that
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.
Let be a compatible flow verifying the constraints in (4), and a vector of such as . Then,
(25) 
In the combinatorial setting, the Green formula gives
Let be a vector of defined for each node as
(26) 
We can provide the following equivalent formulation of the CCMF problem (4) using an incidence matrix of the transport graph deprived of the sinksource edge:
(27) 
As verifies the divergence free constraint ,
and as the equation may be written equivalently , we conclude that
and by weak duality (Property 3.3.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.
4 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 ATCMF for nonsegmentation applications [16, 17, 18, 19]. As in the original ATCMF work, none of these come with a convergence proof. In contrast, the work of [47] and [35] 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 ATCMF 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 [29]), the two problems are related enough to be compared. Specifically, we will also compare with Chambolle and Pock’s recent work [21] which presented an algorithm for optimizing a TVbased 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 4connected 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 ATCMF 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 maxflow algorithm (graph cuts). Figure 3 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 ATCMF 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 [10], 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 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 [10]. The results are displayed in Figure 4, where we show that CCMF approximates the analytical solution of the catenoid with even greater fidelity than the ATCMF example.



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 CPTV. (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.
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. 5, 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 AppletonTalbot’s result using . We recall that ATCMF solves the partial differential equation system (2.1) in order to solve the continuous maximum flow problem (1), but no proof was given for convergence. The potential function approximates an indicator function, with values for the background labels, and for the foreground. It can be difficult to know when to stop the ATCMF algorithm, since the iterations used to solve the ATCMF algorithm may oscillate, as displayed in Figure 6 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 [21] 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 [15] or [10]). In the literature, the continuousdomain duality between TV and maxflow is assumed, but in computational practice they are really different. Nevertheless, even as the method of ChambollePock is solving a different problem, we performed numerical comparative tests. The problem of Chambolle and Pock (now denoted CPTV) for binary segmentation optimizes
(28) 
where is the incidence matrix of the graph of nodes and edges defined on the image, is the metric on the nodes, defined in equation (6).
The CCMF algorithm is faster than both ATCMF and CPTV for 2D image segmentation. We have implemented CCMF in Matlab. We used an implementation of AppletonTalbot’s algorithm in C++ provided by the authors, and a Matlab software implementing CPTV on CPU, also provided by the respective authors. The three implementations make use of multithreaded 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 image is 181 seconds after 21 iterations. For ATCMF, 80000 iterations require 547 seconds. For CPTV, 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 CPTV, 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 [43, 44], 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 [45].
Sometimes, it may not always be necessary to wait so long for the complete convergence of CCMF (or ATCMF) 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 ATCMF 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 [26]. Optimizing TV on a 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 and , taking 4.7 seconds with a Matlab implementation (using a 2threaded solver). We conclude that TVbased 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 CPTV 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.
4.3 Segmentation quality
In this experiment, we compared our CCMF formulation to the ATCMF 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 ATCMF).
Our experiment consists of testing the Graph Cut, ATCMF, CPTV, and Combinatorial Continuous MaxFlow algorithms on a database with the same seeds. We used the Microsoft ‘GrabCut’ database available online [48], 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. 7(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 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 ATCMF algorithm until we were convinced to have reached convergence, i.e. when was nearly binary and did not change significantly when we doubled the number of iterations. For half of the images in the GrabCut database, ATCMF 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 ChambollePock 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 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 and . We noted that for all methods, segmentation quality degraded quickly if these conditions were not respected.
Dice Coeff.  GC  ATCMF  CPTV  CCMF  

First set of seeds  mean  95.2  94.9  95.2  95.3 
median  96.2  96.1  96.3  96.3  
stand. dev.  4.3  4.1  4.1  4.1  
Second set of seeds  mean  89.3  89.7  89.1  89.5 
median  91.2  92.0  91.7  92.3  
stand. dev.  8.4  7.7  8.5  8.6 
Figure 2 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 ATCMF results are really close, and the mean is better than the GC and the CPTV 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.
4.4.1 Unary terms
A simple modification of the transport graph 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 [8]. The placement of unary terms for adding data attachment constraints in the maxflow 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 ATCMF that was pursued by Unger et al. [17]. In our case, we used weighted intermediary nodes simply to keep the CCMF framework consistent. Considering that the original lattice is composed of 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 Fig. 8.
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 and , we can set the capacities of the background nodes to and foreground nodes to . Examples of results using these appearance priors are shown on Fig. 8.
4.4.2 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 [49]. 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. 9). 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 elementbased methods such as ATCMF. 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 10. 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 solver we used, limiting its performances. Future work will address this issue using a dedicated solver.
5 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 primaldual 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 maxflow 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 pointwise 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 ATCMF discretization of the continuous maxflow problem. Contrary to graph cuts, we showed that CCMF estimates the true boundary length of objects. Finally, unlike finiteelement based methods such as ATCMF, 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 AppletonTalbot approach, which uses a system of PDEs, and a C++ implementation. The AppletonTalbot 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 ATCMF 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 splitBregman, CPTV). We have shown that our results are generally better for segmentation than combinatorial TVbased 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 CPTV. 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 maxflow and total variation are different problems with different constraints, yielding different algorithms and different results. One key difference between maxflow and total variation is that maxflow algorithms were developed as segmentation algorithms, following the graph cuts framework. One consequence is that maxflow 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 (graphcuts have always had this constraint). ATCMF 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.
References
 [1] Kass, M., Witkin, A., Terzopoulos, D.: Snakes: Active contour models. IJCV 1 (1988) 321–331
 [2] Sethian, J.A.: Level Set Methods and Fast Marching Methods. Cambridge University Press (1999)
 [3] Paragios, N.: Geodesic Active Regions and Level Set Methods : Contributions and Applications in Artificial Vision. PhD thesis, INRIA Sophia Antipolis (2000)
 [4] Boykov, Y., Veksler, O., Zabih, R.: Fast approximate energy minimization via graph cuts. IEEE Transactions on Pattern Analysis and Machine Intelligence 23 (2001) 1222–1239
 [5] Ishikawa, H.: Exact optimization for Markov random fields with convex priors. IEEE Transactions on Pattern Analysis and Machine Intelligence 25 (2003) 1333–1336
 [6] Darbon, J., Sigelle, M.: Image restoration with discrete constrained total variation part ii: Levelable functions, convex priors and nonconvex cases. J. of Math. Imaging and Vis. 26 (2006) 277–291
 [7] Ford, L.R., Fulkerson, D.R.: Maximal flow through a network. Canadian Journal of Mathematics 8 (1956) 399–404
 [8] Boykov, Y., Jolly, M.P.: Interactive graph cuts for optimal boundary and region segmentation of objects in nd images. ICCV’01 1 (2001) 105–112 vol.1
 [9] Kolmogorov, V., Zabin, R.: What energy functions can be minimized via graph cuts? PAMI 26 (2004) 147–159
 [10] Appleton, B., Talbot, H.: Globally minimal surfaces by continuous maximal flows. PAMI 28 (2006) 106–118
 [11] Boykov, Y., Kolmogorov, V.: Computing geodesics and minimal surfaces via graph cuts. ICCV’03 (2003) 26–33
 [12] Chambolle, A.: An algorithm for total variation minimization and applications. J. Math. Imaging Vis. 20 (2004) 89–97
 [13] Nikolova, M., Esedoglu, S., Chan, T.F.: Algorithms for finding global minimizers of image segmentation and denoising models. SIAM JAM 66 (2006) 1632–1648
 [14] Iri, M.: Theory of flows in continua as approximation to flows in networks. Survey of Mathematical Programming (1979)
 [15] Strang, G.: Maximum flows through a domain. Math. Prog. (1983) 123–143
 [16] Unger, M., Pock, T., Trobin, W., Cremers, D., Bischof, H.: TVSeg  Interactive total variation based image segmentation. In: BMVC. (2008)
 [17] Unger, M., Pock, T., Bischof, H.: Interactive globally optimal image segmentation. Technical Report ICGTR08/02, Graz University of Technology, Gratz, Austria (2008)
 [18] Pock, T., Schoenemann, T., Graber, G., Bischof, H., Cremers, D.: A convex formulation of continuous multilabel problems. In: ECCV ’08. (2008) 792–805
 [19] Zach, C., Gallup, D., Frahm, J.M., Niethammer, M.: Fast global labeling for realtime stereo using multiple plane sweeps. In: Vis. Mod. Vis. (2008)
 [20] Chambolle, A., Cremers, D., Pock, T.: A convex approach for computing minimal partitions. Tech. Rep. 649, Ecole Polytechnique CMAP (2008)
 [21] Chambolle, A., Pock, T.: A firstorder primaldual algorithm for convex problems with applications to imaging. CMAP685 (2010)
 [22] Strang, G.: Maximum flows and minimum cuts in the plane. Journal of Global Optimization (2008)
 [23] D. Shulman, J.H.: Regularization of discontinuous flow fields. Proceedings Workshop on Visual Motion (1989) 81–86
 [24] Rudin, L.I., Osher, S., Fatemi, E.: Nonlinear total variation based noise removal algorithms. Phys. D 60 (1992) 259–268
 [25] Nesterov, Y.: Gradient methods for minimizing composite objective function. CORE Discussion Papers 2007076, Université catholique de Louvain, Center for Operations Research and Econometrics (CORE) (2007)
 [26] Goldstein, T., Osher, S.: The split bregman method for l1regularized problems. SIIMS 2 (2009) 323–343
 [27] Zhang, X., Burger, M., Bresson, X., Osher, S.: Bregmanized nonlocal regularization for deconvolution and sparse reconstruction (2009) [Submitted].
 [28] Nozawa, R.: Maxflow mincut theorem in an anisotropic network. Osaka J. Math 27 (1990) 805–842
 [29] Nozawa, R.: Examples of maxflow and mincut problems with duality gaps in continuous networks. Mathematical Programming 63 (1994) 213–234 10.1007/BF01582067.
 [30] Fefferman, C.L.: Existence and smoothness of the NavierStokes equation (2000) http://www.claymath.org/millennium/NavierStokes_Equations/navierstokes.pdf.
 [31] Hirani, A.N.: Discrete exterior calculus. PhD thesis, Cal. Tech (2003)
 [32] Grady, L.J., Polimeni, J.R.: Discrete Calculus: Applied Analysis on Graphs for Computational Science. Springer (2010)
 [33] Elmoataz, A., Lézoray, O., Bougleux, S.: Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing. IEEE Trans. on Image Processing 17 (2008) 1047–1060
 [34] Grady, L., Alvino, C.: The piecewise smooth MumfordShah functional on an arbitrary graph. IEEE Trans. on Image Proc. 18 (2009) 2547–2561
 [35] Lippert, R.A.: Discrete approximations to continuum optimal flow problems. Studies in Applied Mathematics 117 (2006) 321–333
 [36] Bhusnurmath, A., Taylor, C.J.: Graph cuts via norm minimization. IEEE Transactions on Pattern Analysis and Machine Intelligence 30 (2008) 1866–1871
 [37] Goldfarb, D., Yin, W.: Secondorder cone programming based methods for total variation image restoration. SIAM J. on Scientific Computing 27 (2005) 622–645
 [38] Buades, A., Coll, B., Morel, J.M.: A nonlocal algorithm for image denoising. CVPR 2005 2 (2005) 60–65
 [39] Grady, L., Schwartz, E.L.: Faster graphtheoretic image processing via smallworld and quadtree topologies. In: Proc. of CVPR. Volume 2., IEEE (2004) 360–365
 [40] Caselles, V., Kimmel, R., Sapiro, G.: Geodesic active contours. IJCV 22 (1997) 61–79

[41]
Goldberg, A.V., Tarjan, R.E.:
A new approach to the maximum flow problem.
STOC ’86, Theory of computing (1986) 136–146
 [42] Boyd, S., Vandenberghe, L.: Convex Optimization. CUP (2004)
 [43] Bolz, J., Farmer, I., Grinspun, E., Schröder, P.: Sparse matrix solvers on the GPU: Conjugate gradients and multigrid. In: ACM Transactions on Graphics. Volume 22 of SIGGRAPH. (2003) 917–924
 [44] Kruger, J., Westermann, R.: Linear algebra operators for GPU implementation of numerical algorithms. In: SIGGRAPH. (2003) 908–916
 [45] Grady, L., Schiwetz, T., Aharon, S., Westermann, R.: Random walks for interactive organ segmentation in two and three dimensions: Implementation and validation. In: Proc. In’l Conf. Medical Image Computing and ComputerAssisted Intervention. (2005) 773–780
 [46] Gilboa, G., Osher, S.: Nonlocal linear image regularization and supervised segmentation. Multiscale Modeling and Simulation 6 (2007) 595–630
 [47] Pock, T., Cremers, D., Bischof, H., Chambolle, A.: An algorithm for minimizing the MumfordShah functional. In: Int. Conf. Comp. Vis. (2009)
 [48] Rother, C., Kolmogorov, V., Blake, A.: “GrabCut” — Interactive foreground extraction using iterated graph cuts. In: SIGGRAPH. (2004) 309–314
 [49] Zachary, W.: An information flow model for conflict and fission in small groups. Journal of Anthropological Research 33 (1977) 452–473
Comments
There are no comments yet.