Camille Couprie

Leo Grady

Laurent Najman

Hugues Talbot

Power Watershed: A Unifying Graph-Based Optimization Framework

Keywords: Combinatorial optimization, Image segmentation, Graph cuts, Random walker, Shortest paths, Optimal Spanning Forests, Markov Random Fields !

In this work, we extend a common framework for graph-based image segmentation that includes the graph cuts, random walker, and shortest path optimization algorithms. Viewing an image as a weighted graph, these algorithms can be expressed by means of a common energy function with differing choices of a parameter q acting as an exponent on the differences between neighboring nodes. Introducing a new parameter p that fixes a power for the edge weights allows us to also include the optimal spanning forest algorithm for watershed in this same framework. We then propose a new family of segmentation algorithms that fixes p to produce an optimal spanning forest but varies the power q beyond the usual watershed algorithm, which we term power watershed. In particular when q = 2, the power watershed leads to a multilabel, scale and contrast invariant, unique global optimum obtained in practice in quasi-linear time. Placing the watershed algorithm in this energy minimization framework also opens new possibilities for using unary terms in traditional watershed segmentation and using watershed to optimize more general models of use in applications beyond image segmentation.

G RAPH-BASED segmentation algorithms have be- come quite popular and mature in recent years. The modern variations on graph-based segmentation algorithms are primarily built using a small set of core algorithms -graph cuts, random walker and shortest paths, which are reviewed shortly. Recently these three algorithms were all placed into a common framework that allows them to be seen as instances of a more general seeded segmentation algorithm with different choices of a parameter q [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF]. In addition to these algorithms, the ubiquitous watershed segmentation algorithm [12] shares a similar seeding interface but only recently was a connection made between the watershed algorithm and graph cuts [28]. In this paper, we show how this connection between watershed and graph cuts can be used to further generalize the seeded segmentation framework of [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF] such that watershed, graph cuts, random walker and shortest paths may all be seen as special cases of a single general seeded segmentation algorithm. Our more general formulation has several consequences which form our contributions.

1) This more general formulation reveals a previously unknown family of segmentation algorithms which we term power watershed. In this paper, we give an algorithm for solving the energy minimization problem associated with the power watershed and demonstrate that this new algorithm has the speed of standard watershed but performing almost as well or better as all the other algorithms on our benchmark segmentation tests. 2) Placing watershed in the same framework as graph cuts, random walker and shortest paths allows us to easily incorporate data (unary) terms into conventional watershed segmentation. 3) By placing the watershed algorithm in the same generalized framework as graph cuts, random walker and shortest paths, it is possible to take advantage of the vast literature on improving watershed segmentation to also improve these other segmentation approaches. 4) Defining an energy function for the watershed optimization allows us to provide an MRF interpretation for the watershed. 5) By incorporating unary terms, we can push watershed beyond image segmentation into the area of general energy minimization algorithms which could be applied to any number of applications for which graph and MRF models have become standard. Before proceeding to the exposition of our technique, we first review the graph-based segmentation literature in more detail.

A SHORT REVIEW OF GRAPH-BASED SEG-

MENTATION

The algorithms that are reviewed in this section view the image as a graph with each pixel corresponding to a node and edges weighted to reflect changes in image intensity, color or other features.

Watershed: There exist many possible ways for defining a watershed [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF], [67], [71], [START_REF] Bertrand | On topological watersheds[END_REF], [28], [START_REF] Cousty | Watershed cuts: thinnings, shortest-path forests and topological watersheds[END_REF]. Intuitively, the watershed of a function (seen as a topographical surface) is composed of the locations from which a drop of water could flow towards different minima. The framework allowing the formalization and proof of this statement is the optimal spanning forests relative to the minima [START_REF] Cousty | Watershed cuts. In 7th international symposium on mathematical morphology[END_REF], [28]. For the purpose of seeded image segmentation, the gradient of the image can be considered as a relief map and, instead of minima, seeds may be placed by the user or found automatically to specify the segmentation of the image into desired regions. If the gradient is inverted, the maxima are considered instead of minima, and a thalweg is computed instead of watershed. A thalweg is the deepest continuous line along a valley. In the rest of the paper, we use by convention the term "watershed" instead of "thalweg".

A maximum spanning forest (MSF) algorithm computes trees spanning all the nodes of the graph, each tree being connected to exactly one connected seed component, and the weight of the set of trees being maximum. If the seeds correspond to the maxima, the segmentation obtained by MSF is a watershed [28]. An optimal spanning forest can be computed by Kruskal's or Prim's algorithm [START_REF]On the shortest spanning tree of a graph and the traveling salesman problem[END_REF], [START_REF] Prim | Shortest connection networks and some generalizations[END_REF] among others in quasilinear time. In Kruskal's algorithm, the edges are sorted by decreasing edge weight, and chosen in that order to be added to the forest if they do not create cycles or join trees that are connected to different maxima.

Watersheds are widely used in image segmentation because there exist numerous and efficient algorithms that are easy to implement. However, segmentation results from watershed may suffer from leaks and degeneracy of the solution on the plateaus of the weight function.

Graph cuts: The labeling produced by the graph cuts (GC) algorithm is determined by finding the minimum cut between the foreground and background seeds via a maximum flow computation. The original work on GC for interactive image segmentation was produced by Boykov and Jolly [17], and this work has been subsequently extended by several groups to employ different features [START_REF] Blake | Interactive image segmentation using an adaptive GMMRF model[END_REF] or user interfaces [START_REF] Rother | GrabCut" -Interactive foreground extraction using iterated graph cuts[END_REF], [57]. Although GC is relatively new, the use of minimal surfaces in segmentation has been a common theme in computer vision for a long time [36], [15], [63] and other boundary-based user interfaces have been previously employed [62], [33], [22], [START_REF] Grady | Minimal surfaces extend shortest path segmentation methods to 3D[END_REF]. Two concerns in the literature about the original GC algorithm are metrication error ("blockiness") and the shrinking bias. Metrication error was addressed in subsequent work on GC by including additional edges [START_REF] Boykov | Computing geodesics and minimal surfaces via graph cuts[END_REF], by using continuous max flows [7] or total variation [START_REF] Unger | TVSeg -Interactive total variation based image segmentation[END_REF]. These methods for addressing metrication error successfully overcome the problem, but may incur greater memory and computation time costs than the application of maximum flow on a 4-connected lattice. The shrinking bias can cause overly small object segments because GC minimizes boundary length. Although some techniques have been proposed for addressing the shrinking bias [START_REF] Boykov | Computing geodesics and minimal surfaces via graph cuts[END_REF], [7], [START_REF] Vicente | Graph cut based image segmentation with connectivity priors[END_REF], these techniques all require additional parameters or computation.

Random walker: The random walker (RW) algorithm [START_REF] Grady | Random walks for image segmentation[END_REF] is also formulated on a weighted graph and determines labels for the unseeded nodes by assigning the pixel to the seed for which it is most likely to send a random walker. This algorithm may also be interpreted as assigning the unlabeled pixels to the seeds for which there is a minimum diffusion distance [23], as a semi-supervised transduction learning algorithm [31] or as an interactive version of normalized cuts [START_REF] Shi | Normalized cuts and image segmentation[END_REF], [START_REF] Grady | Fast approximate random walker segmentation using eigenvector precomputation[END_REF]. Additionally, popular image matting algorithms based on quadratic minimization with the Laplacian matrix may be interpreted as employing the same approach for grouping pixels, albeit with different strategies to determine the edge weighting function [START_REF] Levin | A closed form solution to natural image matting[END_REF]. Diffusion distances avoid segmentation leaking and the shrinking bias, but the segmentation boundary may be more strongly affected by seed location than with graph cuts [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF].

Shortest paths (geodesics):

The shortest path algorithm assigns each pixel to the foreground label if there is a shorter path from that pixel to a foreground seed than to any background seed, where paths are weighted by image content in the same manner as with the GC and RW approaches. This approach was recently popularized by Bai and Sapiro [START_REF] Bai | A geodesic framework for fast interactive image and video segmentation and matting[END_REF], but variants of this idea have appeared in other sources [START_REF] Criminisi | GeoS: Geodesic image segmentation[END_REF], [START_REF] Alvino | Efficient segmentation based on Eikonal and diffusion equations[END_REF], [START_REF] Falcão | The image foresting transformation[END_REF]. The primary advantage of this algorithm is speed and prevention of a shrinking bias. However, it exhibits stronger dependence on the seed locations than the RW approach [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF], is more likely to leak through weak boundaries (since a single good path is sufficient for connectivity) and exhibits metrication artifacts on a 4-connected lattice.

All of the above models may be considered as addressing energies comprised of only unary and pairwise (binary) energy terms. However, recent literature has found that the addition of energy terms defined on higher-order cliques can help improve performance on a variety of tasks [START_REF] Kohli | P3 & beyond: Solving energies with higher order cliques[END_REF], [START_REF] Kohli | Robust higher order potentials for enforcing label consistency[END_REF]. Although we do not address higher-order cliques specifically in this work, we note that all recent progress in this area has been through an equivalent construction of pairwise terms. Therefore, our results could also be useful in that context. Despite the recent popularity of energies defined on higher order cliques, pairwise terms (and watershed) are still used ubiquitously in the computer vision literature and any improvement to these models can have a broad impact.

An earlier conference version of this work appeared in [24].

A UNIFYING ENERGY MINIMIZATION FRAME-

WORK

We begin our exposition by reviewing the unity framework of [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF] before showing how to further broaden this framework to provide a general seeded segmentation

H H H H q p 0 finite ∞ 1 Collapse to seeds Graph cuts Power watershed q = 1 2 ℓ 2 norm Voronoi Random walker Power watershed q = 2 ∞ ℓ 1 norm Voronoi ℓ 1 norm Voronoi Shortest Path Forest

TABLE 1

Our generalized scheme for image segmentation includes several popular segmentation algorithms as special cases of the parameters p and q. The power watershed are previously unknown in the literature, but may be optimized efficiently with a maximum spanning forest calculation.

scheme that includes the maximum spanning forest algorithm for watershed as a special case. Examination of the special cases of this general algorithm reveals a new class of watershed segmentation models. We prove several theoretical properties of this new class of watershed and then give an algorithm for minimizing the energy associated with this generalized watershed model.

A review of the existing generalized segmentation framework

In this section, we review the segmentation framework introduced by Sinop and Grady in [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF]. A graph consists of a pair G = (V, E) with vertices v ∈ V and edges e ∈ E ⊆ V × V with cardinalities n = |V | and m = |E|. An edge, e, spanning two vertices, v i and v j , is denoted by e ij . In image processing applications, each pixel is typically associated with a node and the nodes are connected locally via a 4 or 8-connected lattice.

A weighted graph assigns a real value to each edge called a weight. In this work, the weights are assumed to be non-negative. The weight of an edge e ij is denoted by w(e ij ) or w ij . We also denote w F i and w Bi as the unary weights penalizing foreground and background affinity at node v i . In the context of segmentation and clustering applications, the weights encode nodal affinity such that nodes connected by an edge with high weight are considered to be strongly connected and edges with a low weight represent nearly disconnected nodes. One common choice for generating weights from image intensities is to set

w ij = exp(-β(∇I) 2 ), (1)
where ∇I is the normalized gradient of the image I. The gradient for a grey level image is I i -I j . Details on the parameters used are given in the experimental section. We use w to denote the vector of R m that contains the weights w ij of every edge e ij in G.

The generalized energy proposed in [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF] is given by as

min x eij ∈E (w ij |x i -x j |) q + vi∈V (w i |x i -y i |) q (2)
where y represents a measured configuration and x represents the target configuration. In this equation w ij can be interpreted as a weight on the gradient of the target configuration, such that the first term penalizes any unwanted high-frequency content in x and essentially forces x to vary smoothly within an object, while allowing large changes across the object boundaries.

The second term enforces fidelity of x to a specified configuration y, w i being weights enforcing that fidelity.

For an image segmentation in two classes, given foreground F and background B seeds, (2) may be included in the following algorithm

Step 1: x = arg min

x eij ∈E (w ij |x i -x j |) q + vi (w F i |x i |) q + vi (w Bi |x i -1|) q , s.t. x(F ) = 1, x(B) = 0,
Step 2:

s i = 1 if x i ≥ 1 2 , 0 if x i < 1 2 . (3)
In other words, we are looking for an optimum x ⋆ of Eq. (3) that may be interpreted as a probability for a given pixel to belonging to either the foreground or the background, the final decision (hard segmentation) s giving the segmentation being taken by a threshold.

It was shown in [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF] that graph cuts gives a solution to this model when q = 1, random walker gives the solution to this model when q = 2 and shortest paths (geodesics) give a solution to this model as q → ∞. The case of this model with a fractional q was optimized in [START_REF] Singaraju | P-brush: Continuous valued MRFs with normed pairwise distributions for image segmentation[END_REF] via reweighted least squares and shown that intermediate values of q allowed for an algorithm which "interpolated" between the graph cuts, random walker or shortest paths algorithms.

In related work, Strang showed in [START_REF] Strang | l 1 and l ∞ approximation of vector fields in the plane[END_REF] that minimization of the ℓ p norm of the gradients of a potential field with boundary conditions (in continuous space with real-valued potentials) also leads to (continuous) maxflow (for an ℓ 1 norm of the gradients), the Dirichlet problem (for an ℓ 2 norm) and shortest paths (for an ℓ ∞ norm). Therefore, the framework of [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF], which we now extend, may be seen as presenting similar ideas defined on an arbitrary graph, using the bridge between continuous PDEs and graph theory provided by discrete calculus [START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF].

Broadening the framework to watershed

We now broaden the segmentation algorithm in (3) to include watershed simply by separating the exponent on the weights and the variables. Specifically, we introduce parameter p to define a new segmentation model as

Step 1: x = arg min

x eij ∈E w p ij |x i -x j | q + vi w p F i |x i | q + vi w p Bi |x i -1| q , s.t. x(F ) = 1, x(B) = 0,
Step 2:

s i = 1 if x i ≥ 1 2 , 0 if x i < 1 2 . (4)
As before, the final segmentation s being chosen via a threshold.

We observe that (4) can be formulated in a general manner by rewriting it as the minimization of a general energy function E p,q (x) by introducing auxiliary nodes (See [44] for more details):

min x λ eij ∈E w p ij |x i -x j | q smoothness terms + vi∈V w p i |x i -y i | q
data fidelity terms (5) For example, the unary term w p Bi |x i -1| q can also be rewritten as w p i |x iy i | q , where y i is an auxiliary node and the signal at this auxiliary node is fixed at y i = 1.

As with (3), when p is a small finite value, then the various values of q may be interpreted respectively as the graph cuts (q = 1), and random walker (q = 2) algorithms. When q and p converge toward infinity with the same speed, then a solution to (4) can be computed by the shortest path (geodesics) algorithm. Those three algorithms form the underpinning for many of the advanced image segmentation methods in the literature.

It was shown in [START_REF] Allène | Some links between min cuts, optimal spanning forests and watersheds[END_REF], [START_REF] Allène | Some links between extremum spanning forests, watersheds and min-cuts[END_REF] that when q = 1 (graph cuts) and p → ∞ then the solution of ( 4) is given by a maximum spanning forest algorithm. Said differently, as the power of the weights increases to infinity, then the graph cuts algorithm produces a segmentation corresponding to a segmentation by maximum spanning forest. Interpreted from the standpoint of the Gaussian weighting function in (1), it is clear that we may associate β = p to understand that the watershed equivalence comes from operating the weighting function in a particular parameter range. An important insight from this connection is that above some value of β we can replace the expensive max-flow computation with an efficient maximum spanning forest computation. By raising p → ∞ and varying the power q we obtain a previously unexplored family of segmentation models which we refer to as power watershed. An important advantage of power watershed with varying q is that the main computational burden of these algorithms depends on an MSF computation, which is extremely efficient [21]. In the next sections we explore two cases that are, to the best of our knowledge, unexplored. First, we show that case p finite, q → ∞ corresponds to a Voronoi diagram computation from the seeds. Second, we prove that when q is finite, as p → ∞ there exists a value of p after which any of the algorithms (regardless of q) may be computed via an MSF. We then give an algorithm to minimize (4) for any value of q when p → ∞. Table 1 gives a reference for the different algorithms generated by various value of p and q.

The case p finite, q → ∞: Voronoi diagram

Intuitively, we see that when the power over the neighboring differences tends toward infinity, the weights become negligible so that the problem obtained from ( 4) is a Voronoi diagram of the seeds.

A proof showing that solving the minimization problem (4) when p = q and q → ∞ can be achieved by shortest path computations is given in [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF]. Here we use the same idea to prove that the problem (4) in the case p finite, q → ∞ is equivalent to a Voronoi diagram problem.

As q √ . is monotonic, minimizing E p,q is equivalent to minimizing q E p,q . First we may factorize the objective function of our problem (4)

q eij ∈E w p ij |x i -x j | q = q eij ∈E w ij p q |x i -x j | q (6)
Taking the limit lim q→∞ q i X q i of a q-norm yields the maximum norm max i X i .

Therefore, our objective function may be written

lim q→∞ q eij ∈E w p ij |x i -x j | q = lim q→∞ max eij ∈E w ij p q |x i -x j |.
(7) The minimization problem can be written as

min x max eij ∈E lim q→∞ w ij p q |x i -x j |, s.t. x(F ) = 1, x(B) = 0. (8)
When q → ∞ and p is finite, p q converges toward 0, so w ij p q converges toward 1 for every edge of E. Also the case p finite, q → ∞ can be brought back to the case p = 0, q → ∞ which solution is a Voronoi diagram with an ℓ 1 norm (due to the assumed 4-connectivity of the lattice).

The case q finite, p → ∞ leading to watershed

We now generalize the link between GC and MSF established by Allène et al. [START_REF] Allène | Some links between min cuts, optimal spanning forests and watersheds[END_REF], [START_REF] Allène | Some links between extremum spanning forests, watersheds and min-cuts[END_REF] by proving that GC, RW, and generally all cuts resulting of the minimization of E p,q converge to MSF cuts as p tends toward infinity under the condition that all the maxima of the weight function are seeded.

The following properties are presented in the special case of segmentation into two classes, given two sets of labeled nodes F and B. However the following results generalize easily to multilabel segmentation.

Definitions 1: (q-cut, MSF, MSF cut) In a graph G, let F and B be two disjoint nonempty sets of nodes, p and q two real positive values, and s the 4) on the subset (considering the merged nodes as a unique node) with q = 2 (i.e., solution of the combinatorial Dirichlet problem), (e) Another plateau connected to three labeled vertices is encountered, and (f) Final solutions x and s obtained after few more steps. The q-cut, which is also an MSF cut, is represented in dashed lines.

segmentation result defined in Eq 4. The set of edges e ij such that

s i = s j is a q-cut. Let Y be a subgraph of G. We say that Y is an extension of F ∪ B if each connected component of Y contains exactly one vertex of F ∪ B and each vertex of F ∪ B is contained in a connected component of Y .
Consequently it is possible to define a label l on each vertex of Y , 0 to the vertices connected to a vertex of B, and 1 to the vertices connected to a vertex of F .

Examples of extensions appear in Fig. 1, where F and B are displayed in (a), and two possible extensions in bold in (b) and (c), with their corresponding labels.

Let F be a subgraph of G. We say that F is a spanning forest (relative to

F ∪ B) if: (i) F is an extension of F ∪ B, (ii) F contains no cycles, and (iii) V (F) = V (F is spanning all vertices of G).
The weight w F of a forest F for w is the sum of the weight of all edges belonging to F:

w F = eij ∈F w ij .
We say that a spanning forest F is a maximum spanning forest (MSF) for w if the weight of F is maximum, i.e. greater or equal to the weight of any other spanning forest.

Let F be an MSF for w, and l its associated label. An MSF cut for w is the set of edges e ij such that l i = l j .

We call a subgraph M a maximum of w if M is connected, all the edges of M have the same weight w M , and the weight of any edge adjacent to M is strictly lower than w M .

Finally, a plateau is a subgraph of G consisting of a maximal set of nodes connected with edges having the same weight.

Those definitions are compatible with the watershed cut framework of [28]. We may now introduce a general link between the Maximum Spanning Forest segmentation result and the solution of the optimization of (4) when the power of the weights converges toward infinity.

Theorem 1: Let M be the subgraph of G composed of the union of all maxima of the weight function w. If every connected component of M contains at least a vertex of B ∪ F , and q ≥ 1, then any q-cut when p → ∞ is an MSF cut for w.

Proof: The proof is based on the construction of a set of edges that belong to the q-cut when p → ∞. During the construction, we consider the edges of E in decreasing order, following Kruskal's algorithm for maximum spanning forest construction. At the end of the construction, the q-cut obtained is an MSF cut for w. The successive steps of the proof are illustrated on an example on Fig. 1.

At each step, we consider the set E max of edges of maximum weight w max . We normalize all the weights by dividing them by w max , to obtain all the weights between 0 and 1 with the normalized weight of E max equal to 1. The energy to minimize is also

eij ∈E w ij w max p |x i -x j | q , s.t. x(F ) = 1, x(B) = 0. (9)
As all maxima of the weight function contain seeds, each connected component of E max has at least one labeled vertex. For every connected component C max of E max , two cases are possible:

If C max contains no vertices of different labels, the edges of weight w max can not be a part of the minimum qcut energy when p tends toward infinity because all the other normalized weights converge toward 0 and so does any finite sum of these weights. Choosing x i = x j for all edges e ij ∈ C max is the only possibility to eliminate the terms of maximum weight of [START_REF] Audigier | Uniquely-determined thinning of the tie-zone watershed based on label frequency[END_REF]. The edges of C max are not included in the q-cut, and also do not belong to the MSF cut as they have to be merged to labeled nodes to form an MSF (e.g., Figure 1(b)).

If C max contains vertices of different labels, any labeling can be done on the plateau, because adding edges of C max to the q-cut or not will always give an MSF cut on the plateau (e.g., Figures 1(c) and 1(d)).

Repeating the steps recursively until all the vertices are labeled, we find that in building a q-cut, we are also building an MSF cut for w in exactly the same manner as with Kruskal's algorithm.

p = 1 0 0.4 0.6 1 1 1 2 p → ∞ 0 .5 -ǫ .5 + ǫ 1 1 1 2 (a) q-cuts for different p p → ∞ 0 0.5 0.5 1 1 2 1 (b) Cut on limp→∞ x
Fig. 2. Let x = arg min E pq . When the maxima of the weight function are not seeded, the threshold of lim p→∞ x can be different of the limit lim p→∞ of the threshold of x. (a) Labeling x = arg min E pq and corresponding qcut (q = 2) when the weights are at the power p = 1, and below for an arbitrary big value of p. The q-cut in dashed line remains in the center. (b) Labeling x = arg min lim p→∞ E pq and cut (dashed) corresponding to the threshold of x. In this example, the q-cut is not an MSF cut, justifying the condition in Thm. 1. Note that the threshold of lim p→∞ x is an MSF cut, as stated later in Prop. 2.

In theorem 1, the condition for seeds to be the maxima of the weight function is necessary as shown in Fig. 2.

We can note that if the weights are all different, the MSF cut is unique and Thm. 1 is also true without the condition for seeds to be the maxima of the weight function.

The next property states that when the power on the neighboring node differences is strictly greater than one, the minimization of E p,q admits a unique solution.

Property 1: If q is a real number such that 1 < q < ∞, then the solution x to problem (4) is unique.

Proof: Let A be the incidence matrix of the graph G, and x a vector of R n + . We note by | • | the elementwise absolute value operator. The function g : x → Ax is convex. The function h : x → |x| q is convex and nondecreasing. The function f : x → w T x is also convex and non-decreasing. Note that E p,q (x) can be written the following way:

E p,q (x) = f • h • g(x) = w pT |Ax| q . ( 10
)
As h is a non-decreasing convex function, and g is convex, h • g is convex. As f is a non-decreasing convex function, E p,q is convex. If 1 < q < ∞, the function h • g is strictly convex, so E p,q is a strictly convex function, thus the minimization of E p,q subject to the boundary constraints is achieved by a unique x.

Before introducing in Section 3 an algorithm to compute the solution x to the optimization of E p,q when p → ∞, we present an interpretation of the minimization of our general energy as a maximum a posteriori approximation.

Interpretation as a Markov Random Field

An optimum x ⋆ of Eq. 5 may be interpreted as a probability for each pixel to belong to the object (as opposed to background). More rigorously, as shown in [78], this segmentation model can be viewed as an estimation of a continuous valued MRF. By linking the watershed algorithm to this framework, it becomes possible to view the watershed algorithm as the MAP estimation of an MRF. However, note that this analysis allows us to interpret the watershed as the MAP estimate for a particular MRF, which is in contrast to previous efforts to link a probabilistic framework with the watershed (such as [6], who use random placements of seeds to define the most probable locations of the watershed lines).

In this section, we follow the development of [78], with modifications to incorporate the power watershed.

In the interpretation as an MRF, we define the binary segmentation label s i for node v i as a Bernoulli random variable (i.e., s i = 1 if v i is foreground and s i = 0 if v i is background), in which the variable x i denotes the success probability for the distribution of s i , i.e., p(s i = 1|x i ). In this case, the success probability may be written as

p(s i = 1|x i ) = max{min{x i , 1}, 0} =      1 if x i > 1, x i if 0 ≤ x i ≤ 1, 0 if x i < 0.
(11) However, the generalized mean value theorem in [78] guarantees that the optimal solution x ⋆ to (5), assuming that the auxiliary nodes have values comprised between 0 and 1, takes its values between 0 and 1 when the weights are all positive valued. Consequently, in our context, we may simply set p(s i = 1|x i ) = x i without concern that x i will be outside the interval [0, 1].

Our goal is now to infer the hidden variables x i from the image content I. The hidden variables may be estimated in a Bayesian framework by considering the posterior model

p(x, s|I) ∝ p(x)p(s|x)p(I|s) = p(x) vi∈V p(s i |x i ) vi∈V p(I i |s i ), (12)
in which p(x) models how the parameters of the Bernoulli variables vary spatially. The spatial smoothness prior is parameterized by

p(x) ∝ exp   -λ eij ∈E w p ij |x i -x j | q   , (13)
where λ > 0 and the weights are strictly positive. We can estimate the marginalized MAP criterion to obtain the optimum x * by setting

x * = arg max x p(x)p(I|x) = arg max p(x) s p(I|s)p(s|x). ( 14
)
Unfortunately, s p(I|s)p(s|x), is not straightforward to estimate. Therefore, we assume that we can parameterize p(I i |x i ) as

p(I|x) ∝ exp - vi∈V w p i0 |x i -0| q - vi∈V w p i1 |x i -1| q ,
(15) where w i0 ≥ 0 and w i1 ≥ 0, and these terms act to bias the parameters x i toward 0 and 1. Similarly, these terms can be used to encode the user interaction (seeding) by setting a foreground seed v i to have weights (w i0 , w i1 ) = (0, ∞) and a background seed to have weights (w i0 , w i1 ) = (∞, 0). With this parameterization, then the MAP estimate described in ( 14) is equal to our energy minimization problem from (5).

While the use of binary variables (the s variable in our formulation) is more common in recent work which applies MRFs to image segmentation, our focus on estimating a real-valued parameter or variable is far from unique in the computer vision literature. For example, in Gaussian MRFs the variables each have a Gaussian distribution and the goal is often to estimate the (real-valued) parameters of these variables (i.e., mean and/or variance). These kind of MRFs have been applied in image segmentation and other pattern recognition applications [START_REF] Kaufhold | Energy Formulations of Medical Image Segmentations[END_REF], [5], [START_REF] Tappen | Learning Gaussian conditional random fields for low-level vision[END_REF]. Beyond Gaussian MRFs, anisotropic diffusion has been interpreted as a continuous-valued MRF [START_REF] Michael | Image statistics and anisotropic diffusion[END_REF], [START_REF] Krajsek | Diffusion filtering without parameter tuning: Models and inference tools[END_REF] and MRFs requiring continuous-valued estimations have appeared in both early work on computer vision [55], [37], [35], [16] and also recently [START_REF] Li | Learning for optical flow using stochastic optimization[END_REF], [START_REF] Samuel | Learning optimized MAP estimates in continuously-valued MRF models[END_REF], [START_REF] Schmidt | A generative perspective on MRFs in low-level vision[END_REF].

We now introduce an algorithm to optimize E p,q when p → ∞, and show that the threshold s of that solution produces an MSF cut.

ALGORITHM FOR OPTIMIZING THE CASE q FINITE, p → ∞

The algorithm proposed in this section may be seen as Kruskal's algorithm for maximum spanning tree with two main differences -a forest is computed in place of a tree, and the optimization

min x eij ∈ plateau |x i -x j | q (16)
is performed on the plateaus (the maximal set of nodes connected with edges of same weight). The power watershed algorithm is detailed in Alg. 1, and an illustration of different steps on an example is given in Fig. 3.

In Algorithm 1, the merge operation of a set of nodes S consists of removing the nodes in S from the graph and replacing these nodes with a single node such that any edge spanning a node in S to nodes in S now connects the merged node to the same nodes in S. Additionally, in the above algorithm, the unary terms in (4) are treated as binary terms connected to phantom seeds v F and v B , Algorithm 1: power watershed algorithm, optimizing p → ∞, q ≥ 1 Data: A weighted graph G(V, E) and a set of foreground F and background B seeds Result: A potential function x and a labeling s associating a label to each vertex. Set x F = 1, x B = 0 and all other x values as unknown. Sort the edges of E by decreasing order of weight. while any node has an unknown potential do Find an edge (or a plateau) E MAX in E of maximal weight; denote by S the set of nodes connected by E MAX . if S contains any nodes with known potential then Find x S minimizing (4) (using the input value of q) on the subset S with the weights in E MAX set to w ij = 1, all other weights set to w ij = 0 and the known values of x within S fixed to their known values. Consider all x S values produced by this operation as known. else Merge all of the nodes in S into a single node, such that when the value of x for this merged node becomes known, all merged nodes are assigned the same value of x and considered known.

Set

s i = 1 if x i ≥ 1 2 and s i = 0 otherwise. i.e., vi w p F i |x i -0| q + vi w p Bi |x i -1| q = vi w p F i |x i -x B | q + vi w p Bi |x i -x F | q . ( 17
)
We prove in the next section that the labeling x obtained by Algorithm 1 optimizes [START_REF] Alvino | Efficient segmentation based on Eikonal and diffusion equations[END_REF].

An illustration for this section is given in Figure 4. The segmentation was performed with progressively larger values of p, keeping q = 2 and shows that the segmentation result converges to the result given by the above algorithm for the power watershed with q = 2. The value q = 2 was employed for this example since it is known that q = 2 forces a unique minimum to (4) regardless of the value of p.

An implementation of Algorithm 1 when q = 2 can be downloaded from sourceforge [1].

Justification of the power watershed algorithm

We now prove that the algorithm we propose optimizes the energy presented in our framework when q > 1 and p → ∞.

Let us define the labeling x * as the solution x * = arg min x E p,q (x) defined in (4) subject to the boundary constraints. We note the labeling obtained by Algorithm 1 by x.

The two following theorems 2 and 3 state that the energy of the solution computed by the power watershed algorithm converges to the energy which minimizes E p,q when p → ∞.

Theorem 2: Let p, q be real positive numbers. Let w M be the maximum weight of the graph G. For every δ > 0, there exists a real k such that if p > k,

0 ≤ E p,q (x) w M p - E p,q (x * ) w M p ≤ δ. ( 18
)
The proof of this theorem is given in the appendix. Theorem 3: If q > 1, the potential x * obtained by minimizing the energy of (4) subject to the boundary constraints converges toward the potential x obtained by Algorithm 1. as p → ∞.

Proof: We prove that by optimizing (4) we are performing the same steps as Alg. 1. As in Thm. 1, at each step, we consider a set of connected edges of maximum weight E max of E, and we normalize all the weights, minimizing also [START_REF] Audigier | Uniquely-determined thinning of the tie-zone watershed based on label frequency[END_REF].

If E max contains no vertices of different labels, then the weights w max can not be a part of the minimum energy when p tends toward infinity, because all the other normalized weights converge toward 0 and so does any finite sum of these weights. Choosing x i = x j for every edge e ij = e max ∈ E max is the only possibility to eliminate the only term(s) of maximum weight of ( 9). This choice of xi = xj is also performed by Algorithm 1 by the "merge" operation. From the standpoint of energy minimization, having x i = x j in the graph G may be brought back to having one unique node instead of v i , v j , and e ij . We can also replace v i and v j by a unique node.

If E max contains vertices of a different label, as the weights of E max are arbitrarily greater than the weights of the unprocessed edges, minimizing (9) boils down to minimizing

eij ∈Emax |x i -x j | q , ( 19
)
with boundary conditions given by already labeled nodes. It is exactly what is performed by Algorithm 1 in the "If" part.

Repeating the steps recursively until all the vertices are labeled, we find that the Algorithm 1 procedure agrees with the energy minimization of (5).

We can note that even if Algorithm 1 minimizes the energy E p,q in the case p → ∞, several solutions x are possible when q = 1.

Property 2: For any q ≥ 1, the cut C defined by the segmentation s computed by Algorithm 1 is an MSF cut for w.

Proof: At each step of Algorithm 1, we consider a set of connected edges of maximum weight E max .

If E max contains no vertices of different labels, Algorithm 1 chooses x i = x j for the edges e ij ∈ E max . The edges of E max are not included in C, and also do not belong to the MSF cut as they have to belong to an MSF since their weight is maximum.

If E max contains vertices of different labels, any labeling can be done on the plateau, because adding edges of E max to the q-cut or not will always give an MSF cut on the plateau.

Repeating the steps of Algorithm 1 recursively until all the vertices are labeled, we find that we are building an MSF cut for w.

Using mathematical morphology for an efficient preprocessing step

One difficulty in Algorithm 1 is dealing with the set of merged nodes. More precisely, when solving (16) we need to keep track of which nodes have merged (with some nodes merged multiple times). If we look informally at the "emergence" process underlying the algorithm, it will help us to locate those maximal merged nodes. Using topographical references, we view the weights as the surface of a terrain, with the weight of an edge corresponding to its altitude. If the surface were completely covered by water, and the level of water slowly decreases, then islands (regional maxima) would appear that grow and merge. At a given level, when an island that does not contain a seeded pixel meets an island containing one, we can give a value to the (maximal) merged node. Indeed, we can see that any merged node consists of a connected component of a upper-level set of the weights. More precisely, let λ ∈ R + and w be the weight function defined on E. We define

w[λ] = {e ∈ E|w(e) ≥ λ}. (20)
The graph induced by w[λ] is called a section of w.

A connected component of a section w[λ] is called a component of w (at level λ).

The components of w can be used to find merged nodes.

Property 3: Any maximal merged node corresponds to a component of w that:

• does not contain any seed; • and is not contained in a larger unseeded component of w. Conversely, any component of w satisfying these two properties corresponds to a maximal merged node in Algorithm 1. The components of w, ordered by the inclusion relation, form a tree called the max-tree [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF] or the component tree [START_REF] Guillataud | Contribution l'analyse dendroniques des images[END_REF], [START_REF] Hanusse | Sémantique des images par analyse dendronique[END_REF], [START_REF] Breen | Attribute openings, thinnings, and granulometries[END_REF]. Several efficient algorithms exist to compute the component tree, some quasi-linear [START_REF] Najman | Building the component tree in quasilinear time[END_REF] (based on union-find [START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF]) and some parallelized [START_REF] Wilkinson | Concurrent computation of attribute filters on shared memory parallel machines[END_REF], [58]. From Prop. 3, it is easy to see how to use this tree in Algorithm 1. Note that such a tree, which keeps track of all components, can be used when one wants to improve a given segmentation result by adding extra seeds.

Another tool from mathematical morphology [START_REF] Najman | Mathematical Morphology: from theory to applications[END_REF] has been used as a preprocessing step for watershed segmentation with markers. It is called geodesic reconstruction from the markers [START_REF] Meyer | Morphological segmentation[END_REF], [12], and is given as a function w R such that, for every edge e, we set w R (e) to be equal to the level δ of the highest component of w containing e and at least one seed node. Note that any component of w R contains at least one seed.

Property 4: Any maximal merged node corresponds to a connected set of edges e ij that belong to a plateau of w R and that satisfy w ij > w R (e ij ). The converse is also true.

Prop. 4 also suggests that geodesic reconstruction can be used as a preprocessing in Alg. 1. Note that there exist some very efficient and easy to implement algorithms to compute a geodesic reconstruction [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms[END_REF], [START_REF] Najman | Geodesic Saliency of Watershed Contours and Hierarchical Segmentation[END_REF], [START_REF] Géraud | Algorithms for mathematical morphology[END_REF]. Both the component tree and the geodesic reconstruction have the same theoretical complexity, so either approach could be used profitably to reduce the bookkeeping necessary to keep track of merged nodes.

Prop. 4 also suggests links between our framework and the classical watershed-based segmentation framework [12], [START_REF] Meyer | Morphological segmentation[END_REF], [START_REF] Meyer | Segmentation, minimum spanning tree and hierarchies[END_REF]. The framework of watershed cuts [28], [START_REF] Cousty | Watershed cuts: thinnings, shortest-path forests and topological watersheds[END_REF] allows us to make a precise statement about this connection. The cut provided by a maximum spanning forest with one different seed for every maxima is called a watershed cut. Since geodesic reconstruction removes all maxima which are not connected to a seed, then we can state the following:

Property 5: Any q-cut is a watershed cut of the reconstructed weights. This statement ties the cuts produced by our power watershed framework to the concept of watershed cuts in the morphology literature.

Uniqueness of solution

Most of the energy minimization problems in our framework, i.e. the cases optimized by graph cuts, shortest path forests and maximum spanning forest algorithms (and watershed in general [12], [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF], [START_REF] Bieniek | An efficient watershed algorithm based on connected components[END_REF], [71], [34]) have the problem that the optimum solution may not be unique, for example on plateaus. That implies that the result of each one of these algorithms depends on the implementation.

To remove such dependency, two approaches have been proposed:

• A classical approach is to compute a geodesic distance on the plateau [71], and to use that distance as a way to distinguish between points of the plateau. Generally, the cut is located on the "middle of the plateau", but other locations are possible according to the application [START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF], [67]. • Another proposal is the tie-zone watershed [START_REF] Audigier | Uniquely-determined thinning of the tie-zone watershed based on label frequency[END_REF]; it takes into account all the possible solutions derived from a shortest-path based watershed to generate a unique solution: when the multiple solutions disagree with each other on the segmentation result of a region (i.e., the label to be assigned), the region is included to the tie-zone and a specific tie value is assigned to each node, corresponding to the probability of assigning a label to the node according to the number of all possible assignments. A major drawback of that tie-zone approach is that nodes with equal probability of belonging to different label classes can appear.

In contrast to that approaches, the power watershed compute a probability map (consisting of x in ( 4)) minimizing a global energy function, and whenever q is finite and q > 1, the solution is unique.

RESULTS

Generality of the framework

Adding unary terms

We now present an application of the framework to unseeded segmentation. Unary terms were first employed with graph cuts in [44]. Since this initial work, many other applications have used graph cuts with unary terms. Gathering watershed and graph cuts into the same framework allows us to employ unary terms for watershed computation. The unary terms in ( 4) are treated as binary terms connected to phantom seeds v F and v B as in (17).

For the example of image segmentation with two labels, the weights w Bi between v B and v i can be fixed to the absolute difference of the pixel v i intensity with the mean of the gray scales plus the variance, and w F i to the absolute difference of the pixel v i intensity with the mean of the gray scales minus the variance. An example of such a weighted graph is given in Fig. 5. With this construction, we can apply any of the algorithms in our framework to the resulting graph. An example of result is shown at Fig. 6 for purposes of segmenting blood cells. Note that those examples show how to add two phantom seeds, but this idea is extendable to more than two labels as explained in the next section. To the best of our knowledge, this is the first time that the watershed algorithm has been used as an unseeded segmentation method (i.e., without markers or seeds).

Multilabel segmentation

Minimizing exactly the energy E 1,1 is possible by using the graph cuts algorithm in the case of two labels, but is NP-hard if constraints impose more than two different labels. However, the other algorithms presented in our framework can perform seeded segmentation with as many labels as desired efficiently.

We detail the method of multilabel segmentation in the case of the power watershed algorithm. Let N represent the number of different labels l = 1, 2..., N . Instead of computing an x solution of the Foreground/Background as it is done for the two-labels segmentation, N solutions x l have to be computed. In order to perform N -labels segmentation, we may define seeds at a node i by setting x l i = 1 for a given label l and x l = 0 for any other label than l.

The segmentation result is obtained by affecting each node v i to the label which x l i is maximum:

s i = arg max l x l i (21)
An example of result is shown at Fig. 4.1.2.

Seeded segmentation

We now demonstrate the performance of power watershed with respect to the other seeded image segmentation algorithms. In the introduction we discussed how many of the leading graph-based segmentation algorithms (e.g., Grabcut, lazy snapping, closed-form matting) have graph cuts, random walker, shortest paths or watershed as an underlying component. Consequently, we will not compare the Power Watershed to any of the complete segmentation systems listed above, but rather against the comparable (component) algorithms of graph cuts, random walker, shortest paths and watershed. Additionally, to simplify the comparison we will not employ unary terms in our segmentations.

Quantitative assessment

Our experiments consist of testing five algorithms embodying different combinations of p and q, consisting of graph cuts, (GC), random walker (RW), shortest path (SP), watersheds/maximum spanning forest (MSF), and power watershed using the power q = 2. As before, we chose to employ the power watershed algorithm with q = 2 due to the uniqueness of the solution to (4) for this setting. 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. However, the seeds provided by the Grabcut database are generally equidistant from the ground truth boundary. To remove any bias from this seed placement on our comparative results, we produced an additional set of seeds by significantly eroding the original foreground seeds. The weights are set for all algorithms according to equation (1) with the value of β hand-optimized to provide the best results independently for each algorithm. As only the order of the weights is taken into account in the MSF and power watershed algorithms, those two algorithms are independant of β. We used the color gradient given by max((

R i -R j ) 2 , (G i -G j ) 2 , (B i -B j ) 2 )

TABLE 2

Mean errors computed between the segmentation masks and the ground truth images from the GrabCut database.

Symmetrically eroded ground truth are used as seeds. The weight parameter β was set to 600 for Graph cuts, 700 for Random walker, and 900 for Shortest paths in order to maximize the performances of each algorithms.

Tables 2 and 3 display the performance results for these algorithms. We quantify the error in the results using four different standard segmentation measures used in [90] , namely Boundary Error (BE), Rand Index (RI), Global Consistency Error (GCE), and Variation of Information (VoI). Good segmentation results are associ- When segmenting with the first seeding strategy (the seeds contained in the Grabcut database) , the shortest path algorithm is the best performer because this algorithm do well when the seeds are placed roughly equidistant from the desired boundary [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF] as they are with the first set of seeds.

The experiment on the second set of seeds shows that shortest paths is not robust to the seeds number and centering, it is with this set of seeds the worst performer. Graph cuts performs the best under this second seeding strategy but was the worst performers on the first one. Power watershed is in second position under the second seeding strategy, showing a good robustness to both seed quantity and location. It is interesting to note that with the first set of seeds, power watershed and maximum spanning forest results are quite similar, but with the asymetrically eroded seeds, the power watershed results outperform the standard maximum spanning forest (watershed) results. The second set of seeds contained many areas where several contours could possibly be found, given the seeds. The merging operation of the power watershed gathers undetermined areas and in performing the random walker in these ambiguous regions, often generates a better labeling than the arbitrary labeling produced by Prim's or Kruskal algorithms when computing the maximum spanning forest (watershed).

Computation time

For all MSF algorithms, including the power watershed algorithm, only the order of the weights in taken in consideration for the segmentation. Also, there is no parameter choice to make for β, and no exponential to take in the weight function, so it is possible to use a linear sort of the weights.

The worst-case complexity of the power watershed algorithm (obtained if all the edges weights are equal) is given by the cost of optimizing (4) for the given q. In best-case scenario (all weights have unique values), the power watershed algorithm has the same asymptotic complexity as the algorithm used for MSF computation, that is to say quasi-linear. In practical applications where the plateaus are have size less than some fixed value K, then the complexity of the power watershed algorithm matches the quasi-linear complexity of the standard watershed algorithm. In our experiments in Section 4 with practical image segmentation tasks, the dependence of the computation time on image size of the power watersheds is very similar to the dependence in standard watersheds. For generating the computation time for the graph cuts algorithm, we used the software provided at http://www.cs.ucl.ac.uk/staff/V.Kolmogorov/software.html and described in [18]. Our implementation of the shortest path is performed with a Fibonacci heap using double precision weights. For the implementation of Prim's algorithm, weights with integer precision were used, and red and black tree as a sorting data structure. Finally, the random walker algorithm was implemented following the multigrid method described in [40] for 2D image segmentation, and by a conjugate gradient descent method for 3D image segmentation. Unlike most watershed algorithms, the power watershed algorithm (with q = 2) has the property of providing a unique segmentation. Fig. 11 shows the behavior of the algorithm of our framework in presence of a plateau. Additionally, the power watershed (with q = 2) is not subject to the same shrinking bias exhibited by graph cuts segmentation. Fig. 12 compares the results of graph cuts and the power watershed on an example in which the shrinking bias could substantially affect the result. The power watershed is an MSF and therefore it inherits the standard properties of MSF, among others contrast invariance and scale invariance [START_REF] Allène | Some links between extremum spanning forests, watersheds and min-cuts[END_REF]. The contrast invariance property means that if a strictly monotonic transformation is applied to the weights of the graph, then the algorithm produces exactly the same result. This property is due to the fact that only the order or the weights is used to build a maximum spanning forest. The scale invariance property means that if we extend the image or graph in a way that does not change the relative ordering of weights, for example by linear interpolation, the result is invariant.

We summarize the performance of the algorithms of the framework:

• GC is a good fit for 2D image segmentation into two labels when the seeds are far away from the boundary (asymmetric seeding), but is too slow to be used for 3D segmentation. • SPF (geodesics) may be used if the object to segment is well centered around foreground and background seeds. • The RW is efficient and performs well for both seeding strategies (equidistant seeds and strongly asymmetric seeds). • Maximum Spanning Forest (watershed) algorithms provide better segmentations than SPF when seeds are not centered, and their fast computation time makes the algorithm suitable for 3D segmentation. • The power watershed algorithm when q = 2 has the additional property of a well-defined behavior in presence of plateaus improving also the quality of the segmentation compared to standard MSF. As an MSF, it is still sensitive to leaking, but less so than traditional algorithms due to the random walk behavior. The computational speed of the power watershed is faster than all of the algorithms except the pure MSF.

CONCLUSION

In this paper we clarified, simplified and extended the recent work connecting graph cuts and watershed [START_REF] Allène | Some links between min cuts, optimal spanning forests and watersheds[END_REF], [START_REF] Allène | Some links between extremum spanning forests, watersheds and min-cuts[END_REF]. Extending the framework of [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF], we have proposed a general framework encompassing graph cuts, random walker, shortest-path segmentation and watersheds. This connection allowed us to define a new family of optimal spanning forest for watershed segmentation algorithms using different exponents, which we termed the "power watershed". We produced an algorithm for computing the power watershed and our experiments showed that the power watershed with q = 2 retains the speed of the MSF algorithm while producing improved segmentations. In addition to providing a new image segmentation algorithm, this work also showed how unary terms could be employed with a standard watershed algorithm to improve segmentation performance. Viewed as energy minimization algorithms, graph cuts, random walker and shortest paths have found many different applications in the computer vision field that go beyond image segmentation, such as stereo correspondence, optical flow and image restoration (e.g., [82], [76], [53]). By placing the optimal spanning forest algorithm for watersheds in the same energy minimization framework as these other algorithms, watershed algorithms may find new uses and applications within the computer vision field beyond its traditional domain of image segmentation. Due to the relative speed of the optimal spanning forest algorithms, we believe that it may be an attractive alternative to current systems in these other applications of energy minimization.

Future work will develop along several directions. One direction is the further improvement of image segmentation algorithms using power watersheds as a component to larger systems in a similar manner as graph cuts, random walker and shortest paths have been used. Additionally, we hope to use the common framework for these algorithms to leverage existing ideas from the watershed literature into these other algorithms. In particular, hierarchical schemes [START_REF] Najman | Geodesic Saliency of Watershed Contours and Hierarchical Segmentation[END_REF], [START_REF] Najman | Ultrametric watersheds[END_REF], [START_REF] Najman | Ultrametric watersheds: a bijection theorem for hierarchical edge-segmentation[END_REF], [8], [START_REF] Guigues | Scale-sets image analysis[END_REF] looks like an interesting topic that can take advantage of the power watershed uniqueness. A second direction for future work will be to characterize the limits of the watershed algorithm as an energy minimization procedure [25]. Ultimately, we hope to employ power watersheds as a fast, effective alternative to the energy minimization algorithms that currently pervade the wide variety of applications in computer vision.

APPENDIX

Proof of Theorem

2 E p,q (x) w M p = e M |x i -x j | q + eij =e M w ij w M p |x i -x j | q . E p,q (x) w M p - E p,q (x * ) w M p = e M | xi -xj | q - e M |x * i -x * j | q + eij =e M w ij w M p | xi -xj | q - eij =e M w ij w M p |x * i -x * j | q . ( 22
)
The first part of ( 22) is bounded by 0,

e M | xi -xj | q - e M |x * i -x * j | q ≤ 0, (23)
because the energy obtained with x can not be greater than the one obtained by the optimal solution x * . More precisely, if there are no plateaus with different labels, the x * i , x * j computed on the edges e M with Algorithm 1 are equal, leading to a sum equal to 0. Else (if there are plateaus with different labels), e M |x ix j | q subject to the boundary constraints is minimized on the plateaus, so the solution is optimal.

The last part of ( 22) is also negative,

- eij =e M w ij w M p |x * i -x * j | q ≤ 0. ( 24
)
It only remains to bound the middle part of ( 22)

eij =e M w ij w M p | xi -xj | q ≤ eij =e M w ij w M p ≤ M 2 w M2 w M p , (25)
with M 2 the number of edges of weight inferior to w M , and w M2 the second maximum weight.

Thus we have

E p,q (x) w M p - E p,q (x * ) w M p ≤ M 2 w M2 w M p . (26)

Fig. 1 .

 1 Fig.1. Illustration of different steps in the proof of Thm. 1 for q = 2. The values on the nodes correspond to x, their color to s. The bold edges represents edges belonging to an MSF. (a) A weighted graph with two seeds, all maxima of the weight function are seeded, (b) First step, the edges of maximum weight are added to the forest, (c) After several steps, the next largest edge set belongs to a plateau connected to two labeled trees, (d) Minimize (4) on the subset (considering the merged nodes as a unique node) with q = 2 (i.e., solution of the combinatorial Dirichlet problem), (e) Another plateau connected to three labeled vertices is encountered, and (f) Final solutions x and s obtained after few more steps. The q-cut, which is also an MSF cut, is represented in dashed lines.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Example of behaviour of the power watershed algorithm for q = 2 with the formation of a plateau that was not present in the original graph. (a) Initialization: A weighted graph with two seeds. (b),(c),(d) First steps: The nodes of edges of maximum weight are merged , (e) The next largest edge set belongs to a plateau connected to two different labels, (f) Minimize (4) on the subset with q = 2 (i.e., utilize the random walker algorithm on the plateau), and (g) Final segmentation obtained after one more step.

Fig. 5 .

 5 Fig. 5. Example of the unseeded segmentation of a 3 × 3 image computed with a maximum spanning forest (watershed).

Fig. 6 .

 6 Fig. 6. Unseeded segmentation using unary terms. (a) Original image of blood cells, (b) graph cuts (c) Maximum Spanning Forest (watershed).

Fig. 7 .

 7 Fig. 7. Example segmentations with more than two labels. (a) Seeds, (b, c) power watershed result (q = 2).

Fig. 8 .

 8 Fig. 8. Example segmentations using the provided (top images) and skeletonized (bottom images) set of seeds on the Grabcut database images: (a) Seeds, (b) Graph cuts, (c) Random walker, (d) Shortest path, (e) Maximum spanning forest (standard watershed), and (f) Power watershed (q = 2).

  for a color image of red, green, and blue components R, G, B. The normalization is achieved by dividing the gradient by the maximum value of the gradient over every edge in the graph G. Example seeds and segmentations for the five algorithms with the first seeding strategy are shown at the top of Figure 8(a) and with the second seeding strategy at the bottom of Figure 8(a).

Fig. 11 .

 11 Fig. 11. (a) Image with foreground (red) and background (blue) seeds, (b) A segmentation obtained with graph cuts, (c) Segmentation obtained with Prim's algorithm for maximum spanning forest, and with a shortest path algorithm, (d) Segmentation obtained with random walker as well as power watershed with q = 2.

Fig. 9 .

 9 Fig. 9. Computation time for 2D and 3D seeded image segmentation. For each dimension, the times were generated by segmenting the same image scaled down.

Fig. 10 .

 10 Fig. 10. Example of 3D image segmentation. The foreground seed used for this image is a small rectangle in one slice of each lung, and the background seed is the frame of the image.

Fig. 12 .

 12 Fig. 12. Example comparison of Graph cuts and Power watershed faced to a weak foreground seeds quantity. (a) Seeds (the foreground seed is in red, indicated by an arrow), (b) Graph cuts segmentation result, (c) Power watershed (q = 2) result.

  images via lossy data compression. Computer Vision and Image Understanding, 110(2):212-225, May 2008. 11 Camille Couprie received in 2008 both the engineer degree at ESIEE Paris, majoring in computer science and graduating with highest honors, and the master degree in imaging science at the Universit é Paris Est, graduating with honors. Since Autumn 2008, she has been a PhD student supported by the French Direction G én érale de l'Armement MRIS program and the Centre National de la Recherche Scientifique. Her research focuses on image segmentation, filtering, combinatorial optimization techniques, PDEs and mathematical morphology. Other interests include stereovision, image registration, medical imaging and applied topology. She is a member of the IEEE. Leo Grady earned a BSc in Electrical Engineering at the University of Vermont in 1999 and a PhD in the Cognitive and Systems department at Boston University in 2003. Since Autumn 2003, he has worked as a Principal Research Scientist at Siemens Corporate Research (Princeton) in the Image Analytics and Informatics department. His research focuses on image segmentation, data clustering, biomedical imaging, complex networks, learning, compressive sensing and filtering using techniques from graph theory, discrete calculus, optimization and PDEs. Recent work has been devoted to writing his book "Discrete Calculus: Applied Analysis on Graphs for Computational Science", (Springer, 2010) which focuses on how the structure of connections between data points may be used to improve the analysis and understanding of the data. He is a member of the IEEE and the IEEE Computer Society. Laurent Najman received the Habilitation Diriger les Recherches in 2006 from University the University of Marne-la-Valle, the Ph.D. degree in applied mathematics from Universit é Paris-Dauphine in 1994 with the highest honor (Flicitations du Jury) and the engineering degree from the Ecole Nationale Sup érieure des Mines de Paris in 1991. He worked in the central research laboratories of Thomson-CSF for three years after his engineering degree, working on infrared image segmentation problems using mathematical morphology. He then joined a start-up company named Animation Science in 1995, as director of research and development. The particle systems technology for computer graphics and scientific visualization developed by the company under his technical leadership received several awards, including the "European Information Technology Prize 1997" awarded by the European Commission (Esprit programme) and by the European Council for Applied Science and Engineering as well as the "Hottest Products of the Year 1996" awarded by the Computer Graphics World journal. In 1998, he joined OC É Print Logic Technologies, as senior scientist. There he worked on various image analysis problems related to scanning and printing. In 2002, he joined the Informatics Department of ESIEE, Paris, where he is professor and a member of the Gaspard-Monge computer science research laboratory (LIGM), Universit é Paris-Est. His current research interest is discrete mathematical morphology. Hugues Talbot received the PhD degree in mathematical morphology from Ecole Nationale Sup érieure des Mines de Paris (ENSMP) with highest honors in 1993, under the guidance of professors Linn W. Hobbs (MIT), Jean Serra, and Dominique Jeulin (ENSMP). He was affiliated with CSIRO, Mathematical and Information Sciences, Sydney, Australia, between 1994 and 2004. He is currently a computer science professor at ESIEE, affiliated with Universit é Paris-Est, France. He has worked on a number of applied projects in image analysis with various companies, earning the Australian Institute of Engineers Award in 2004 and the DuPont Australian and New-Zealand Innovation Award in 2005 for work related to melanoma diagnosis. He has contributed more than 80 publications to international journals and conferences. His research interests include image segmentation, optimization and algorithms with application to image analysis and computer vision. He is a member of the IEEE and the IEEE Computer Society.

TABLE 3

 3 Mean errors computed between the segmentation masks and the ground truth images from the GrabCut database. Asymmetrically eroded ground truth are used as seeds.ated with low BE, high RI, low GCE and low VoI.

		BE	RI	GCE	VoI	Average
						rank
	Graph cuts	4.70	0.953	0.0380	0.284	1
	Power wshed	4.93	0.951	0.0407	0.297	2.5
	(q = 2)
	Random walker	5.12	0.950	0.0398	0.294	2.75
	MSF (Prim)	5.11	0.950	0.0408	0.298	3.5
	Shortest paths	5.33	0.947	0.0426	0.308	5

ACKNOWLEDGMENT

We would like to thank Jean Cousty for the numerous interesting and fruitful discussions.