Leo Grady

Laurent Najman

Hugues Talbot

Antonin Chambolle

Jean-Michel Morel

Thomas Pock

Gilles Bertrand

Jean Serra

Jean Cousty

Yukiko Ken- Mochi

Denis Bureau

Les Professeurs

John Chaussard

Benjamin Raynal

Roland Levillain

Olena Tankyevych

Yohan Thibault

Vincent Bismuth

Imen Melki

Vincent Nozick

Venceslas Biri

Nadine Dommanget

Anthony Giroud

Adrien Herubel

Keywords: variational problems, discrete calculus, graphs, regularization, convex optimization, image processing, segmentation, restoration, surface reconstruction, mathematical morphology traitement d'images, optimisation convexe, méthodes variationnelles, graphes, segmentation, restauration, vision par ordinateur, morphologie mathématique

I particularly want to thank him for captivating me with many

Chapter 1

Introduction

Computer vision encompasses a large diversity of important and fascinating problems that have revolutionized many endeavors in medicine, industrial vision, visualization, and multimedia. The computer vision community was able to contribute to these applications by formulating their corresponding problems in a solvable fashion. Particularly, energy minimization methods have been very successful in this regard.

Problems statement

Distinguishing objects from their background may sometimes appear like an easy and intuitive task for humans, but it is in fact complex and difficult to reproduce on a computer.

Effectively removing noise or reducing blur on images, although better defined, also remains an open problem. Example of such inverse problems are shown in Figure 1.1.

With respect to image segmentation, early optimization methods were formulated in terms of active contours and surfaces [START_REF] Kass | Snakes: Active contour models[END_REF] and then later level sets [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF]. These formulations were used to optimize energies of varied sophistication (e.g., using regional, texture, motion or contour terms [START_REF] Paragios | Geodesic Active Regions and Level Set Methods : Contributions and Applications in Artificial Vision[END_REF]) but generally converged to a local minimum, generating results that were sensitive to initial conditions and noise levels. Consequently, focus has shifted on energy formulations (and optimization algorithms) for which a global optimum can be found.

Many energy-based problems in vision involve the resolution of Partial Differential Equations (PDEs). Since the development of computer science, there has been a need to translate the classical mathematical models into a computable form to obtain numerical solution.

While classical solutions to PDEs involve either finite difference or finite elements schemes, in imaging data is available only in discrete form, typically as a regular graph. Designing formulations of PDE problems on graphs [START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF] has numerous advantages over previous schemes, such as removing the necessity to discretize the data again in a different form. In this way, it is possible to reuse existing efficient combinatorial optimization tools to minimize energies (network flows, shortest paths, minimum spanning trees) and to develop new ones. Finally, this approach generalizes these discrete algorithms to arbitrary graph, which is useful for modeling other kinds of data than planar images (e.g. triangulated or parametric surfaces, unevenly sampled images, etc).

→ ? → ?

→ ? → ? We begin our exposition by introducing our graph notations as well as a key operator used in this manuscript.

Discrete calculus

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

Discrete calculus has to be differentiated from traditional finite elements discretization.

Such discretization aims to produce approximate solution to problems defined in the continuum. The approximation becomes better as the discretization becomes finer and finer. In contrast, discrete calculus does not refer to the continuum. For example, as noted in [START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF],

social networks represent data that do not rely on a continuous domain. They are not a sampling of an underlying continuum. However, the tools of discrete calculus may be applied as we shall see in Chapter 2.

We introduce here the graph notations that are used in this manuscript.

A graph consists of a pair G = (V, E) with vertices v ∈ V and edges e ∈ E ⊆ V × V .

Let n represent the number of nodes, i.e. n = |V |, and m the number of edges of G, i.e. m = |E|. An edge, e, spanning two vertices, v i and v j , is denoted by e ij . In this thesis, we deal with weighted graphs that include weights on both edges and nodes. An edge weight is a value assigned to each edge e ij , and is denoted by w ij . We assume that w ij ∈ R * + and use w to denote the vector of R m containing the w ij for every edge e ij of G. In addition to edge weights, we may also assign weights to nodes. The weight of node v i is denoted by g i . In this work, we also assume that g i ∈ R * + . We use g to denote the vector of R n containing the g i for every node v i of G.

The incidence matrix of a graph is a key operator for defining combinatorial formulations of variational problems. Specifically, the incidence matrix A ∈ R n×m is known to define the discrete calculus analogue of the gradient, while A ⊤ is known to define the discrete calculus analogue of the divergence (see [START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF] and the references therein). The incidence matrix maps functions on nodes (a scalar field) to functions on edges (a vector field) and may be defined as As explained in [START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF], several important theorems such as the fundamental theorem of calculus or the Helmholtz decomposition hold in this discrete framework.

A e ij v k =          -1 if i = k, +1 if j = k,
In the next section, we briefly introduce some classical regularization methods used as references, as well as some of the optimization methods employed in this thesis.

Regularization methods

Numerous problems in signal processing and computer vision involve the optimization of an energy function composed of a regularization term and a data fidelity term. For example, image denoising problems may be solved by considering that the optimal recovered image is somehow regularized, that is to say that the variations (the local differences) of intensity are limited. For image segmentation problems, the sum of variations between neighboring labels may be penalized. In stereovision variational approaches, estimation of depth maps also assume a local consistency for the depth, and limit the sum of variations in order to obtain a piecewise smooth estimation.

In all these cases, the problem may be expressed as follows. A labeling x is estimated, x being the solution to the minimization of its total variations added to a term enforcing data fidelity. Remaining questions include, in particular "How do we translate the norm of the gradient of x numerically?", and "How to optimize the resulting the energy efficiently?". Also the assumption of local smoothness may be relaxed by considering of non-local gradient, for instance by extending graphs with long distance edges [START_REF] Buades | A non-local algorithm for image denoising[END_REF]. Before presenting the different regularization terms explored in this thesis, we shortly review two classical and efficient models for this matter.

Total variation

The Total Variation problem (TV) was introduced originally in computer vision by Shulman

and Hervé [START_REF] Shulman | Regularization of discontinuous flow fields[END_REF] as a regularizing criterion for optical flow computation and later Rudin, Osher

and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] as a regularizing criterion for image denoising. It has been shown to be quite efficient for smoothing images without affecting contours so much. Moreover, a major advantage of TV is that it is a convex problem, making it possible to find a global minimizer.

When applied to a 2D discrete image (x i,j ), 1 ≤ i, j ≤ N , the total variation minimization (x i+1,j -x i,j ) 2 + (x i,j+1 -x i,j ) 2 1 2 + D(x).

(1.3) using a finite-difference discretization popularized in [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. We note that other discrete formulation exist, for example a recent upwind variant of this model preserving more efficiently isotropic features was introduced in [START_REF] Chambolle | Upwind and multiscale finitedifference methods for total variationbased image smoothing[END_REF]. The Total variation problem has been expressed in weighted and non-local graphs in [START_REF] Bougleux | Discrete regularization on weighted graphs for image and mesh filtering[END_REF][START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF], leading to better penalization of the gradient.

A number of convex optimization techniques suitable for solving this problem have been employed for several decades. The most recent and efficient approaches are compared in [START_REF] Chambolle | An introduction to Total Variation for Image Analysis[END_REF].

Among the most efficient methods, one can cite Nesterov's algorithm [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF], Split-Bregman / Douglas-Rachford methods [START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF][START_REF] Goldstein | The split Bregman method for ℓ 1 -regularized problems[END_REF], and Chambolle-Pock's Primal-dual algorithm [54].

Most methods minimizing TV focus on image filtering as an application, and even if those methods are remarkably fast in denoising applications, we observe in Chapters 2 and 3 that segmentation problems require more iterations for those algorithms to converge.

Max flow -Min cut

The max-flow/min-cut problem on a graph is a classical problem in graph theory, for which the earliest solution algorithm goes back to Ford and Fulkerson [START_REF] Ford | Maximal flow through a network[END_REF]. Initial methods for global optimization of the boundary length of a region formulated the energy on a graph and relied on max-flow/min-cut methods for solution [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in N-D images[END_REF]. In the context of image segmentation, the seeded max-flow / min-cut model is also know as "Graph cuts". Graph cuts algorithms provide a mechanism for the discrete optimization of an energy functional [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF], which have been used in a variety of applications such as image segmentation, stereo, image stitching and texture synthesis (See Fig. 1.4). More specifically, the labeling x produced by graph cuts algorithms is an optimal (not necessarily unique) solution to min x∈{0,1} n e i,j ∈E w ij |x i -x j | + D(x). (1.4) This energy may be optimized using the non-polynomial but experimentally fast method of [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF]. As explained in [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF], more general -submodular -energies than (1.4) may be optimized using Graph cuts. Consequently, the graph cuts technique has been adapted to multilabel problems using various graph constructions. One notable construction is that of Ishikawa [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF] presenting the multilabel problem as a segmentation problem. The bold edges represent saturated edges, and constitute a minimal cut.

Combinatorial Dirichlet Problem

The Laplace equation arises in many physical situations : heat propagation, fluid dynamic, and electronics among others. The Dirichlet problem aims at finding a solution satisfying the Laplace equation subject to some boundary constraints.

In computer vision and image processing problems, solving the Dirichlet problem provide effective solution to labeling problems, given some markers. In this context, the combinatorial Dirichlet problem is written min

x∈R n e i,j ∈E w ij (x i -x j ) 2 + D(x).

(1.5)

In image processing, this problem has been applied for example to interpolation [START_REF] Grady | Anisotropic interpolation on graphs: The combinatorial dirichlet problem[END_REF],

inpainting [START_REF] Auroux | Fast algorithms for image processing and data assimilation[END_REF], image filtering [START_REF] Black | Robust anisotropic diffusion[END_REF], and seeded image segmentation [START_REF] Grady | Random walks for image segmentation[END_REF]. Figure 1.5 shows an application of the "random walker" [START_REF] Grady | Random walks for image segmentation[END_REF], an adaptation of Dirichlet problem to image segmentation. A summary of the different graph-based regularization terms previously introduced is given in Table 6.1, presenting as well the different novel regularization terms introduced in this thesis. Before the exposition of our contributions, we recall in the next section some combinatorial optimization methods useful in Chapter 4 and 5, as well as three algorithms for convex optimization that are used in Chapters 2 and 3.

Some optimization tools 1.3.1 Combinatorial methods

A good overview of combinatorial optimization techniques is presented in [START_REF] Papadimitriou | Combinatorial optimization, algorithms and complexity[END_REF]. We focus our exposition in this section on the combinatorial methods futher studied in Chapter 4.

Maximum Flow/Minimum cut algorithms

A transport graph G is as a graph with exactly one source s, with no entering edge, and one sink t, with no exiting edge.

A flow in a transport graph is a function associating a value to every edge according to some rules. The value of the flow in each edge must be less than the capacity (weight) of the edge, and the sum of the flow entering in each node -except s and t -must be equal to the sum of flow exiting the node. The max-flow problem consists of maximizing the amount of flow sent from the source toward the sink. This process creates a bottleneck for some edges.

If the flow of an edge is equal to its capacity, we say that the edge is saturated.

A cut is a partition of the nodes into two sets S and T , such that s is in S and t is in T . The capacity of a cut is the sum of the weights of all the edges crossing the cut, from the region S to the region T . Ford and Fulkerson [START_REF] Ford | Maximal flow through a network[END_REF] proved that the maximum flow is equals to the capacity of a minimal cut. As explained in Section 1.2.2, this result has been employed in various applications. Figure 1.4 shows image segmentation and stereovision results obtained using max-flow algorithms.

Most max-flow algorithms can be sorted in two categories, augmenting path [START_REF] Edmonds | Theoretical Improvements in Algorithmic Efficiency for Network Flow Problems[END_REF][START_REF] Ford | Maximal flow through a network[END_REF] and push-relabel methods [START_REF] Goldberg | A new approach to the maximum flow problem[END_REF][START_REF] Goldenberg | Fast geodesic active contours[END_REF]. Also min-cut algorithms without computing a maximum flow exist, a well known one being the Stoer-Wagner algorithm [START_REF] Stoer | A simple min-cut algorithm[END_REF]. Its complexity of O(nm+n 2 log(m)) is one of the most effective for general purpose min-cut/max-flow methods, but in practice the runtime on image data is not the fastest.

The regular lattice structure of images was exploited in the dedicated max-flow algorithm of [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF], resulting in faster runtimes.

Shortest Path methods

The shortest path problem consists of finding, from two nodes of a graph, a path of minimal cost (such that the sum of weights of component edges is minimized) in the graph. When the graph does not contain negative weights, Dijkstra's algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF] is very efficient. When implemented using Fibonacci heaps, its complexity is O(m + n log(n)). First introduced as a byproduct of reinitializing level set energies, the fast marching algorithm [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF] has been used

for shortest path computation [START_REF] Cohen | Minimal paths and fast marching methods for image analysis[END_REF] over sampled continuous domains. Shortest path methods are popular in image processing and particularly for thin object segmentation. The "Live wire", as known as "Intelligent scissors" [START_REF] Mortensen | Intelligent scissors for image composition[END_REF] consists for example of computing shortest paths between several points placed on the contour of objects to be extracted. Improvements have been proposed to limit the user interaction, and to define better weighting strategies [START_REF] Peyre | Geodesic Methods for Shape and Surface Processing[END_REF]. Other strategies for image segmentation using shortest path based algorithms include the construction of shortest path forests. Given foreground and background seeds, each pixel is assigned to the foreground label if there is a shorter path from that pixel to a foreground seed than to any background seed. 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 [6,[START_REF] Criminisi | GeoS: Geodesic image segmentation[END_REF][START_REF] Falcão | The image foresting transformation[END_REF].

Maximum Spanning Trees and Forests

A spanning tree of a graph is a tree (a connected graph without cycles) connecting all the nodes of the graph. The weight of a spanning tree is defined by the sum of edge weights composing the tree. There exists several greedy algorithms minimizing or maximizing this weight, called respectively minimum or maximum spanning tree algorithms. Borůvka [START_REF] Borůvka | O jistém problému minimálním[END_REF] was the first to propose an algorithm for this problem. The Prim [START_REF] Prim | Shortest connection networks and some generalizations[END_REF], and Kruskal [START_REF]On the shortest spanning tree of a graph and the traveling salesman problem[END_REF] algorithms are the most commonly used. When using a union-find data structure [START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF] for cycles detection, the complexity of maximum spanning tree algorithms is quasi-linear [START_REF] Chazelle | A minimum spanning tree algorithm with inverse-ackermann type complexity[END_REF]. We shall see in Chapter 4 that these greedy procedures, which are very fast, may be used for optimizing meaningful and useful energies. The following section introduces a few convex optimization methods employed in Chapters 2 and 3.

Interior Point Method

Interior point methods have been successfully applied to large size problems, and in particular in image processing [START_REF] Bhusnurmath | Graph cuts via ℓ 1 norm minimization[END_REF][START_REF] Chan | A nonlinear primal-dual method for total variation based image restoration[END_REF][START_REF] Kim | An efficient method for compressed sensing[END_REF].

As detailed in [START_REF] Boyd | Convex Optimization[END_REF], interior point methods solve problems of the form min

x f 0 (x), s. t. f i (x) ≤ 0 i = 1, ..., M, Hx = b, (1.6)
where x is a vector ∈ R N , H a matrix of R P ×N such that rank(H) = P < N , b a vector ∈ R P , and f 0 , ..., f M are convex, twice differentiable functions of R N → R.

In Chapter 2 we shall be interested in the dual solution of the considered optimization problem, this is why we focus our exposition on a primal-dual interior point method.

The dual of Problem (1.6) can be computed by mean of a Lagrangian function. The basic idea is to take the constraints into account by augmenting the objective function with a weighted sum of the constraint functions.

The Lagrangian associated with the problem (1.6) is defined by

L(x, λ, ν) = f 0 (x) + M i=1 λ i f i (x) + M i=1 ν i (H i x -b i ),
with λ ∈ R M and ν ∈ R M two Lagrange multipliers. Then the dual objective function is defined by

h(λ, ν) = inf x L(x, λ, ν).
Taking the best lower bound of the Lagrangian function leads to the dual problem max λ,ν h(λ, ν),

s. t. λ ≥ 0.
The primal-dual interior point (PDIP) algorithm computes iteratively the primal x and dual λ, ν variables so that the Karush-Kuhn-Tucker (KKT) optimality conditions are satisfied. The KKT conditions guarantee that if the functions f i are convex, and x, λ and ν are satisfying

f i (x) ≤ 0, i = 1, ..., M Hx = b λ i ≥ 0 i = 1, ..., M λ i f i (x) = 0 i = 1, ..., M ∇f 0 (x) + M i=1 λ i ∇f i (x) + H ⊤ ν = 0, (1.7)
then x and (λ, ν) are primal-dual optimal.

The PDIP algorithm solves the problem (1.6) by applying the Newton method to a sequence of a slightly modified version of KKT conditions. We define the vector f (x) of R M as

f (x) = [f 1 (x), ..., f M (x)]. Algorithm 1: Primal-dual interior point algorithm Input x that satisfies f 1 (x) < 0, ..., f M (x) < 0, λ > 0, µ > 1, ǫ > 0 Output x such that f 0 (x) is minimized under the constraints of (1.6) Repeat       1.
Compute the surrogate duality gap η = -f (x) ⊤ λ, and set t = µM/ν. The second step of the PDIP algorithm is called the Newton step and consists of solving

the M ∆ y = r system    ∇ 2 f 0 (x) + N i=1 λ i ∇ 2 f i (x) Df (x) ⊤ H ⊤ -diag(λ)Df (x) -diag(f (x)) 0 H 0 0       ∆ x ∆ λ ∆ ν    = -    r d = ∇f 0 (x) + Df (x) ⊤ λ + H ⊤ ν r c = -diag(λ)f (x) -(1/t)1 r p = Hx -b    (1.8)
with r d , r c , and r p the dual, central, and primal residuals, and

Df (x) =     ∇f 1 (x) ⊤ . . . ∇f M (x) ⊤     .
Although the complexity of the given primal-dual interior point method is not studied in [START_REF] Boyd | Convex Optimization[END_REF], the authors observe that the growth of the number of iterations in function of the problem size is approximately logarithmic for a linear programming problem. We note that for solving the nonlinear problem addressed in Chapter 2, the number of iterations with respect to the problem size is observed to be linear.

Some proximal algorithms

Proximal methods provide efficient solutions to convex optimization problems where the objective function can be split into a sum of convex functions. This body of methods are applicable to large size problems as well as to non smooth or non-finite functions.

Let us denote Γ 0 (R N ) the class of convex functions f of R N to ]-∞, +∞] that are proper (such that dom f = z ∈ R N f (z) < +∞ = ∅) and lower semi-continuous.

Specifically, proximal methods may be employed to solve problems of the form min

x∈R N M i=1 f i (L i x) (1.9) given M linear operators L i ∈ R Q i ×N , i ∈ {1, . . . , M }, where M, (∀i ∈ {1, . . . , M })Q i , and
N are natural positive integers and

f i ∈ Γ 0 (R Q i ).
The proximal operator of a function is a natural extension of the notion of projection onto a convex set, and is defined as follows.

Definition 1.3.1. For every u ∈ R N , the function v → ϕ(v)+ u-v 2 /2, where ϕ ∈ Γ 0 (R N ),
achieves its minimum at a unique point denoted by prox ϕ u, i.e.,

(∀u ∈ R N ), prox ϕ u = arg min v∈R N ϕ(v) + 1 2 u -v 2 . ( 1
ı C : x →    0 if x ∈ C, +∞ otherwise, (1.11)
and the projection

P C x of x ∈ R N onto C is the solution y to min y∈R N ı C (y) + 1 2 x -y 2 .
(1.12)

Among classical existing proximal methods, one can cite the augmented Lagrangian methods, Simultaneous Direction Method of Multipliers, and Douglas Rachford algorithms [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF] among others. More recently parallel proximal methods have been introduced in order to more easily take advantage of parallel architectures [START_REF] Arias | A monotone+skew splitting model for composite monotone inclusions in duality[END_REF][START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF].

1.3.3.1 Parallel Proximal Algorithm Algorithm 2: Parallel ProXimal Algorithm (PPXA) optimizing (1.13) γ ∈ ]0, +∞[ , λ ∈]0, 2[(relaxation parameter) (ω 1 , . . . , ω M ) ∈]0, 1[ M , M i=1 ω i = 1, k = 1 x 1 = M i=1 ω i t i,1 Repeat until convergence                For (in parallel) i = 1, . . . , M p i,k = prox γ ω i f i t i,k c k = M i=1 ω i p i,k For (in parallel) i = 1, . . . , M t i,k+1 = t i,k + λ 2c k -x k -p i,k x k+1 = x k + λ(c k -x k ) k = k + 1
The Parallel ProXimal Algorithm, also called PPXA, has recently been introduced by Combettes and Pesquet [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]. It computes a sequence of (x k ) k∈N * converging toward the minimizer of min

x∈R N M i=1 f i (x) (1.13) for f i ∈ Γ 0 (R M i ), i ∈ {1, . . . , M }.
The algorithm, shown in Algorithm 2 computes in parallel the proximal operators of the functions f i of (1.13), merges them, and updates the iterates also in parallel. Examples of applications, and comparison with existing proximal methods can be found in [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Pustelnik | Proximal methods for the resolution of inverse problems : application to positron emission tomography[END_REF]. A generalization of the algorithm optimizing (1.9), and also including the simultaneous direction of multipliers method is present in [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF]. A simple example illustrating the behavior of PPXA is given in Figure 1.8.

Example: find a point x ∞ in set C at the closest distance from D.

min x∈R N ι C (x) f 1 (x) + min v∈R N ι D (v) + 1 2 ||x -v|| 2 f 2 (x) λ = γ = k = 1, ω 1 = ω 2 = 1 2 , t 1,1 = t 2,1 = x 1 . p 1,k = P C t 1,k p 2,k = 1 2 (t 2,k + P D t 2,k ) c k = 1 2 (p 1,k + p 2,k ) t 1,k+1 = t 1,k + p 2,k -x k t 2,k+1 = t 2,k + p 1,k -x k x k+1 = c k prox f 1 t 1 = P C t 1 prox f 2 t 2 = 1 2 (t 2 + P D t 2 )

M+SFBF algorithm

The "Monotone + Skew Forward Backward Forward" (M+SFBF) algorithm of Briceño and

Combettes [START_REF] Arias | A monotone+skew splitting model for composite monotone inclusions in duality[END_REF][Prop. 4.4] may be employed to perform the joint optimization of the two following problems min

x∈R N M i=1 ω i f i (L i x), (1.14)
given M with i ∈ {1, . . . , M } linear operators .15) where for all i ∈ {1, . . . , M }, f * i denotes the conjugate function of f i , whose definition is recalled below.

L i ∈ R Q i ×N , f i ∈ Γ 0 (R Q i ), ω i ∈]0, 1] such that M i=1 ω i = 1, and min v 1 ∈R Q 1 ,...,v M ∈R Q M i=0 M ω i L * i v i =0 M i=1 ω i f * i (v i ). ( 1
Definition 1.3.2. Let f : R N → R. The conjugate function f * : R N → R of f is defined as f * (x) = sup y∈dom f (x ⊤ y -f (y)). Algorithm 3: M+S Forward Backward Forward algorithm Initialization     For every i ∈ {1, . . . , M }, ω i ∈]0, 1] such that M i=1 ω i = 1 For every i ∈ {1, . . . , M }, x i ∈ R N , v i ∈ R Q i β = max 1≤i≤M L i , let ǫ ∈]0, 1/(β + 1)[, let in [ǫ, (1 -ǫ)/β] Repeat until convergence                      x = M i=1 ω i x i For i = 1, . . . , M y 1,i = x i -γL * i v i y 2,i = v i + γL i x i p 1 = M i=1 ω i y 1,i For i = 1, . . . , M         p 2,i = prox γf * i y 2,i q 1,i = p 1 -γL * i p 2,i q 2,i = p 2 + γL i p 1 x i = x i -y 1,i + q 1,i v i = v i -y 2,i + q 2,i
The weak convergence of the sequence of iterates x toward a minimizer of M i=1 ω i f i (L i ) is established (and further generalized) in [START_REF] Arias | A monotone+skew splitting model for composite monotone inclusions in duality[END_REF] under the assumption that

∃x ∈ R N | i ω i L * i • (δf i ) • L i x = 0.
An advantage of the M+SFBF algorithm is the absence of matrix inversion for the linear operator handling, contrary to the generalized PPXA+ method [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF]. We discuss in Chapter 3 cases where the M+SFBF algorithm is more efficient than PPXA in an application to image restoration.

Manuscript organization

The manuscript is organized as follows. Chapter 2

Combinatorial Continuous Maximum Flow

Maximum flow (and minimum cut) algorithms have had a strong impact on computer vision.

In particular, as briefly reviewed in Chapter 1, graph cut algorithms provide a mechanism for the discrete optimization of an energy functional which has been used in a variety of applica- discrete formulation is different from the classical max-flow problem. We then apply an appropriate combinatorial optimization technique to this combinatorial continuous max-flow (CCMF) problem to find a null-divergence solution that exhibits no metrication artifacts and may be solved exactly by a fast, efficient algorithm with provable convergence. Finally, by exhibiting the dual problem of our CCMF formulation, we clarify the fact, already proved by Nozawa in the continuous setting, that the max-flow and the total variation problems are not always equivalent.

Introduction

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 chapter, we focus on image segmentation as the example application on which to base our exposition. Indeed, segmentation has played a prominent (and early) role in the development of the use of global optimization methods in computer vision, and often forms the basis of many other applications [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF][START_REF] Darbon | Image restoration with discrete constrained total variation part II: Levelable functions, convex priors and non-convex cases[END_REF][START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF].

Metrication errors from the min-cut approach

Briefly described in Section 1.2.2, the max-flow/min-cut problem on a graph is a classical problem in graph theory. It was soon realized that these methods introduced a metrication error. Metrication errors are clearly visible when gradient information is weak or missing, for example in the case of medical image segmentation or in some materials science applications like electron tomography, where weak edges are unavoidable features of the imaging modality.

Metrication errors are also far more obvious in 3D than in 2D and increasingly so, as the dimensionality increases (see for example the 3D-lungs example in [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF]). Furthermore, metrication artifacts are even more present in other applications such as surface reconstruction and rendering, where the artifacts are a lot worse than in image segmentation. Various solutions for metrication errors were proposed. One solution involved the use of a highly connected lattice with a specialty edge weighting [START_REF] Yuri | Computing geodesics and minimal surfaces via graph cuts[END_REF], but the large number of graph edges required to implement this solution could cause memory concerns when implemented for large 2D or 3D images [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF].

Continuous approaches

To avoid the metrication artifacts without increasing memory usage, one trend in recent years has been to pursue spatially continuous formulations of the boundary length and the related problem of total variation [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF][START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Nikolova | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF]. Historically, a continuous max-flow (and dual min-cut problem) was formulated by Iri [START_REF] Iri | Theory of flows in continua as approximation to flows in networks[END_REF] and then Strang [START_REF] Strang | Maximum flows through a domain[END_REF]. Strang's continuous formulation provided an example of a continuization (as opposed to discretization) of a classically discrete problem, but was not associated to any algorithm. Work by Appleton and Talbot [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF] provided the first PDE-based algorithm to find Strang's continuous max-flows and therefore optimal min-cuts. This same algorithm was also derived by Unger et al. from the standpoint of minimizing continuous total variation [START_REF] Unger | TVSeg -Interactive total variation based image segmentation[END_REF]. Adapted to image segmentation, this algorithm is shown to be equivalent [START_REF] Unger | Interactive globally optimal image segmentation[END_REF] to the Appleton-Talbot algorithm and has been demonstrated to be fast when implemented on massively parallel architectures. Unfor- Links and differences with total variation minimization G. Strang [START_REF] Strang | Maximum flows and minimum cuts in the plane[END_REF] has shown the continuous max-flow problem for the l 2 norm to be the dual of the total variation (TV) minimization problem under some assumptions. Most methods minimizing TV focus on image filtering as an application, and even if those methods are remarkably fast in denoising applications, segmentation problems require a lot more iterations for those algorithms to converge. As stated previously, continuous max-flow is dual with total variation in the continuous setting under restrictive regularity conditions [START_REF] Nozawa | Max-flow min-cut theorem in an anisotropic network[END_REF]. In fact, Nozawa [START_REF] Nozawa | Examples of max-flow and min-cut problems with duality gaps in continuous networks[END_REF] showed that there is a duality gap between weighted TV and weighted max-flow under some conditions in the continuous domain. Thus, it is important to note that weighted TV problems are not equivalent to the weighted maximum flow problem studied in this chapter. Furthermore, many works assume that the continuous-domain duality holds algorithmically, but we show later in this chapter that at least in the combinatorial case this is not true.

The previous works for solving the max-flow problem illustrate two difficulties with continuous-based formulation, that are (1) the discretization step, which is necessary for deriving algorithms, but may break continuous-domain properties; and (2) the convergence, both of the underlying continuous formulation, and the associated algorithm, which itself depends on the discretization. It is very well known that even moderately complex systems of PDEs may not converge, and even existence of solutions is sometimes not obvious [START_REF] Fefferman | Existence and smoothness of the Navier-Stokes equation[END_REF].

Even when existence and convergence proofs both exist, sometimes algorithmic convergence may be slow in practice. For these reasons, combinatorial approaches to the maximal flow and related problems are beginning to emerge.

Combinatorial approaches

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

Motivation and contributions

In this chapter, we pursue a combinatorial reformulation of the continuous max-flow problem initially formulated by Strang, which we term combinatorial continuous maximum flow (CCMF). Viewing our contribution differently, we adopt a discretization of continuous 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 max-flow problem produces a (divergence-free) flow problem defined on an arbitrary graph. Strikingly, CCMF are not equivalent to the discretization of continuous max-flows produced by Appleton and Talbot or that of Lippert (if for no other reason than the fact that CCMF is defined on an arbitrary graph). Moreover, CCMF is not equivalent to the classical max-flow on a graph. In particular, we will see that the difference lies in the fact that capacity constraints in classical max-flow restrict flows along graph edges while the CCMF capacity constraints restrict the total flow passing through the nodes.

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

The CCMF problem has several desirable properties, including:

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

2. in contrast to continuous max-flow algorithms of Appleton-Talbot (AT-CMF) and equivalent, the solution to the CCMF problem can be computed with guaranteed convergence in practical time;

3. the CCMF problem is formulated on an arbitrary graph, which enables the algorithm to incorporate nonlocal edges [START_REF] Buades | A non-local algorithm for image denoising[END_REF][START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF][START_REF] Grady | Faster graph-theoretic image processing via small-world and quadtree topologies[END_REF] or to apply it to arbitrary clustering problems defined on a graph; and 4. the algorithm for solving the CCMF problem is fast, easy to implement, compatible with multi-resolution grids and is straightforward to parallelize.

Our computation of the CCMF dual further reveals that duality between total variation minimization and maximum flow does not hold for CCMF and combinatorial total variation (CTV). The comparison between those two combinatorial problems is motivated by several interests:

1. for clarification: in the continuous domain, the duality between TV and MF holds under some regularity constraints. In the discrete anisotropic domain (i.e. Graph Cuts), this duality always holds. However, In the isotropic weighted discrete case, i.e. in the CTV, AT-CMF or CCMF cases, we are not aware of any discretization such that the duality holds. It is not obvious to realize that the duality does not hold in the discrete setting even though it may in the continuous case. We present clearly the differences between the two problems, theoretically, and in term of results;

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

3. for efficiency: in image segmentation, the fastest known algorithms to optimize CTV are significantly slower than CCMF. Therefore, there is reason to believe that CCMF may be used to efficiently optimize energies for which TV has been previously shown to be useful (e.g., image denoising).

In the next section, we review the formulation of continuous max-flow, derive the CCMF formulation and its dual and then provide details of the fast and provably convergent algorithm used to solve the new CCMF problem.

Method

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

The continuous max-flow (CMF) problem

First introduced by Iri [START_REF] Iri | Theory of flows in continua as approximation to flows in networks[END_REF], Strang presents in [START_REF] Strang | Maximum flows through a domain[END_REF] an expression of the continuous maxi-

mum flow problem max F st , s.t. ∇ • - → F = 0, || - → F || ≤ g. (2.1)
Here we denote by F st the total amount of flow going from the source location s to the sink location t, -→ F is the flow, and g is a scalar field representing the local metric distortion.

The solution to this problem is the exact solution of the geodesic active contour (GAC) (or surface) formulation [START_REF] Caselles | Geodesic active contours[END_REF]. In order to solve the problem (2.1), the Appleton-Talbot algorithm (AT-CMF) [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF] solves the following partial differential equation system:

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

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

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

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

A discrete calculus formulation

Before establishing the discrete calculus formulation of the continuous max-flow problem, we specify our notation. We refer to Section1.1 for general definitions on graphs. In this chapter, a transport graph G(V, E) comprises two additional nodes, a source s and a sink t, and additional edges linking some nodes to the sink and some to the source. Including the source and the sink nodes, the cardinalities of G are given by n = |V | and m = |E|. We define a flow through edge e ij as F ij where F ij ∈ R and use the vector F ∈ R m to denote the flows through all edges in the graph . Each edge is assumed to be oriented, such that a positive flow on edge e ij indicates the direction of flow from v i to v j , while a negative flow indicates the direction of flow from v j to v i .

The incidence matrix is a key operator for defining a combinatorial formulation of the continuous max-flow problem. In our formulation of continuous max-flow, we use the expression |A| to denote the matrix formed by taking the absolute value of each entry individually.

Given these definitions, we now produce a discrete (combinatorial) version of the continuous max-flow of (2.1) on a transport graph. As in [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF][START_REF] Grady | The piecewise smooth Mumford-Shah functional on an arbitrary graph[END_REF][START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF], the continuous vector field indicating flows may be represented by a vector on the edge set, F . Additionally, the combinatorial divergence operator allows us to write the first constraint in (2.1) as A ⊤ F = 0.

The second constraint in (2.1) involves the comparison of the point-wise norm of a vector field with a scalar field. Therefore, we can follow [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF][START_REF] Grady | The piecewise smooth Mumford-Shah functional on an arbitrary graph[END_REF][START_REF] Grady | Discrete Calculus: Applied Analysis on Graphs for Computational Science[END_REF] to define the point-wise ℓ 2 norm of the flow field F as |A ⊤ |F 2 . In our notations here, as in the rest of the chapter, F 2 = F • F denotes a element-wise product, "•" denoting the Hadamard (element-wise) product between the vectors, and the square root of a vector is also here and in the rest of the chapter the vector composed of the square roots of every elements.

Putting these pieces together, we obtain max F st ,

s.t. A ⊤ F = 0, |A ⊤ |F 2 ≤ g 2 . (2.3) Written differently, the capacity constraint is ∀v i ∈ V, e ij ∈E F 2 ij + e ji ∈E F 2 ji ≤ g 2 i .
Compare this formulation to the classical max-flow problem on a graph, given in our notation as max F st ,

s.t. A ⊤ F = 0, |F | ≤ w. (2.4)
By comparing these formulations of the traditional max-flow with our combinatorial formulation of the continuous max-flow, it is apparent that the key difference is in the capacity constraints. In both formulations, the flow is defined through edges, but in the classical case the capacity constraint restricts flow through an edge while the CCMF formulation restricts the amount of flow passing through a node by taking in account the neighboring flow values. This contrast-weighting applied to nodes (pixels) has been a feature of several algorithms with a continuous formulation, including geodesic active contours [START_REF] Caselles | Geodesic active contours[END_REF], CMF [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF], TVSeg [START_REF] Unger | TVSeg -Interactive total variation based image segmentation[END_REF], and [54]. On the other hand, the problem defined in (2.3) 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 max-flow formulation and therefore our graph formulation of Strang's formulation carries over these same capacities. The most important point here is not having expressed the capacity on the nodes of the graph but rather how the flow is enforced to be bounded by the metric in each node. Bounding the norm of neighboring flow values in each node simulates a closer behavior of the continuous setting than if no contribution of neighboring flow was made as in the standard max-flow problem. Figure 2.2 illustrates the relationship of the edge capacity constraints in the classical max-flow problem to the node capacity constraints in the CCMF formulation. In particular, the minimal cut is localized on saturated edges in the classic max flow case, whereas the minimal cut appears on saturated nodes in the CCMF case. More formally,

Definition 2.2.1. A node v i is saturated if |A ⊤ | i F 2 = g 2 i , where |A ⊤ | i indicates the i th row of |A ⊤ |.
In the context of image segmentation, the vector g varies inversely to the image gradient.

We propose to use, as in [START_REF] Unger | TVSeg -Interactive total variation based image segmentation[END_REF],

g = exp(-β||∇I|| 2 ), (2.5)
where I indicates the image intensity. For simplicity, this weighting function is defined for greyscale images, but g could be used to penalize changes in other relevant image quantities, such as color or texture. Before addressing the solution of the CCMF problem, we consider the dual of the CCMF problem and, in particular, the sense in which it represents a minimum cut. Since the cut weights are present on the nodes rather than the edges, we must expect the minimum cut formulation to be different from the classical minimum cut on a graph. The difference between classical max-flow and CCMF is that classical max-flow uses edge-weighted capacities while CCMF uses node-weighted capacities. This difference is manifest in the different solutions obtained for both algorithms and the algorithms required to find a solution. Specifically, the solution to the CCMF problem on a lattice does not exhibit metrication bias.

The CCMF dual problem

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

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

Property 1. In a transport graph G with m edges, n nodes, we define a vector c of R m , composed of zeros except for the element corresponding to the source/sink edge which is -1.

Let λ and ν be two vectors of R n . The CCMF problem max

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

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

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

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

The proof, using the Lagrange duality theorem, is given in Appendix A.

(ν i -ν j ) 2 λ i + λ j + source/sink enforcement 1 4 (ν s -ν t -1) 2 λ s + λ t s. t. λ i ≥ 0 ∀i ∈ V.
(2.10)

Interpretation:
The optimal value λ * is a weighted indicator of the saturated vertices (See Def. 2.2.1):

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

ν * (v i ) =      0 + δ if v i ∈ S,
a number between (0 + δ) and

(1 + δ) if |A ⊤ | i F 2 = g(v i ) 2 , 1 + δ if v i ∈ T.
The expression (2.8) of the CCMF dual shows that the problem is equivalent to finding a minimum weighted cut defined on the nodes.

Finally, the "weighted cut" is recovered in (2.10), 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.3.

Several optimization methods have been tested for solving the CCMF problem, its dual and both jointly. We report in Appendix A.1.3 a proximal method for solving the primal The dual problem to CCMF is a node-weighted minimum cut in which the variable λ is a weighted indicator vector labeling boundary nodes and the variable ν is a nearly binary vector indicating the source/sink regions. As a result, the contours of ν are slightly blurry. This is due to the equilibrium effect between the two dual variables. In practice, as λ is nonzero only in presence of a contour, ν is binary almost everywhere, except on a very thin line.

problem and a first order method for solving the dual problem. The most efficient method we found, solving both problems jointly, is the primal dual interior point method detailed in the next section.

Primal Dual Interior Point Method

When considering how to solve the CCMF problem (2.3), the first key observation is that the constraints bound a convex feasible region.

The real valued functions f i : R m → R defined for every node i = 1, . . . , n as f i (F ) =

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

R n as f (F ) = [f 1 (F ), . . . , f n (F )] ⊤ .
Since the constraints are convex, the CCMF problem may be solved with a fast primal dual interior point method (see [START_REF] Boyd | Convex Optimization[END_REF]), which we now review in the specific context of CCMF. The primal dual interior point (PDIP) algorithm iteratively computes the primal F and dual λ, ν variables so that the Karush-Kuhn-Tucker (KKT) optimality conditions are satisfied. This algorithm solves the CCMF problem (2.6) 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 4. Specifically for CCMF, the system M ∆ y = r system may be written

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

∇f i (F ) = 2|A ⊤ | i • F, Df (F ) =     ∇f 1 (F ) ⊤ . . . ∇f n (F ) ⊤     , ∇ 2 f i (F ) = diag     2|A ⊤ | i1 . . . 2|A ⊤ | im     ,
with "•" denoting the Hadamard (element-wise) product between the vectors.

Consequently, the primary computational burden in the CCMF algorithm is the linear system resolution required by (2.12). In the result section we present execution times obtained in our experiments using a conventional CPU implementation.

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

Comparison between CCMF and existing related approaches 2.3.1 Min cut

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

u gT |Au + c|. (2.13)
We note that the ℓ 1 norm of the gradient of the solution u of the above min-cut expression turns into an ℓ 2 norm of the gradient of the solution ν in the CCMF dual (2.7).

Primal Dual Total variations

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

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

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

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 (in 2D)

min u Ω g ||∇u|| dx dy, s. t. Γ uf ds = 1, (2.16)
where Ω represents the image domain and Γ the boundary of Ω. The source and sink membership is represented by f , such that f (x, y) > 0 if (x, y) belongs to the sink, f (x, y) < 0 if (x, y) belongs to the source, and f = 0 otherwise. Considering a transport graph G, and a vector u defined on the nodes, this continuous problem may be written with combinatorial operators

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

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

u v i ∈V g i e ij ∈E (u i -u j ) 2 , s. t. u s -u t = 1.
(2.18)

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

CCMF and CTV are not dual

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

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

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

Transport graph with weights g on the nodes

Combinatorial Continuous Max-Flow (CCMF)

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

min u g ⊤ |A ⊤ |(Au) 2 , s. t. u s -u t =
1 0 0 1 1 0 0 1 0 0 0.2 F * st = 1.63 g ⊤ |A ⊤ |(Au * ) 2 = 8.16
Table 2.1: Example illustrating the difference between optimal solutions of combinatorial continuous maximum flow problem (CCMF), and combinatorial total variation (CTV).

Even if CCMF is not dual with CTV, the two problems are weakly dual (more details are given in Appendix A.1.2).

Results

We now present applications of Combinatorial Continuous Maximum Flow in image segmentation. To be used as a segmentation algorithm, three solutions are possible: the dual variable ν may be used directly if the user needs a matted result, otherwise ν may be thresholded, or finally an isoline or isosurface may be extracted from ν. In the introduction we discussed how several works are related to CMF. Some are equivalent or slight modifications of AT-CMF for non-segmentation applications [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF][START_REF] Unger | Interactive globally optimal image segmentation[END_REF][START_REF] Unger | TVSeg -Interactive total variation based image segmentation[END_REF][START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF].

(a) (b) (c) (d) GC CCMF GC CCMF GC CCMF
As in the original AT-CMF work, none of these come with a convergence proof. In contrast, the work of [START_REF] Pock | An algorithm for minimizing the Mumford-Shah functional[END_REF] and [START_REF] Lippert | Discrete approximations to continuum optimal flow problems[END_REF] are provably convergent but both are very slow, and do not generalize easily to 3D image segmentation. Consequently, it seems reasonable to compare our segmentation method to the original AT-CMF as representative of the continuous approach, as well as graph cuts, which represent the purely discrete case. Finally, even if the total variation problem is different from the CMF problem (existence of a duality gap in the continuous [START_REF] Nozawa | Examples of max-flow and min-cut problems with duality gaps in continuous networks[END_REF]), the two problems are related enough to be compared. Specifically, we will The difference between those results is due to the fact that the CMF algorithm enforces exactly the source and sink points, leading to discretization around the disks. In contrast, the boundary localized around the seeds of ν is smooth, composed of grey levels. Thus the resulting isosurface computed by CCMF is more precise.

Metrication artifacts and minimal surfaces

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

In the continuous setting, the maximum flow computed in a 3D volume produces a minimal surface. The CCMF formulation may be also recognized as a minimal surface problem. In the dual formulation, the objective function is equivalent to a weighted sum of surface nodes.

In [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF], Appleton and Talbot compared the surfaces obtained from their algorithm with the analytic solution of the catenoid problem to demonstrate that their algorithm was a good approximation of the continuous minimal surface and was not creating discretization artifacts.

The catenoid problem arises from consideration of two circles with equal radius whose centers lie along the z axis. The minimal surface which forms to connect the two circles is known as a catenoid. The catenoid appears in nature, for example by creating a soap bubble between two rings. In order to demonstrate that CCMF is also finding a minimal surface, we performed the same catenoid experiment as was in [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF]. The results are displayed in Figure 2.5, where we show that CCMF approximates the analytical solution of the catenoid with even greater fidelity than the AT-CMF example. In order to verify that the solution obtained by CCMF approximates perimeters of planar objects in 2D, we segmented several shapes -squares, discs -with known perimeters, scaled at different size, and compared the analytic perimeters with the energies obtained by CCMF.

The results are reported in Fig. 2.6, where we observe that the true boundary length is proportional to the energy obtained with CCMF.

Stability, convergence and speed

We may compare the segmentation results using ν to Appleton-Talbot's result using P . We recall that AT-CMF solves the partial differential equation system (2.3) in order to solve the continuous maximum flow problem (2.1), but no proof was given for convergence. The potential function P approximates an indicator function, with 0 values for the background labels, and 1 for the foreground. It can be difficult to know when to stop the AT-CMF algorithm, since its iterations may oscillate, as displayed in Figure 2.7 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 [54] optimizing a total variation based energy using a dual variable. In comparison to our work, the dual variable is not a flow (no constraint on the divergence and no proof of convergence toward a minimal surface as defined in [START_REF] Strang | Maximum flows through a domain[END_REF] or [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF]). In the literature, the continuous-domain duality between TV and max-flow is assumed, but in computational practice they are really different. Nevertheless, even as the method of Chambolle-Pock is solving a different optimization problem, we performed numerical comparative tests. The problem of Chambolle and Pock (now denoted CP-TV) applied to binary segmentation op-

timizes min x∈[0,1] n max F ∈R m ,||F ||∞≤1 F ⊤ (Ax) + h ⊤ x + hard constraints, (2.19)
where A ∈ R m×n is the incidence matrix of the graph of n nodes and m edges defined on the image, h ∈ R n is the metric on the nodes, defined in equation (2.5).

The CCMF algorithm is faster than both AT-CMF and CP-TV for 2D image segmentation. We have implemented CCMF in Matlab. We used an implementation of Appleton-Talbot's algorithm in C++ provided by the authors, and a Matlab software implementing CP-TV on CPU, also provided by the respective authors. The three implementations make use of multi-threaded parallelization, and the CPU times reported here were computed on a Intel Core 2 Duo (CPU 3.00GHz) processor, with 2 Gb of RAM. The CCMF average computation time on a 321 × 481 image is 181 seconds after 21 iterations. For AT-CMF, 80000 iterations require 547 seconds. For CP-TV, the average number of iterations needed for convergence is 36000, requiring an average time of 1961 seconds on CPU, and the median computation time is 1406 seconds. Note that a GPU implementation of CP-TV, also provided by the authors, achieves a speed gain factor of about 100 on a NVidia Quadro 3700. Although it was not pursued here, it should be noted that the CCMF optimization approach could also fully benefit from parallelization (e.g., on a GPU) because the core computation is to solve a linear system of equations. Realizing the benefit of a GPU solution would require the use of an iterative algorithm such as conjugate gradients or multigrid [START_REF] Bolz | Sparse matrix solvers on the GPU: Conjugate gradients and multigrid[END_REF][START_REF] Kruger | Linear algebra operators for GPU implementation of numerical algorithms[END_REF], which would make comparisons to the direct solver used here less immediate. However, previous examples in the literature have suggested that this change in solver has not created difficulties in order to realize the benefits of a GPU architecture [START_REF] Grady | Random walks for interactive organ segmentation in two and three dimensions: Implementation and validation[END_REF].

Sometimes, it may not always be necessary to wait so long for the complete convergence of CCMF (or AT-CMF) to obtain acceptable results. In cases where images exhibiting sufficiently strong gradients, one iteration may be enough to obtain a satisfying segmentation.

On such an image, one iteration of CCMF, and 100 iterations of AT-CMF show acceptable approximate results reached only after about 2 seconds for either algorithm. However, one cannot always rely on strong edges, and the power of our CCMF formulation is to behave in a predictable fashion even when edges are weak.

We may also compare the computation time of CCMF to the optimization of total variation using the Split Bregman method [START_REF] Goldstein | The split Bregman method for ℓ 1 -regularized problems[END_REF]. Optimizing TV on a 100 × 100 image requires 5000 iterations and takes 23 seconds with a sequential implementation of Split Bregman in C. On the same image, CCMF required only 17 iterations to reach the convergence criteria ||r d || < 1 and || η|| < 2, taking 4.7 seconds with a Matlab implementation (using a 2-threaded solver). We conclude that TV-based methods appear to be quite slow in the context of image segmentation, although we employed the Split Bregman algorithm which is known as one of the fastest algorithms for TV optimization for image denoising, and the very recent CP-TV algorithm. This difference in speed between the denoising and segmentation applications is due to the very large number of iterations required to convergence to a binary segmentation.

Segmentation quality

In this experiment, we compared our CCMF formulation to the AT-CMF formulation, CP-TV approach for segmentation and conventional graph cuts on the problem of image segmentation, to determine if there were strong performance differences between the formulations.

We expect that there would not be, since in principle all three formulations are trying to "minimize the cut", but have slightly different definitions of the cut length. The primary advantage that we expect from our formulation is in the reduction of metrication artifacts (as compared to conventional graph cuts) and speed/convergence (as compared to AT-CMF).

Our experiment consists of testing the Graph Cut, AT-CMF, CP-TV, and Combinatorial

Continuous Max-Flow algorithms on a database with the same seeds. We used the Microsoft 'GrabCut' database available online [START_REF] Rother | GrabCut" -Interactive foreground extraction using iterated graph cuts[END_REF], which is composed of fifty images provided with seeds. From the seed images provided, it is possible to extract two different possible markers for the background seed. Examples of such seeds are shown in Fig. 2.8(a). Different convergence criteria are available for CCMF, such as the duality gap and norms of the residuals. However, for the CMF algorithm, we do not have any satisfying criteria. Bounding the number of non binary occurrences of P does not mean that the convergence is reached, because of possible oscillations. An intermediate result after ten thousand iterations may be significantly different from the result reached after one hundred thousand. Consequently, we have run the AT-CMF algorithm until we were convinced to have reached convergence, i.e.

(a) (b) (c) (d) (e) (f) (g) (h)
when P was nearly binary and did not change significantly when we doubled the number of iterations. For half of the images in the GrabCut database, AT-CMF algorithm required more than 20000 iterations to reach convergence, for a third of the images, more than 80000 iterations, and for 1/4 of the images, more than 160000 iterations. Binary convergence was still not reached after even 500000 iterations for the rest of the images (1/4), but we stopped the computation anyway. The Chambolle-Pock TV based algorithm for binary segmentation is provably convergent and benefits also from several possible convergence criteria. The lines above the double bar show results for the first set of seeds, and the lines below the double bar show results obtained with the second set of seeds.

vergence criteria used is a ratio of changed pixels from two successive intervals of iterations being smaller than an epsilon that we fixed to 10 -7 in order to obtain satisfying results. In practice, 1/4 of the images converged in less than 60000 iterations and for the rest of them we increased the maximum number of iterations up to 200000. In contrast, CCMF only needed 21 iterations on average, and never more than 27, to reach the convergence criteria ||r d || < 1 and || η|| < 2. We noted that for all methods, segmentation quality degraded quickly if these conditions were not respected.

Table 2.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 AT-CMF results are really close, and better in mean than the GC and the CP-TV results.

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.

Unary terms

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

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

g i = exp(-β(F C i -I i ) 2
). Examples of results using these appearance priors are shown on Fig. 2.9. Each node is connected to a source node and sink node through an intermediary node which is weighted to reflect the strength of the positive and negative unary term. This construction is then applied with a simple appearance prior to automatically segment two images.

Classification

As the general CCMF formulation is defined on arbitrary graphs, it is possible to employ it in tasks beyond image segmentation. For instance CCMF may be employed to solve general problems of graph clustering, even on networks with no meaningful embedding such as social networks. We show a possible application in a classification example, considering the social network studied by Zachary [START_REF] Zachary | An information flow model for conflict and fission in small groups[END_REF]. After the split of a university karate club due to a conflict between its two leaders, Zachary's goal was to study if it was possible to predict the side joined by the members, based only on the social structure of the club. Classically, the graph is built by associating each member to a node, and edges link two members when special affinities are known between them. As weights are associated here to the strength of coupling between members, different strategies are possible for assigning the weights onto the nodes, which is necessary for using CCMF. For example, a given node/member may be assigned by the mean of affinity with its neighboring members in the graph. This weighting strategy works well in this example (See Fig. The two leaders in conflict are represented by the left hand side and right hand side nodes. Just one node is misclassified with CCMF. However, this node is known classically to be an unusual situation within the social network (see [START_REF] Zachary | An information flow model for conflict and fission in small groups[END_REF]) which is not captured by any known algorithm. 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.

3D segmentation

Conclusion

In this chapter, we have presented a new combinatorial formulation of the continuous maximum flow problem and a solution using an interior point convex primal-dual algorithm. This formulation allows us to optimize the maximum flow problem as well as its dual for which we provide an interpretation as a minimal surface problem. This new combinatorial expression of continuous max-flow avoids blockiness artifacts compared to graph cuts. Furthermore, the formulation of CMF on a graph reveals that it is actually the fact that the capacity constraints are applied to point-wise norms of flow through nodes that allows us to avoid metrication errors, as opposed to the conventional graph cut capacity constraint through individual edges.

Additionally, it was shown that our CCMF formulation provides a better approximation to the analytic catenoid than the conventional AT-CMF discretization of the continuous maxflow problem. Contrary to graph cuts, we showed that CCMF estimates the true boundary length of objects. Finally, unlike finite-element based methods such as AT-CMF, the CCMF formulation is expressed for arbitrary graphs, and thus can be employed in a large variety of tasks such as classification.

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

We have also compared our algorithm to known weighted combinatorial TV minimization methods.

We have shown that our results are generally better for segmentation than combinatorial TV-based methods, and that our implementation on CPU is much faster and more predictable. Specifically, the number of iterations required for convergence does not vary much for CCMF, whereas it can vary of more than a factor of ten for CP-TV. The deep study of the relationships between CCMF and combinatorial TV reveals that, in contrast with expectations from their duality under restrictive assumptions in the continuous domain, their duality relationship is only weak in the combinatorial setting. In fact max-flow and total variation are different problems with different constraints, yielding different algorithms and different results. One key difference between max-flow and total variation is that max-flow algorithms were developed as segmentation algorithms, following the graph cuts framework.

One consequence is that max-flow formulations have a null divergence objective, which is not present in total variation formulations. Null divergence can be viewed as a consequence or a necessity in order to obtain constant partitions almost everywhere, in particular binary ones. This difference is important because in the proposed CCMF framework we impose a tight null divergence constraint, which to our knowledge is novel for an isotropic formulation (graph-cuts have always had this constraint). AT-CMF and similar frameworks achieve null divergence if and when convergence is achieved. There is no null divergence constraint at all in total variation frameworks. As a consequence, while in practice it is possible to compare total variation and max flow formulations, they should be treated differently. However, the strong computational performance and segmentation quality results of CCMF as compared to TV (using the strongest known optimizations methods) suggests that it may be advantageous to apply our CCMF formulation to problems for which TV has proven effective (such as filtering).

Several further optimizations of CCMFs are possible, for instance using multi-grid 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.

In the next chapter, inspired by the CCMF formulation, we propose a combinatorial formulation for the weighted Total Variation primal dual formulation. In contrast to CCMF, the new TV-based model is adapted to multi-label problems.

Chapter 3

Dual Constrained TV-based regularization framework

Total Variation minimization based algorithms are prevalent in computer vision. They are used in a variety of applications such as image denoising, compressive sensing and inverse problems in general. In this work1 , we extend the TV dual framework that includes Chambolle's and Gilboa-Osher's projection algorithms for TV minimization. We use a flexible graph data representation that allows us to generalize the constraint on the projection variable. We show how this new formulation of the TV problem may be solved by means of fast parallel proximal algorithms. We show on denoising and deblurring examples that the proposed approach performs better than classical TV. We also show that the proposed method is applicable to a variety of other inverse problems including image fusion and mesh filtering.

Introduction

The Total Variation (TV) model was introduced as a regularization criterion for image denoising by Rudin, Osher and Fatemi (ROF) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], and has been shown to be quite efficient for smoothing images while preserving contours. Total variation minimization is used for solving a large number of inverse problems, including deconvolution [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF], inpainting, MRI reconstruction, optical flow estimation [START_REF] Bagnato | Optical flow and depth from motion for omnidirectional images using a TV-ℓ 1 variational framework on graphs[END_REF], stereovision [START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF], among others. A major advantage of the ROF-TV minimization problem is that it is a convex problem. However, the non-differentiability of the involved objective function requires the use of specific optimization techniques, particularly with respect to speed and efficiency. Much progress has been achieved by employing primal-dual approaches [START_REF] Carter | Dual Methods for Total Variation-based Image Restoration[END_REF][START_REF] Chan | A nonlinear primal-dual method for total variation based image restoration[END_REF], the dual approach of Chambolle [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF],

and fast first order methods [START_REF] Aujol | Some first-order algorithms for total variation based image restoration[END_REF].

Qualitative improvements have also been obtained by introducing a weighted model [START_REF] Bougleux | Discrete regularization on weighted graphs for image and mesh filtering[END_REF][START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF][START_REF] Gilboa | Nonlocal convex functionals for image regularization[END_REF][START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF]. In this model, a discrete TV energy is optimized in edge-weighted graphs.

The methods described in [START_REF] Gilboa | Nonlocal convex functionals for image regularization[END_REF] and [START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF] employ a non-local weighting strategy following the non-local means concept [START_REF] Buades | A non-local algorithm for image denoising[END_REF]. Furthermore, Gilboa and Osher provided an efficient dual algorithm in [START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF], based on Chambolle's [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] method to address this problem. Currently, one of the fastest methods for optimizing weighted TV is an alternating direction of multipliers method [START_REF] Goldstein | The split Bregman method for ℓ 1 -regularized problems[END_REF][START_REF] Zhang | Bregmanized nonlocal regularization for deconvolution and sparse reconstruction[END_REF], related to augmented Lagrangian approaches [3], and also called split-Bregman. Another (node) weighting strategy, and a fast primal dual algorithm, is proposed in [54].

Recently, several approaches were proposed to reduce the staircaising effect of total variations methods. The Total Generalized Variation method of Bredies et al. [START_REF] Bredies | Total generalized variation[END_REF][START_REF] Bredies | Properties of l 1 -TGV 2 : The onedimensional case[END_REF] employs higher order derivatives of the function to regularize. The benefit is to replace the staircase effect by smoothed variations, while preserving contours. The work of Jalalzai and Chambolle aims at achieving the same goal [START_REF] Jalalzai | Enhancement of blurred and noisy images based on an original variant of the total variation[END_REF]. In the hybrid TV approach of Combettes and Pesquet [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF](See [START_REF] Pustelnik | Parallel proximal algorithm for image restoration using hybrid regularization[END_REF] for further investigations) an energy combining wavelets and TV is optimized to jointly reduce the staircase effect of both TV and wavelet artifacts.

In addition to the staircase effect, TV approaches typically result in a loss of contrast. A few method have been proposed [START_REF] Ciril | Image denoising with a constrained discrete total variation scale space[END_REF][START_REF] Darbon | Image restoration with discrete constrained total variation part I: Fast and exact optimization[END_REF] for contrast preservation.

The projection algorithms in [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF] are based on a relatively simple local constraint on the norm of the projection variable. In this work, we extend the constraint on this variable.

This extension allows us to better adapt the optimization procedure to the local information.

We name this new approach "Dual-Constrained Total Variation" (DCTV) regularization.

Our main contributions are the following:

1. By using an appropriate reformulation and decomposition of our energy functional, we are able to optimize it efficiently, owing to the use of fast parallel proximal algorithms.

2. By formulating the primal dual total variation problem on arbitrary graphs, we are able to use our regularization method on graphs that are not uniformly sampled, for example on a graph representations of the fovea or an arbitrary mesh.

3. By providing a formulation taking into account weights both on nodes and edges of the graph, we obtain a more flexible regularization.

4. These improvements lead to better results, in particular a better contrast preservation and a more detailed recovery in image processing examples.

Proximal methods provide efficient solutions to convex optimization problems where the objective function can be split into a sum of convex functions [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF]. This core of methods are applicable to large size problems, as well as to non-smooth and non-finite functions.

Note that the constrained TV-based energy optimized in this work is different than the one considered in previous works [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF][START_REF] Malgouyres | Minimizing the total variation under a general convex constraint for image restoration[END_REF]. Combettes and Pesquet [START_REF] Combettes | Image restoration subject to a total variation constraint[END_REF] introduce a total variation bound as a constraint in a convex energy minimization problem. Malgouyres [START_REF] Malgouyres | Minimizing the total variation under a general convex constraint for image restoration[END_REF] defines a constraint in order to combine the ROF-TV model with wavelet-based soft thresholding methods.

The chapter is organized as follows: in Section 3.2, we generalize TV models on a graph by extending the dual formulation under a constrained form. In Section 3.3, we demonstrate how our constrained TV-based optimization problem can be efficiently solved by using proximal algorithms. Finally, results obtained with the proposed approach are presented in four applications and compared with classical weighted TV.

Graph extension of TV models

Given an original corrupted image f , the purpose of variational methods for image restoration is to deduce a restored image u close to the original image f under the assumption of smooth variations of intensity values inside objects. Let λ ∈]0, +∞[ be a real positive value representing a regularization parameter. In a functional setting, given a planar domain Ω, and denoting by u and v two arbitrary points of Ω, the weighted TV model [START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF] is given by min

x Ω Ω w u,u ′ (x u ′ -x u ) 2 du ′ 1/2 du + 1 2λ Ω (x u -f u ) 2 du, (3.1)
where w is a non-negative valued function defined on Ω 2 . As shown by Chan et al. in [START_REF] Chan | A nonlinear primal-dual method for total variation based image restoration[END_REF] the TV minimization problem (3.1) is equivalent to the min-max problem min

x max ||p||∞≤1 Ω 2 w 1/2 u,u ′ (x u ′ -x u )p u,u ′ du du ′ + 1 2λ Ω (x u -f u ) 2 du, (3.2)
where p a two-variable function and

p ∞ = sup u∈Ω Ω p 2 u,u ′ du ′ 1/2 .
Let us now establish the formulation of the DCTV energy within a discrete framework, using the notations specifyied in the introduction chapter. Now, f = (f i ) 1≤i≤n and x = (x i ) 1≤i≤n are vectors in R n . The combinatorial TV model as defined by [START_REF] Bougleux | Discrete regularization on weighted graphs for image and mesh filtering[END_REF][START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF] is

min x∈R n i∈J j∈N i w i,j (x j -x i ) 2 1/2 + 1 2λ n i=1 (x i -f i ) 2 (3.3)
where, for every i ∈ {1, . . . , n}, N i = {j ∈ {1, . . . , n} | e i,j ∈ E}, J = {i ∈ {1, . . . , n} | N i = ∅}. The dual formulation, which can be optimized by a projection algorithm [START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF], is given by min

x∈R n max p ∞≤1 p ⊤ (Ax) • √ w + 1 2λ x -f 2 , (3.4)
where • denotes the Euclidean norm, w ∈ R m is a vector with components w i,j for every couple of nodes (i, j) with i ∈ J and j ∈ N i and similarly denoting by p i,j the components of

p ∈ R m , p ∞ = max i∈J j∈N i p 2 i,j 1/2 . (3.5)
By introducing a vector F ∈ R m with components F i,j such that, for every i ∈ J and j ∈ N i ,

F i,j = p i,j √ w i,j
, the problem can be reformulated as min

x∈R n max F ∈C F ⊤ Ax + 1 2λ x -f 2 , (3.6)
C = {(F ∈ R m | (∀i ∈ J) j∈N i F 2 i,j w i,j ≤ 1}. (3.7)
The objective of this work is to extend the discrete weighted TV variational formulation in (3.6) by investigating the following optimization problem min

x∈R n sup F ∈C F ⊤ Ax + 1 2 (Hx -f ) ⊤ Λ -1 (Hx -f ) (3.8)
where f ∈ R q is an observed vector of data, H ∈ R q×n and Λ is a weighting symmetric definitepositive matrix in R q×q . The metrix H may simply be the identity matrix for image denoising, a convolution operator in restoration tasks, or a projection matrix in reconstruction problems.

The matrix Λ may be proportional to the covariance matrix of the noise corrupting the data, as commonly used in weighted least squares approaches. In addition, the main contribution of this work is to consider more general convex sets C than those given by (3.7).

Considered class of constraint sets

The proposed optimization approach will allow us to address nonempty convex sets C that can be decomposed as an intersection of closed convex subsets (C r ) 1≤r≤s of R m , the projections onto which take closed forms. An example of a set C of interest is given by

C = s r=1 C r (∀r ∈ {1, . . . , s}) C r = {F ∈ R m | (∀i ∈ S r ) θ (i) • F α ≤ g i }. (3.9)
where (S r ) 1≤r≤s is a partition of {1, . . . , n},

• α is the ℓ α norm of R m with α ∈ [1, +∞]
and, for every i ∈ {1, . . . , n}, θ (i) ∈]0, +∞[ m is a vector of multiplicative constants for every couple of node (i ′ , j) with i ′ ∈ J and j ∈ N i ′ .

The form of C in (3.9) obviously includes (3.7) as a particular case where α = 2 and

(∀i ∈ {1, . . . , n}) (∀i ′ ∈ J)(∀j ∈ N i ′ ) θ (i) i ′ ,j =    1 √ w i ′ ,j if i ′ = i and j ∈ N i ′ 0 otherwise g i = 1. (3.10)
However, the proposed approach offers much more flexibility.

In this chapter, we are mainly interested in the case when for every i ∈ {1, . . . , n}, θ (i) is the i-th line vector of |A ⊤ |. Specifically, C may be defined as

C = {F ∈ R m | g •2 -|A ⊤ |F •2 ∈ [0, +∞[ n }, (3.11)
where g = (g i ) 1≤i≤n and, for every vector a, a to multi-label problems.

Concerning the choice of the node weight g i at the i-th node, a simple strategy consists of considering a monotonically decreasing function of the data gradient. More specifically, given positive reals ǫ and χ, we suggest using g defined as

(∀i ∈ {1, . . . , n}) g i = exp(-χ ∇x i 2 ) + ǫ, (3.12)
where x is some reference data defined on a graph, for instance image data, which corresponds to some rough estimate of x, and

∇x i 2 =    ( j∈N i (x i -x j ) 2 ) 1 2 if i ∈ J 0 otherwise (3.13)
corresponds to the Euclidean norm of its discrete gradient ∇x i at node i. In the absence of a contour, g i takes large values, so are the components of F corresponding to nonzero values of θ (i) , preventing large local variations of x in the minimization (3.8). Conversely, in the presence of a contour, g i ≃ ǫ, and the components of F corresponding to nonzero values of θ (i) are small, thus allowing large local variations of x. In image filtering applications, in addition to intensities information, g may be used to penalize changes in other relevant image quantities such as color or texture.

Proposed algorithms

We show in this section that it is possible to efficiently solve Problem (3.8) by proximal methods [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. To do so, define the support function σ C of the closed convex constraint set C as

σ C : R m →] -∞, +∞] : a → sup F ∈C F ⊤ a. (3.14)
This is a proper lower-semicontinuous convex function, the conjuguate of which is the indi-

cator function of C, ı C : F →    0 if F ∈ C, +∞ otherwise. (3.15)
This leads us to consider the following optimization problem:

min x σ C (Ax) + 1 2 (Hx -f ) ⊤ Λ -1 (Hx -f ) + η 2 Kx 2 (3.16)
where η ∈]0, +∞[ and K ∈ R n×n is the projection matrix onto the nullspace of H, specifically,

K = I -H ⊤ (HH ⊤ ) -1 H. When H is injective (rank H = n),
min F ϕ(F ) + ı C (F ), (3.17)
where ϕ :

F → 1 2 F ⊤ AΓA ⊤ F -F ⊤ AΓH ⊤ Λ -1 f and Γ = (H ⊤ Λ -1 H + ηK) -1 .
The optimal solution to the primal problem (3.16) is deduced from any optimal solution F of the dual problem by the relation

x = Γ H ⊤ Λ -1 f -A ⊤ F . (3.18)
Note that the dual forward-backward algorithm proposed in [START_REF] Combettes | Dualization of signal recovery problems[END_REF] is not applicable to Problem 3.17 since the projection onto C is not explicit in general. Dykstra's algorithm [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in hilbert spaces[END_REF] to compute exact projections on an intersection of convex sets would not be applicable either, as the operator AΛA ⊤ is singular. In order to numerically solve (3.17), recall that C = s r=1 C r , so that ı C can be decomposed into the sum of the indicator functions of the convex subsets (C r ) 1≤r≤s . Hence, the problem is equivalent to solving min

F ∈R m s r=1 ı Cr (F ) + ϕ(F ).
(3.19)

Parallel proximal algorithm (PPXA)

The above sum of (s + 1) convex functions can be efficiently optimized by resorting to the Parallel ProXimal Algorithm (PPXA) proposed in [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]. PPXA was generalized in [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF],

where links with augmented Lagrangian methods [3] were established.

As shown in

π r,k = P Cr (y r,k ) if r ≤ s (γAΓA ⊤ + I) -1 (γAΓH ⊤ Λ -1 f + y s+1,k ) otherwise z k = 2 s+1 (π 1,k + • • • + π s+1,k ) -F k for r = 1, . . . , s + 1 do in parallel y r,k+1 = y r,k + ν(z k -p r,k ) F k+1 = F k + ν 2 (z k -F k ) k = k + 1 until convergence
Let us now turn our attention to convex sets as described in (3.9). The main difficulty for applying PPXA is to find a partition (S r ) 1≤r≤s for which the projections (P Cr ) 1≤r≤s can be computed easily. To ensure this property, the partition can be defined in such a way that

(∀r ∈ {1, . . . , s})(∀(i, i ′ ) ∈ S 2 r ) i = i ′ ⇒ #(θ (i) • θ (i ′ ) ) = 0, (3.20)
where #(a) denotes the number of nonzero components in a vector a. Then, it is easy to see that the projection onto C r with r ∈ {1, . . . , s} of a vector F can be decomposed as a finite number of projections of subvectors of dimensions #(θ (i) ) i∈Sr extracted from F onto weighted ℓ α balls.

For a 4-connected lattice, when (∀i ∈ {1, . . . , n}) θ (i) is the i-th line of matrix |A ⊤ |, a decomposition of C can be performed by setting s = 2. The partition (S r ) 1≤r≤2 corresponds to two spatial disjoint sets each with a checkered pattern (also called red black checkerboard pattern). The projection onto C r with r ∈ {1, 2} for α = 2 (resp. α = 1 or α = +∞) reduces to simple projections onto hyperspheres (resp. hypercubes [START_REF] Van Den | Probing the Pareto frontier for basis pursuit solutions[END_REF]).

For example, the decomposition on a small lattice when α = 2 is illustrated by Figure 3.1.

g 1 g 2 g 3 g 4 g 5 g 6 F1,2 F2,3 F4,5 F5,6 F1,4 F2,5 F3,6 g 1 g 2 g 3 g 4 g 5 g 6 C 1 ={F | F 2 1,2 + F 2 2,3 + F 2 2,5 ≤ g 2 2 , F 2 1,4 + F 2 4,5 ≤ g 2 4 , F 2 3,6 + F 2 5,6 ≤ g 2 6 } C2 ={F | F 2 1,2 + F 2 1,4 ≤ g 2 1 , F 2 2,3 + F 2 3,6 ≤ g 2 3 , F 2 2,5 + F 2 4,5 + F 2 5,6 ≤ g 2 5 } Figure 3.1: Decomposition of C = C 1 ∩ C 2
Note also that, at each iteration k ∈ N, the computation of π s+1,k requires a matrix inversion. In cases where H and Λ are (or can be well approximated by) circulant-block circulant matrices, and the graph is regular, this step can be performed efficiently [START_REF] Mazancourt | The inverse of a block-circulant matrix[END_REF][START_REF] Fischer | Fast inversion of matrices arising in image processing[END_REF]. For this purpose, an equivalent formulation of matrix (γAΓA ⊤ + I) -1 is computed using the matrix inversion lemma, also known as Woodbury's matrix identity. The new resulting linear system is then composed of a circulant block circulant matrix and is solvable efficiently using a Fast Fourier Transform.

Decomposing arbitrary graphs

In arbitrary graphs, the graph decomposition of Figure 3.1 in a checkerboard pattern is no longer valid, so we use a graph coloring algorithm for the decomposition. Coloring a graph with a minimum number of colors for the nodes is an NP-complete problem. However, there exist heuristic algorithms that find a non-optimal, but feasible solution in polynomial time. In our experiments, we employed the DSAT (Degree SATuration) algorithm [START_REF] Brélaz | New methods to color the vertices of a graph[END_REF],

which is guaranteed to color triangulated graphs in linear time. An illustration is provided in Figure 3.2. Experimentally, this greedy algorithm typically decomposes of the convex set C in five sets or fewer.

M+SFBF algorithm

Starting from Problem (3.16), we deduce an alternative formulation of the energy to minimize.

The formulation we obtain shall allow us to employ another parallel proximal algorithm, the Monotone + Skew Forward Backward Forward (M+SFBF) algorithm, which differs from PPXA. Interestingly, M+SFBF [START_REF] Arias | A monotone+skew splitting model for composite monotone inclusions in duality[END_REF] does not require us to solve a linear system at each iter-

min x∈R n ,a 1 ∈R m ,...,as∈R m a 1 +•••+as=Ax s r=1 σ Cr (a r ) + 1 2 (Hx -f ) ⊤ Λ -1 (Hx -f ) + η 2 Kx 2 . (3.21)
Proof. Noticing that

σ * C = ι C = s r=1 ι Cr , (3.22)
we deduce that

σ C = ι * C 1 • • • ι * Cs = σ C 1 • • • σ Cs (3.23)
where denotes the infimal convolution. The infimal convolution of two convex functions ϕ and ψ from R m to ] -∞, +∞] is the convex function defined as

(∀a ∈ R m ) (ϕ ψ)(a) = inf a ′ ∈R m ϕ(a ′ ) + ψ(a -a ′ ). (3.24)
We have then

(∀a ∈ R m ) σ C (a) = inf a 1 ∈R m ,...as∈R m a 1 +•••+as=a σ C 1 (a 1 ) + • • • + σ Cs (a s ). (3.25)
This allows us to deduce that (3.16) can be reformulated as the following multivariate optimization problem:

min x∈R n ,a 1 ∈R m ,...,as∈R m a 1 +•••+as=Ax s r=1 σ Cr (a r ) + 1 2 (Hx -f ) ⊤ Λ -1 (Hx -f ) + η 2 Kx 2 .
Algorithm 6: M+S Forward-backward-forward algorithm solving (3.21)

Fix ε ∈]0, (s + 1)β + 1 -1 [ where β = max{1, A }. Set k = 0. Choose (F r,0 ) 1≤r≤s+1 and (a r,0 ) 1≤r≤s in R m and x 0 ∈ R n . repeat F k = 1 s+1 s+1 r=1 F r,k γ k ∈ [ε(s + 1), (1 -ε)/β] for r = 1, . . . , s + 1 do in parallel G r,k = F r,k -γ k a r,k if r ≤ s F s+1,k + γ k Ax k otherwise. b r,k = a r,k + γ k F r,k if r ≤ s x k -γ k A ⊤ F s+1,k otherwise. Π k = 1 s+1 s+1 r=1 G r,k for r = 1, . . . , s + 1 do in parallel π r,k = b r,k -γ k P Cr (γ -1 k b r,k ) if r ≤ s γ k (H ⊤ Λ -1 H + ηK) + I -1 (γ k H ⊤ Λ -1 f + b s+1,k ) otherwise G ′ r,k = Π k -γ k π r,k if r ≤ s Π k + γ k Aπ s+1,k otherwise. b ′ r,k = π r,k + γ k Π k if r ≤ s π s+1,k -γ k A ⊤ Π k otherwise. F r,k+1 = F r,k -G r,k + G ′ r,k a r,k+1 = a r,k -b r,k + b ′ r,k if r ≤ s x k+1 = x k -b s+1,k + b ′ s+1,k otherwise k = k + 1 until convergence
We see that the proposed approach amounts to decomposing the gradient Ax of a candidate solution x in a sum of components (a r ) 1≤r≤s onto which a simpler regularization is performed. In particular, as in the previous section, we choose the decomposition of C in such a way that each C r is the product of convex sets in a lower-dimensional space, i.e. × i∈Sr C i r (for example, C i r is a hypersphere of R 4 ). In this case, σ C i r with i ∈ S r is a "local" sparsity measure, in the sense that the proximity operator of each σ C i r is a proximal thresholder [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]: One additional advantage of this primal-dual algorithm is that the convergence of the sequence (a 1,k , . . . , a r,k , x k ) k∈N toward a solution to Problem 3.21 is guaranteed even in the limit case when η = 0.

prox σ C i r a i r = 0 ⇔ a i r ∈ C i r . ( 3

Results

We now demonstrate the behavior of DCTV with respect to different choices of convex sets, weighting strategies, and graph constructions.Due to the generality of the proposed framework, we can employ DCTV in various applications such as image denoising, image deconvolution, image fusion, and mesh filtering. In our first application, we demonstrate the performance of DCTV with respect to weighted TV [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF] for image denoising, and show that the contrast of images is better preserved using DCTV. In this section 3. The weights are set according to (3.12) with χ = 0.04 for DCTV. The stopping criteria for the algorithms is x k -x k-1 < 5.10 -3 x k . The results obtained with the weighted TV [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF] look slightly better than using the classical, unweighted TV model [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. Performing the DCTV optimization using convex sets with an ℓ ∞ norm results in more artifacts, however an improvement is obtained using ℓ 2 -norm based convex sets of (3.11) with PPXA.

In a second experiment, we compare quantitatively DCTV with weighted TV. We used four standard test images that we corrupted with synthetic white zero-mean Gaussian noise with variance σ 2 . The value of the image fidelity parameter λ was set according to an empirical rule depending on the variance of the noise, and ǫ was set to λ. Signal to Noise Ratio (SNR) is used as the performance measure in our quantitative evaluation. Xeon 2.5GHz 8-core system.

Non local denoising

Rather than using locally connected graphs, non-local strategies [START_REF] Bresson | A short note for nonlocal TV minimization[END_REF][START_REF] Buades | A non-local algorithm for image denoising[END_REF] have been shown

to achieve denoising improvements. A non-local strategy may naturally be employed in the DCTV framework. The weights between non-neighbor nodes are computed following the main idea of [START_REF] Buades | A non-local algorithm for image denoising[END_REF]: for each pixel p of the image, the squared sum of differences (SSD) between the intensities in a block around p and all other blocks in a large neighborhood around that block is computed. Then, edges are added between p and the nodes producing the best The noise was AWGN with variance σ 2 = 35. The graphic shows the profile gray levels on one vertical line of the original image (in black), the TV result (in green), and DCTV result (in red).

matchs, with weights corresponding to the normalized SSD scores. As a result, the Laplacian -related to the matrix to be inverted at each iteration in PPXA -of this non-local graph is no longer circulant block circulant, resulting in loss of efficiency. We propose to employ the M+SFBF algorithm (Algorithm 6) in this case. The nonlocal weights are computed using [START_REF] Bresson | A short note for nonlocal TV minimization[END_REF]. Regularization parameters used to obtain the best results: (d) λ = 0.017 (e) λ = 0.015 (f) λ = 0.85.

Comparison of applicability of PPXA and M+SFBF

The proximity operator associated with the data fidelity term in PPXA requires solving a linear system. In the case of an image recovery task using a local, regular graph (for example a 4-connected grid) as in the first part of Section 3.4.1, the special configuration of the matrix in the system may be taken into account. As this resulting matrix is circulant block-circulant, the system is efficiently solved using a Fast Fourier Transform. However, in the case of an arbitrary graph, for example a non-local graph as in Section 3.4.1.2, the matrix of the linear system is no longer circulant, so the system may be solved using a conjugate gradient method for instance.

When using the M+SFBF algorithm, the proximity operator of the function enforcing the data fidelity does not require us to solve a linear system. Instead, four matrix multiplications by the incidence matrix are performed per iteration. The speed of convergence using M+SFBF and PPXA are compared in Figure 3.7, for denoising an image using a local, and a non-local graph. In the figure, only computation times are presented because the required number of iterations to reach convergence is the same for either algorithms. In a local graph, the matrix multiplications in M+SFBF make the algorithm slower than PPXA. Conversely, in a non-local graph, where the linear system of PPXA is solved using a conjugate gradient, M+SFBF converges much faster.

Image deconvolution

We now give some quantitative comparison results for joint denoising and deblurring tasks.

We compare the DCTV results to Wiener based deconvolution using the Matlab function ''deconvwrn''. The comparison also includes the hybrid TV/wavelet regularization method of Combettes and Pesquet [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]. We report in Table 3.4.2 the SNR values for restoration of images corrupted with additive white zero-mean Gaussian noise with variance σ 2 = 5, and σ 2 = 10, and convolved with uniform blur kernels of size 5 × 5 and 7 × 7. We observe that DCTV unsurprisingly outperforms the standard Wiener filter. More importantly, DCTV is competitive with a state-of-the-art method [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF], both quantitatively in term of SNR, and qualitatively, without presence of checkerboard artifacts observed for the Hybrid TV method in Figure 3 (3.12), where the estimate of the image x has been precomputed using DCTV with the original corrupted image f . This trick allows us to obtain an SNR improvements of about 0.5 dB, while computation times are still twice as fast than the hybrid TV method [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]. We observe in the TV hybrid result the presence of checkerboard artifacts due to the use of discrete filters for approximating the gradient. DCTV reduces the staircase effect of TV while preserving more details. Parameter used: Hybrid TV with regularization parameters [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF]α = 0 and β = 0.025, and DCTV with λ = 0.005 and η = 0.04.

Image fusion

We show here how to apply our Dual Constrained TV approach to image fusion. From N degraded images f 1 , . . . , f N , the problem is to find a restored image x. The degradation model is the following:

    f 1 . . . f N     =     H 1 . . . H N     x +     b 1 . . . b N     , (3.27)
where H 1 , ..., H N are degradation matrices (i.e. generating blur), and b 1 , ..., b N represent white additive zero-mean noises possibly with different variances, assumed to be uncorrelated.

The problem can be formulated as in (3.8) by setting The M+SFBF algorithm was used to optimize DCTV with λ = 0.25.

f =     f 1 . . . f N     , H =     H 1 . . . H N     , Λ =        λ 1 I 0 . . . 0
On a rectangular grid, where the convex C is decomposed as described in Figure 3.1, the parallel proximal algorithm given in Algorithm 5 can be employed efficiently in this case.

Graph filtering

3D mesh filtering

The dual-constrained TV-based formulation being defined on arbitrary graphs, x is not limited to only represent image pixels. In Figure 3.10, we present an example of mesh denoising, where x is a vector composed of the spatial coordinates x X , x Y and x Z of the mesh nodes. In this experiment, we added a randomly oriented with white zero-mean Gaussian magnitude noise vector to the original nodes of a mesh. This results in noisy mesh nodes of coordinates f X , f Y and f Z . The degradation model is the following: For this application, we used the M+SFBF algorithm, which is more efficient than PPXA as the graph is not regular. This application shows that the DCTV framework is well suited for regularizing various type of data.

f =    f X f Y f Z    =    x X x Y x Z    + σ 2    b X b Y b Z    , (3.29

Biologically sampled image filtering

This final example demonstrates the ability of DCTV to efficiently filter data on arbitrary graphs. In some situations, images are acquired from nonuniform samples. Specifically, most biological systems are known to acquire light nonuniformly. Following the work of [START_REF] Grady | Space-Variant Computer Vision: A Graph-Theoretic Approach[END_REF], in Results obtained with DCTV are compared to results obtained using the mean filtering and lowpass filtering method of Taubin et al. [START_REF] Taubin | Optimal surface smoothing as filter design[END_REF]. We note that DCTV improves contrast preservation and thin objects rendering compared to Taubin's method.

Conclusion

In this chapter we have extended existing TV models by generalizing the constraint on the projection variable of the dual TV formulation. This new approach shows improved results compared with the weighted TV approach in image restoration applications. More generally, the proposed proximal algorithms make it possible to efficiently solve convex minimization problems involving the support function of an intersection of convex sets as a penalty term.

It is also worth emphasizing that this approach can be applied on any graph data structures, in particular in those frequently employed in 3D modeling. By taking into account both edge and node weights, we have shown that the proposed framework brings improvements in denoising and deconvolution tasks. Its versatility also allows us to exploit extra information such as non local patch similarity into weight computation that makes it even more appealing.

Chapter 4

A unifying graph-based optimization framework: Power watershed

Introduction

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 [START_REF] Beucher | The morphological approach to segmentation: The watershed transformation[END_REF] shares a similar seeding interface but only recently was a connection made between the watershed algorithm and graph cuts [4,5]. In this chapter 1 , 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 chapter, 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.

1 Published in Proc. of ICCV, 2009 [START_REF] Couprie | Power watersheds: a new image segmentation framework extending graph cuts, random walker and optimal spanning forest[END_REF], and IEEE Trans. on PAMI, 2011 [START_REF] Couprie | Power Watersheds: A Unifying Graph Based Optimization Framework[END_REF] 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.

Defining an energy function for the watershed optimization allows us to provide an

MRF interpretation for the watershed.

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. We refer to Chapter 1, and references therein for a more complete description.

A short review of graph-based segmentation

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] Bertrand | On topological watersheds[END_REF][START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF][START_REF] Cousty | Watershed cuts: thinnings, shortest-path forests and topological watersheds[END_REF][START_REF] Najman | Watershed of a Continuous Function[END_REF][START_REF] Roerdink | The watershed transform: Definitions, algorithms, and parallellization strategies[END_REF][START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[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][START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF]. 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 chapter, 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 [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF]. 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 quasi-linear 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 [START_REF]Mathematical Morphology: from theory to applications[END_REF] 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 [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in N-D images[END_REF], 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] Li | Lazy snapping[END_REF][START_REF] Rother | GrabCut" -Interactive foreground extraction using iterated graph cuts[END_REF]. Although GC is relatively new, the use of minimal surfaces in segmentation has been a common theme in computer vision for a long time [START_REF] Blake | Visual Reconstruction[END_REF][START_REF] Geman | Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images[END_REF][START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] and other boundary-based user interfaces have been previously

employed [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF][START_REF] Falcão | User-steered image segmentation paradigms: Live wire and live lane[END_REF][START_REF] Grady | Minimal surfaces extend shortest path segmentation methods to 3D[END_REF][START_REF] Mortensen | Interactive segmentation with intelligent scissors[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 [START_REF] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF] 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] Appleton | Globally minimal surfaces by continuous maximal flows[END_REF][START_REF] Boykov | Computing geodesics and minimal surfaces via graph cuts[END_REF][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 [START_REF] Coifman | Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps[END_REF], as a semi-supervised transduction learning algorithm [START_REF] Duchenne | Segmentation by transduction[END_REF] or as an interactive version of normalized cuts [START_REF] Grady | Fast approximate random walker segmentation using eigenvector precomputation[END_REF][START_REF] Shi | Normalized cuts and image segmentation[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 4.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.

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 [6,[START_REF] Criminisi | GeoS: Geodesic image segmentation[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.

A unifying energy minimization framework

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 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 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]. 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. We recall that 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 ), (4.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 e ij ∈E (w ij |x i -x j |) q + v i ∈V (w i |x i -y i |) q (4.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, (4.2) may be included in the following algorithm

Step 1: x = arg min

x e ij ∈E (w ij |x i -x j |) q + v i (w F i |x i |) q + v i (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 . (4.3)
In other words, we are looking for an optimum x ⋆ of Eq. ( 4.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) max-flow (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 (4.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 e ij ∈E w p ij |x i -x j | q + v i w p F i |x i | q + v i 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.4)
As before, the final segmentation s being chosen via a threshold.

We observe that (4.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 [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF] for more details):

min x λ e ij ∈E w p ij |x i -x j | q smoothness terms + v i ∈V w p i |x i -y i | q
data fidelity terms (4.5)

For example, the unary term w p B i |x i -1| q can also be rewritten as w p i |x i -y i | q , where y i is an auxiliary node and the signal at this auxiliary node is fixed at y i = 1.

As with (4.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.4) can be computed by the shortest path (geodesics) algorithm. These three algorithms form the underpinning for many of the advanced image segmentation methods in the literature.

It was shown in [4,5] that when q = 1 (graph cuts) and p → ∞ then the solution of (4.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 (4.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 [START_REF] Chazelle | A minimum spanning tree algorithm with inverse-ackermann type complexity[END_REF]. In the next sections we explore two cases that are, to the best of our knowledge, unexplored. First, we show in Appendix A.2.1 that the 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.4) for any value of q when p → ∞. Table 4.1 gives a reference for the different algorithms generated by various value of p and q.

The case q finite, p → ∞ leading to watershed

We now generalize the link between GC and MSF established by Allène et al. [4,5] 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. 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.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.

of Y , 0 to the vertices connected to a vertex of B, and 1 to the vertices connected to a vertex of F . Let F be a subgraph of G. We say that F is a spanning forest (relative to F ∪ B) if:

Examples of extensions appear in

(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 = e ij ∈F w ij .
We say that a spanning forest F is a maximum spanning forest (MSF) (relative to F ∪ B)

for w if the weight of F is maximum, i.e. greater or equal to the weight of any other spanning forest (relative to F ∪ B).

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.

These definitions are compatible with the watershed cut framework of [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF]. We may now introduce a general link between the Maximum Spanning Forest segmentation result and the solution of the optimization of (4.4) when the power of the weights converges toward infinity. 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. 4.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

e ij ∈E w ij w max p |x i -x j | q , s.t. x(F ) = 1, x(B) = 0. (4.6)
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 q-cut 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 (4.6). 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 4.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., 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.

In theorem 4.3.1, the condition for seeds to be the maxima of the weight function is necessary as specified in Appendix A.2.2.

We can note that if the weights are all different, the MSF cut is unique and Thm. 4.3.1 is also true without the condition for seeds to be the maxima of the weight function. We can also note that in the case q = 1, the authors of [5] explicitely provide a lower bound for p such that the min-cut corresponds to an MSF-cut.

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, and the proof is given in Appendix A.2.3.

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

An interpretation of the minimization of our general energy as a maximum a posteriori approximation is presented in Appendix A.2.5.

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.

4.4 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 e ij ∈ plateau |x i -x j | q (4.7)
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. 4.2.

In Algorithm 7, 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.4) are treated as binary terms connected to Theorem 4.4.1. 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 ≤ δ. (4.9)
The proof of this theorem is given in Appendix A.2.4.

Theorem 4.4.2. If q > 1, the potential x * obtained by minimizing the energy of (4.4)

subject to the boundary constraints converges toward the potential x obtained by Algorithm 7.

as p → ∞.

Proof. We prove that by optimizing (4.4) we are performing the same steps as Alg. 7. As in Thm. 4.3.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 (4.6).

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 (4.6). This choice of xi = xj is also performed by Algorithm 7 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 (4.6) boils down to minimizing Repeating the steps recursively until all the vertices are labeled, we find that the Algorithm 7 procedure agrees with the energy minimization of (4.5).

e ij ∈Emax |x i -x j | q , ( 4
We can note that even if Algorithm 7 minimizes the energy E p,q in the case p → ∞, several solutions x are possible when q = 1.

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

Proof. At each step of Algorithm 7, we consider a set of connected edges of maximum weight

E max .
If E max contains no vertices of different labels, Algorithm 7 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 7 recursively until all the vertices are labeled, we find that we are building an MSF cut for w.

Uniqueness of the 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 [START_REF] Beucher | The morphological approach to segmentation: The watershed transformation[END_REF][START_REF] Bieniek | An efficient watershed algorithm based on connected components[END_REF][START_REF] Couprie | Quasi-linear algorithms for the topological watershed[END_REF][START_REF] Falcao | The image foresting transform: theory, algorithms, and applications[END_REF][START_REF] Roerdink | The watershed transform: Definitions, algorithms, and parallellization strategies[END_REF]) have the problem that the optimum solution may not be unique, for example on plateaus. This 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 [START_REF] Roerdink | The watershed transform: Definitions, algorithms, and parallellization strategies[END_REF], 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][START_REF] Najman | Watershed of a Continuous Function[END_REF].

• 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.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 [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF]. 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.4) are treated as binary terms connected to phantom seeds v F and v B as in (4.8).

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. 4.4. 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. 4.4 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 generally NP-hard if constraints impose more than two different labels. (a) Seeds, (b, c) power watershed result (q = 2).

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

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.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 (4.1) with the value of β hand-optimized to provide the best results independently 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 asymmetrically 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 reassembles undetermined areas and in performing the random walker on 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). 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 sort the weights with a linear sort algorithm.

Computation time

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.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.5 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 a software package provided at http://www.cs.ucl.ac.uk/staff/V.Kolmogorov/software.html and described in [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF]. Our implementation of the shortest path is performed with a Fibonacci heap using double precision weights. For the implementation of Prim's algorithm, we used weights with integer precision, and red black tree as a sorting data structure. Finally, the random walker algorithm was implemented following the multigrid method described in [START_REF] Grady | A lattice-preserving multigrid method for solving the inhomogeneous Poisson equations used in image analysis[END_REF] for 2D image segmentation, and a conjugate gradient descent method for 3D image segmentation.

Qualitative assessment

Unlike most watershed algorithms, the power watershed algorithm (with q = 2) has the property of providing a unique segmentation. Fig. 4.9 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. 4.10 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 [5]. 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 welldefined 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 higher than all of the algorithms except the pure MSF.

Conclusion

In this chapter we clarified, simplified and extended the recent work connecting graph cuts and watershed [4,5]. 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., [START_REF] Lempitsky | Fusionflow: Discrete-continuous optimization for optical flow estimation[END_REF][START_REF] Shen | Stereo matching using random walks[END_REF][START_REF] Szeliski | A comparative study of energy minimization methods for Markov random fields with smoothness-based priors[END_REF]). 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. The next chapter follows that idea, presenting three original applications of power watershed in computer vision.

Chapter 5

Power watershed applied to filtering, stereovision and surface reconstruction

In this chapter, we present three applications of the power watershed method to image restoration1 , stereovision and surface reconstruction from noisy sets of points 2 . This study paves the way for using the power watershed as a useful general-purpose minimizer in many different computer vision contexts.

Anisotropic diffusion

Not all problems in computer vision are convex: For instance, L 0 norm optimization such as seen in compressive sensing is not. In Chapter 4, a novel discrete framework encompassing many known segmentation methods was proposed: power watershed. We are interested in exploring the possibilities of this minimizer to solve problems other than segmentation, in particular with respect to unusual norm optimization. In this section we reformulate the problem of anisotropic diffusion as an L 0 optimization problem, and we show that power watersheds are able to optimize this energy quickly and effectively.

Introduction

The most common assumption about image noise is that the noise exhibits high frequency.

Therefore, a natural choice of image filter is a lowpass filter which dampens power in the highfrequency range. Unfortunately, a lowpass filter has the undesirable effect of blurring object boundaries, which also feature high frequencies. This dilemma was recognized early in the image processing literature and various approaches were proposed to solve it, with the goal of smoothing the image content internal to an object, but preserving image discontinuities across boundaries.

An approach to discontinuity-preserving filtering was proposed by Perona and Malik (PM) [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] who modeled image filtering as an anisotropic diffusion process that smoothed image intensities inside an object while preserving the intensity discontinuity between objects.

Later, the anisotropic diffusion approach was shown by Black et al. [START_REF] Black | Robust anisotropic diffusion[END_REF] to be interpreted as a gradient descent method for optimizing a robust error function model which had been proposed previously. Anisotropic diffusion has been a very successful filtering algorithm, in part because it is easy to implement. However, this algorithm does have the problem of setting two parameters, the robust estimator parameter and the diffusion time. Practical usage of anisotropic diffusion requires a choice between long computation times or blurry boundaries.

Although the power watershed technique was introduced in the context of image segmentation in Chapter 4, the method can be used as an optimization method for some functionals. In this chapter we show that PW is well-suited to address the robust estimator filtering model and therefore provides an alternative to anisotropic diffusion for the optimization of this model. There are several advantages of using PW as compared to anisotropic diffusion in the context of optimizing the robust estimator filtering model. First, we do not require the robust estimator parameter. Second, we preserve discontinuities without blurring images.

Third, the PW method is a very fast optimization procedure. Fourth, no time step needs to be determined (eliminating the risk of a divergent solution). Therefore, our optimization via PW removes the unfortunate compromise presented by anisotropic diffusion between fast optimization and sharp object boundaries.

Other models for discontinuity-preserving filtering are also well-known, including total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], the Mumford-Shah functional [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] and the piecewise constant Mumford-Shah (Chan-Vese) model [START_REF] Vese | A multiphase level set framework for image segmentation using the Mumford and Shah model[END_REF]. The total variation model is convex, allowing for a relatively efficient solution [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF][START_REF] Darbon | Image restoration with discrete constrained total variation part I: Fast and exact optimization[END_REF] while the Mumford-Shah and Chan-Vese models are more difficult to optimize. However, recent work [START_REF] Bresson | Fast global minimization of the active contour/snake model[END_REF][START_REF] Grady | The piecewise smooth Mumford-Shah functional on an arbitrary graph[END_REF] has provided more efficient optimization methods.

Despite this progress in the optimization of other discontinuity-preserving filtering methods, very little work has been done to improve the speed/blur tradeoff that is necessary in the anisotropic diffusion algorithm. Rather than promoting one discontinuity-preserving model over another, we simply show here how the PW may be used to very efficiently optimize the robust estimator model underlying anisotropic diffusion such that we no longer have to choose between fast optimization and blurred boundaries. Our optimization will be achieved at the cost of a few iterations of the efficient PW, which is simple to implement. where the tradeoff between gradient smoothness and data fidelity is governed by λ = 1 dt . Therefore, by setting a fixed time for the anisotropic diffusion, we may view the solution obtained for this time as a steady-state optimization of our second energy functional (5.6) with a corresponding λ parameter. Viewing the time parameter as a loss function for the fidelity of the filtered image with the original, we may freely alter the loss function.

Anisotropic diffusion and L 0 norm

As said in the previous section, it is possible to alter the loss function. For example, we could optimize the fully robust energy

E(x) = σ(Ax) + λσ (x -f ) ,
with gradient given by

dE dx = 2αA ⊤ σ ′ (Ax)Ax + 2αλσ ′ (x -f )x.
If, at any iteration, we fix the values of x inside the robust error function, then we have

dE dx k+1 = A ⊤ σ ′ (Ax k )Ax k+1 + λσ ′ (x k -f )x k+1 .
(5.8)

This energy may be written as the steady-state optimization of the energy functional

E k+1 = x k+1 ⊤ A ⊤ σ ′ (Ax k )Ax k+1 + λx k+1 ⊤ σ ′ (x k -f )x k+1 ,
that may also possibly be written

E k+1 = e ij σ ′ (Ax k ) x k+1 i -x k+1 j 2 + λ v i σ ′ (x k -f ) x k+1 -f 2 .
(5.9) Use PW with y = f to optimize (5.10) to obtain

f x k
x k+1 . k = k + 1; until ||x k+1 -x k || 2 < ǫ
This expression for the energy to compute a minimum step is of a form that may be optimized by the power watershed (PW) if the parameter α → ∞ for g(x).

We recall that the generalized PW energy is given by min

x e ij ∈E w ij p |x i -x j | q + v i ∈V w i p |x i -y i | q (5.10)
where y represents a measured configuration and x represents the target configuration. In this equation, the first term essentially forces x to vary smoothly within an object, and the second term enforces data fidelity. We can see that Eq. (5.9) is of the same form as Eq. 5.10, with q = 2, p = α. In the initialization step of the algorithm, in image filtering applications, we need to build a graph in linking each pixel node with its neighbors, for example its 4 neighbors for 4-connectivity. The values of those nodes will be the x k . Noting N the number of pixels in the image, we add N nodes to this graph in linking each pixel node v i by an edge to the new node v i+N . Those additional nodes are set to constant values of f . A description of the graph on a example is shown in Figure 5.1. The algorithm using power watershed for anisotropic diffusion is summarized in Alg. 8.

Results

We now demonstrate the performance of the power watershed algorithm for anisotropic diffusion in presenting results on synthetic and real images.

A first, the test in Figure 5.2 is performed on a synthetic image corrupted with a Gaussian noise of standard deviation σ 2 = 16. For each result, we may compute the peak signal-tonoise ratio (PSNR) relatively to the original image. Typical values for the PSNR in denoising lie between 20 and 40 dB where higher is better. In comparison with PM algorithm, the power watershed algorithm tends to produce piecewise constant results. The result obtained at Figure 5.2(d) shows a good compromise between noise removing and edge preservation.

However there is a choice to make for PM algorithm ( Here noise and small vessels are both removed, leading to a result which may be used as a first step before segmentation.

result as shown in Figure 5.5. We give some timing with necessary number of iterations to reach convergence at Table 5.1.

Conclusion

The power watershed (PW) algorithm was originally presented as a technique applied to image segmentation. In this section, we showed that PW provides an optimization procedure that allows us to optimize robust error measures which are operated in a particular parameter range.

Black et al. showed how the anisotropic diffusion method could be viewed as optimization of a robust error filtering model. We showed that PW could be used to optimize the same robust error filtering model, leading to an alternative optimization procedure.

An aspect of the PW optimization in this context is that it applies to the optimization of the robust error filtering model when operated in a particular parameter range which effectively models the image as piecewise constant. For cases in which this model was reasonable (such as the synthetic images), the algorithm exhibited strong denoising capabilities. When this denoising algorithm was applied to real images, the effect of imposing a piecewise con- stant model was to effectively quantize images into a small number of grayscale levels (Figure 5.4 and 5.5). This transformation of an image into a piecewise constant form may be helpful as a preprocessing step for segmentation or recognition [START_REF] Chaux | A spatial regularization approach for vector quantization[END_REF].

Stereovision

Dense stereo is a technique for estimating a depth map of a scene using the point of view of several images. In effect, as illustrated in Figure 5.6, the disparity, which is the distance between two corresponding points of the same object in two different views, is a function of the depth. Among the large diversity of two-frame dense stereo techniques [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF], global optimization methods such as graph cuts [START_REF] Woodford | Global stereo reconstruction under second-order smoothness priors[END_REF] and total variations minimization methods [START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF] have been quite successful. We present in this section an example showing the ability of power watershed to perform dense stereo as well.

In order to compute the disparity map from two aligned images, we can build a graph as illustrated in Figure 5.7. The nodes of the graph correspond to the disparity labels x to be recovered.

The nodes marked l 1 , ..., l k correspond to labeled nodes, with values 0, ...k -1 corresponding to distances (in pixels) for matching pixels. The unary weights w i,l j between a nodes v i and a node l j are given by similarity coefficients between blocks centered in the position i in the first image and blocks in position i shifted of l j in the second image. The pairwise

Surface reconstruction

Surface reconstruction from a set of noisy point measurements has been a well studied problem for several decades. Recently, variational and discrete optimization approaches have been applied to solve it. These methods demonstrate good robustness to outliers thanks to a global energy minimization scheme. In this chapter, we use the power watershed algorithm, which is fast and robust to markers placement to produce smooth surfaces. Experiments show that the algorithm compares favorably in terms of speed, memory requirement and accuracy with existing algorithms.

Introduction

This chapter develops a watershed-based algorithm providing a global optimal solution to the surface reconstruction problem from a set of scattered points, under the hypothesis that the surface to recover is a close surface. Surface fitting is a challenging problem when dealing with data containing sparse noise, gaps, and outliers. The initial set of points may be for example acquired by several scans of an object (range scanning). In this context, regularizationbased methods have been shown to be robust when the points lack connectivity, ordering information, and may be contaminated by noise. While there exist numerous explicit surface extraction techniques that estimate the exact positions of surface points, in this work we will focus on implicit surface representation. Implicit surfaces may be represented by level sets (e.g. [START_REF] Zhao | Implicit, nonparametric shape reconstruction from unorganized points using a variational level set method[END_REF]) or binary partitions (e.g. [START_REF] Lempitsky | Global Optimization for Shape Fitting[END_REF]).

Local methods for surface reconstruction including the MPU method [START_REF] Braude | Contour-based surface reconstruction using mpu implicit models[END_REF] are sensitive to noise, as shown in the experiments of [START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF]. Among the recent global approaches, the

Poisson method [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] and FFT-method [START_REF] Kazhdan | Reconstruction of solid models from oriented point sets[END_REF] are more robust to noise, however they require orientation information.

The method of Jalba and Roerdink [START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF] makes it possible to avoid having to estimate orientation information. This is achieved by computing approximations of Coulomb potentials in a grid as an input for the convection method of [START_REF] Zhao | Implicit, nonparametric shape reconstruction from unorganized points using a variational level set method[END_REF]. We propose a different approach that allows us to perform a global optimization of the surface reconstruction problem without requiring the computation of such field as Coulomb potentials.

A generic optimization-based regularization formulation for shape fitting minimizes the total variation functional weighted by the distance function from the set of points P . More generally, given two positive numbers p and q, we consider an object indicator partition u where Ω in is the set of labels inside the surface and Ω out is the set of labels outside the surface.

The weight function w is defined at every point z in a grid as w(z) = d P (z), where d P (z) is the distance map from the points. When p is finite and q = 1, Eq. (5.11) leads to a binary solution u [START_REF] Strang | Maximum flows through a domain[END_REF]. A solution can be deduced in the discrete setting using e.g. the network flow technique [START_REF] Ford | Maximal flow through a network[END_REF], also known as Graph cuts [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF]. Augmenting path max flow implementations are fast and efficient in 2D but memory consuming as the connectivity increases, for instance in 3D. Lempitzky and Boykov [START_REF] Lempitsky | Global Optimization for Shape Fitting[END_REF] have overcome this problem by limiting the size of the search for a solution in the grid while still guaranteeing a global optimum. However their solution is based on restrictive assumptions assuming that the data points are provided with an estimate of the surface orientation. Furthermore, at low resolution, results exhibit metrication artifacts and look blocky, so a high resolution is essential for getting smooth results using the Graph cuts method.

We propose using the power watershed approach in order to provide a way to quickly obtain smooth surfaces at a high resolution without any need to pre-estimate the surface orientation. The power watershed method can be seen as an anisotropic discretization of (5.11) with p → ∞ and q = 2. Although this technique was introduced in the context of image segmentation, as described in previous chapters, it may be employed as an optimization method for various functionals shown in the previous sections 5.1 and 5.2. In the present section we show that the power watershed method is well-suited to address the surface reconstruction problem. The idea of using watersheds for surface reconstruction from a set of points is quite natural. It can be seen as a an extension of the classical "coffee bean" segmentation example, where a watershed is applied on a filtered distance function to separate overlapping convex objects [START_REF] Beucher | Micromorph reference manual, applications and solutions[END_REF].

For further comparison, we also propose to examine a well-known weighted isotropic discretization of (5.11) (for q = 1) known as "Total Variation" (TV). We mention that there exist other TV-based approaches developed to overcome metrication artifacts, e.g. [START_REF] Goldstein | Geometric applications of the split Bregman method: Segmentation and surface reconstruction[END_REF][START_REF] Ye | A fast variational method for surface reconstruction from sets of scattered points[END_REF].

Further comparison with these techniques will be the subject of future work. Finally, we will also show that our algorithm compares favorably with existing surface reconstruction algorithms [7,[START_REF] Braude | Contour-based surface reconstruction using mpu implicit models[END_REF][START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF][START_REF] Kazhdan | Poisson surface reconstruction[END_REF][START_REF] Kazhdan | Reconstruction of solid models from oriented point sets[END_REF].

Method

For solving the surface reconstruction problem, we first place the points cloud onto a regular grid. Our method aims to label the nodes of the grid as an indicator of the object to reconstruct. Since the power watershed is defined on a graph, we begin by casting the surface reconstruction problem formulation in discrete terms.

Given foreground F and background B node values (also called seeds), and p, q two real positive values, the energy presented for binary segmentation in [START_REF] Couprie | Power watersheds: a new image segmentation framework extending graph cuts, random walker and optimal spanning forest[END_REF] is a discretization of (5.11) The definition of the weights for surface reconstruction from a set of points P is based on the construction of a discrete distance map d P , in a grid bounding the set of points.

w ij = min (d P (i), d P (j)) , (5.12)
where d P (i) is the discrete Euclidean distance between the node i and the set of points P .

We recall that exact Euclidean discrete distance map may be obtained in linear time using the algorithm of Hirata [START_REF] Hirata | A unified linear-time algorithm for computing distance maps[END_REF], and that high-quality ordered algorithms also exist [START_REF] Ragnemalm | The euclidean distance transform in arbitrary dimensions[END_REF] . The background seeds may simply correspond to the frame, or bounding box of the lattice. The foreground seeds can be given by the maxima of the distance function that are not connected to the frame. The distance map may be previously filtered, for example using an attribute filter to obtain more robust markers [START_REF] Breen | Attribute openings, thinnings and granulometries[END_REF][START_REF] Najman | Building the component tree in quasi-linear time[END_REF].

As we illustrate in the remainder, the energy defined in (5.10) essentially forces x to remain smooth within the object, while allowing it to vary quickly close to point clusters near the boundary of the object. The data constraints enforce fidelity of x to a specified configuration, taking the values zero and one as the reconstructed object indicator. Observe that the values of x may not necessarily be binary when the value of q is strictly greater than one, which is a positive point for the surface reconstruction problem as we will further explain in this work.

The different values of p and q lead to different algorithms for optimizing the energy.

When the power of the weight, p, is finite, and the exponent q = 1, we recover the Graph cuts energy which can be optimized by a max flow algorithm. When q = 2, we obtain a combinatorial Dirichlet problem also known as the Random walker problem [START_REF] Grady | Random walks for image segmentation[END_REF]. As described in [START_REF] Couprie | Power watersheds: a new image segmentation framework extending graph cuts, random walker and optimal spanning forest[END_REF][START_REF] Couprie | Power Watersheds: A Unifying Graph Based Optimization Framework[END_REF], when the exponent p tends toward infinity, the cut obtained when minimizing the energy is a watershed cut [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF], which has been shown to be equivalent to Maximum Spanning Forests (MSF). Furthermore, an algorithm is presented to optimize the unique watershed that optimizes the energy for q = 2 and p → ∞.

This set of parameters q = 2 and p → ∞ is particularly interesting :

1. The power watershed algorithm has a well-defined behavior in the absence or lack of (a) Three dots in a 4 × 5 lattice. (b) Associated lattice weighted by a distance map from the dots according to (5.12) and squared, with Foreground and Background seeds. Note that only the ordering of the weights counts in the power watershed algorithm, thus the squared distance map can be used directly and produces the same solution as if the weights in (5.12) were used. (c) First steps of the power watershed algorithm to optimize (5.10) in the case q = 2 and p → ∞. Nodes having a maximum weight are merged. (d) A plateau of weight 2 (in green) including different seeded nodes is encountered. The Random walker algorithm is applied to label the nodes on the plateau. (e,f) New plateaus of weight 1 and 0 are encountered, the Random walker algorithm is applied, (g) Final labeling x solution of (5.10). (h) The isocontour is represented in red.

weight information (presence of plateaus). An example is shown at Figure (5.9).

2. The worst-case complexity of the power watershed algorithm in the case p → ∞ is given by the cost of optimizing (5.10) for the given q. In the best-case scenario (all weights have unique values), the power watershed algorithm has the same asymptotic complexity as the algorithm used for a MSF computation (quasi-linear) (See [START_REF] Chazelle | A minimum spanning tree algorithm with inverse-ackermann type complexity[END_REF] for more details). In practical applications where the plateaus have a 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.

Although the PW algorithm is fast, it can be further accelerated for the specific case of weight defined according to a distance map (5.12). Given the foreground and background seeds, we define a narrow band S as the set obtained by thresholding the distance map d P with the the smallest threshold T S such that the connected components are divided between at least an interior (foreground) and an exterior (background). We then compute the PW only on this incomplete distance map, which saves both time and memory. In practice, it is possible to avoid computing an exact distance map on the full grid, using for instance an ordered algorithm propagating from the point cloud [START_REF] Ragnemalm | The euclidean distance transform in arbitrary dimensions[END_REF].

Recall that the exterior seed is connected to the frame of the image, so this is a simple unambiguous connectivity criterion. Applying PW on S is guaranteed to provide the same (unique) global optimizer of the energy as the solution on the full grid, because the connectivity criterion ensures this computation yields a Jordan curve (2D) or surface (3D) separating foreground and background seeds and passing through only already-computed distance values. Performing the PW computation on the full map would not change this result because no distance weight would be lowered, and so no new surface with a lower weight could be found on this full map. Conversely, thresholding the distance map at a lower value than T S would change the result, as any Jordan curve or surface separating foreground and background seeds and computed on this map would necessarily cross some nodes where the distance map had not been computed. Note this does not mean that this thresholding criterion is necessary and sufficient for optimal computation and least memory usage, as an adaptive threshold depending on the local point density could be used instead. However, we leave it as future work to find a better criterion than the simple global threshold.

Results and comparative evaluation

Comparison with Graph cuts and Total Variation

We now demonstrate the performance of the power watershed algorithm for surface reconstruction with respect to two graph-based methods discretizing the energy defined in equation (5.11), namely the weighted total variation (TV) and the Graph cuts (GC) method.

Our first experiment consists of finding a contour fitting sparse and noisy dots in a 2D plane. Figure 5.10 compares the result of TV minimization, the Graph cuts result, a maximum spanning forest result, and the result obtained with power watershed algorithm (q = 2, p → ∞). We observe that all resulting contours are excluding outliers. The Graph cuts results demonstrate that this algorithm is less sensitive to noise, but the contours are blocky because the obtained object indicators are binary. We note that a post-processing step for producing a smooth isosurface from such a binary object reconstruction has recently been proposed by Lempitzky [START_REF] Lempitsky | Surface extraction from binary volumes with higher-order smoothness[END_REF]. The TV contour is the smoothest one, which may be an unwanted effect for rendering details in surfaces. Thus, the TV method requires us to adapt the parameter p in the exponent of the weights to the desired smoothness of the surface. The maximum spanning forest result passes through most points and the contour looks noisy. The power watershed result demonstrates good performance for fitting the dots, while avoiding both blocky contours and noise. This good performance is due to the presence of interwoven plateaus in the distance map. During the execution of the algorithm, the Random walker (RW) algorithm is called several times around the dots, resulting in a smooth output x.

An isocontour or isosurface computation at the 0.5 level is thus providing smooth contours compared to the binary results obtained with Graph cuts or maximum spanning forests.

Figure 5.11 shows surfaces reconstructed from noisy scanned dot sets using the power watershed (PW) algorithm. We used the coordinates of points acquired from scans of several 3D shapes (bunny, Buddha) from the Stanford database available online [2] In our experiments, we embedded those points in 6-connected grids. Quantitative comparisons for the fitting quality are difficult because not all data points are required to be part of the surface.

However, we show in Figure 5.12 that the power watershed method is producing reasonable results even if the point cloud is corrupted by Gaussian noise. We also compared our results to the results obtained using TV minimization and Graph cuts at Figure 5.14. We can observe that the surface obtained with the Graph cuts method is quite blocky. Using the same rendering method to render the output x minimizing the power watershed energy, the power watershed algorithm obtains a smoother surface showing significantly more details.

Figure 5.13 shows that the power watershed method performs well even with a small amount of seeds. In contrast, the total variation method requires a large amount of seeds placed close to the searched surface.

Computation times and memory requirements

In our comparisons, we used the C++ software library of Lemptitzky and Boykov available online which implements the touch-expand algorithm that minimizes the GC energy. For that purpose, the touch-expand algorithm calls a max flow algorithm on a partial graph, a band around the points which is extended when the solution touches the boundary of the band.

We also compare to a fast TV-based solver often called split-Bregman algorithm [START_REF] Goldstein | The split Bregman method for ℓ 1 -regularized problems[END_REF] but identical in principle to the Alternating Direction Method of Multipliers [START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF][START_REF] Setzer | Split Bregman algorithm, douglas-rachford splitting and frame shrinkage[END_REF]. This TV solver is implemented in C and may be called from Matlab. Finally, we also implemented the PW method in C. Computation times of the compared algorithms for the bunny and Buddha data sets are provided on Table 5.2. The three methods have large memory requirements when applied to the full grid. However, the graph cuts and power watershed banded methods are less memory intensive. For the bunny dataset, the touch-expand algorithm optimizing the Graph cuts energy only needs to allocate 3.6% of the full grid size. The solution space for the bunny dataset is much larger for our power watershed algorithm and reaches 33% of the full grid, because some points are widely spaced out. Thus, the touch expand algorithm is faster than the power watershed method. However, considering the ratio time per number of grid node, the power watershed approach is faster. The TV approach has currently no guarantee to produce a global optimum when called on a banded graph, so the TV algorithm was implemented only on a full grid. Table 5.2: Timing experiments (in seconds) for reconstructing the bunny surface using an PC with a 3GHz Intel dual-core processor and 2G of RAM. The memory requirements are given in Mega Octets, and the computation times in seconds. Note that the seeds used are different for the three methods. For TV and PW, the seeds are imposed as hard constraints and are located far from the point cloud (see Figure 5.13). The GC touch expand method uses soft constraints seeds computed from the normals given in input with the point cloud.

We observe that Graph cuts are faster than PW and TV. However, note that the solution is computed on only 3.6% of the full grid, ten times less than PW. On the contrary, the ratio time per voxel shows that the PW could be obtained in a short time using a similar graph reduction.

Comparison with some other approaches

Following the quantitative information given in [START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF], we compare the performances of power watershed for the Stanford bunny reconstruction with different methods [7,[START_REF] Braude | Contour-based surface reconstruction using mpu implicit models[END_REF][START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF][START_REF] Kazhdan | Poisson surface reconstruction[END_REF][START_REF] Kazhdan | Reconstruction of solid models from oriented point sets[END_REF].

All these methods, including the power watershed, reconstruct a surface close to the Bunny set of dots, with an error comprised between 2 × 10 -4 and 6 × 10 -4 . In terms of computation time, the Power Crust method [7] is about 10 times slower, and the Poisson method [START_REF] Kazhdan | Poisson surface reconstruction[END_REF] 3 times slower than our PW approach. Although the FFT [START_REF] Kazhdan | Reconstruction of solid models from oriented point sets[END_REF], MPU [START_REF] Braude | Contour-based surface reconstruction using mpu implicit models[END_REF] and Hoppe et al [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF] methods are fast, the FFT method suffers from large memory requirements limiting the grid resolution, and Hoppe et al and MPU methods produces artifacts in the presence of noise (See [START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF]). The method of Jalba and Roerdink [START_REF] Jalba | Efficient surface reconstruction using generalized coulomb potentials[END_REF], based on Coulomb potentials, uses 4 times less memory than our power watershed method, but is 5 to 10 times slower on a CPU and is still slower using a GPU. Furthermore our PW is more flexible in the choice of the markers. In our experiments, the amount of markers is not very large as shown in Figure 5.13, but a strategy using larger seeds could be employed to reduce the size of the solution space and the computation time.

Conclusion

The power watershed method can be used to efficiently produce surfaces fitting noisy measurements. Contrary to standard watershed algorithms and the Graph cuts approach, the unique solution provided by the power watershed is not binary, resulting in the reconstruction of both smooth and detailed surfaces. In addition, this method is fast and not limited by large memory requirements, when using a restriction of the solution space. Finally, in com-parison with other methods, our power watershed is robust to seed placement, and requires fewer parameters to be set. In practice, when using the power watershed method, closefitting markers are not as mandatory as in other methods. Future work will follow several directions: memory and computation times improvements are still achievable, in particular by improving the touch-expand idea used in the Graph cuts optimization and adapting it for power watershed. The power watershed energy could also be modified to add some priors to the reconstructed surface, such as local orientation. Finally, we hope to demonstrate the efficiency of the power watershed technique for solving related problems such as multiview reconstruction.

We used the energy minimization aspect of PW to expand the traditional use of watersheds from a segmentation algorithm to image filtering, stereovision, and surface reconstruction algorithms. Although our algorithm optimizes the filtering objective or disparities as a realvalued optimization problem, we showed that we still retain the speed of the watershed algorithm. However, the limit of PW as an energy optimization strategy remains unclear -What energy functions can be minimized via power watersheds? Future work will address this issue and continue to demonstrate applications which can benefit from the speed of power watershed.

Chapter 6

Conclusion

In this thesis we have proposed several graph-based variational approaches and have provided theoretical links of significant importance between existing and new methods. We summarize in Table 6.1 the diverse regularization terms studied and developed in this work. It is not exhaustive but it illustrates the different energies from a unified point of view. The methods we developed were successfully applied to such diverse problems as image segmentation, restoration, stereovision and surface reconstruction.

Contributions

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

We have also compared our algorithm to known weighted combinatorial TV minimization methods. We have shown that our results are generally better for segmentation than TV-based methods, and that our implementation on CPU can be 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 the fastest TV-based methods.

Method Energy formulation Details

Max flow / Graph cuts arg min which we term "Power watershed". The algorithm produced for computing power watershed is a novel, globally optimal optimization method. Our experiments showed that the case q = 2 retains the speed of efficient watershed algorithms while providing improved segmentations. Furthermore, this case exhibits the interesting properties of uniqueness of the solution, and possible application for fast, exact global optimization of multi-label problems.

Placing the maximum spanning forest for watershed in the same energy minimization framework than state-of-the-art optimization algorithms -used in a large variety of applications in computer vision -opens the way to new use of watershed within the computer vision field. Indeed, we have shown that it may be used beyond its traditional domain of image segmentation. Specifically, we have demonstrated how to use power watershed as a solver for image filtering optimizing an L 0 norm energy. We also presented an example of stereoreconstruction using power watershed, and performed fast and smooth surface reconstruction from a noisy cloud of 3D points.

Future work

The work presented in this thesis offers many avenues for further investigation.

Flow based methods -The current optimization of the CCMF problem -performed by an interior point method -is faster than the existing continuous methods. However, a drawback of our current implementation of the interior point method is that the resolutions of linear systems can be expensive for 3D tasks. Future improvements could involve the implementation of preconditioned conjugate gradient for the linear system resolution, or the use of fast first order methods. In the context of restoration using DCTV, diverse formulations of the constraints on the projection variable remain to be explored.

Power watersheds -Future work 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] Arbeláez | A metric approach to vector-valued image segmentation[END_REF][START_REF] Guigues | Scale-sets image analysis[END_REF]163,[START_REF] Najman | Ultrametric watersheds: a bijection theorem for hierarchical edgesegmentation[END_REF][START_REF] Najman | Geodesic Saliency of Watershed Contours and Hierarchical Segmentation[END_REF] look like an interesting topic that can take advantage of the power watershed uniqueness. A second direction for future work will also include the characterization of the class of energy that may possibly be optimized in this framework. For example, the possibility of using negative weights could be explored, with potential applications in texture analysis, or segmentation with a curvature constraint [START_REF] El-Zehiry | Fast global optimization of curvature[END_REF].

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. The efficiency of the power watershed for large datasets is suitable to many problems in vision (and beyond), such as stereo-vision, inpainting, video segmentation, optical flow, and more.

Finally, further developments may be achieved by using extra information to obtain better results. For example shape constrained segmentation [START_REF] Veksler | Star shape prior for graph-cut image segmentation[END_REF], co-segmentation methods [START_REF] Hochbaum | An efficient algorithm for co-segmentation[END_REF],

and segmentation of moving objects. Extra information coming from machine learning could bring benefits. For example, the convolutional networks of Lecun et al. [START_REF] Yann Lecun | Convolutional networks and applications in vision[END_REF] simulate the learning of a neuronal network by mixing thousands of simple convolution operations, leading to and effective real-time object recognition and classification. However, there are still possible improvements for these methods, which provide only rough estimates of objects contours. It will be highly interesting to study the potential benefit of graph-based approaches for image recognition in the learning process.

Taking the sum of right hand side and left hand side,

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

F c ⊤ F = c ⊤ F * = 2λ * ⊤ g 2 .

A.1.2 Weak duality of CCMF and CTV

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

In the combinatorial setting the weak duality property is given by

||F || ≤ g ⇒ F ⊤ (Au) ≤ g ⊤ ||Au||. (9)
The next property shows that the norm ||F || = |A ⊤ |F 2 verifies weak duality.

Property 6. Let G be a transport graph, F a flow in G verifying the capacity constraint |A ⊤ |F 2 ≤ g. Let u be a vector of R n (defined on nodes of G). Then

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

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

F ⊤ Au ≤ (|A ⊤ |F 2 ) ⊤ |A ⊤ |(Au) 2 .
Using summation notation, then by the Cauchy-Schwartz inequality,

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

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

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

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

a i =           
-1 if v i belongs to the sink, 1 if v i belongs to the source, 0 otherwise. [START_REF] Arbeláez | A metric approach to vector-valued image segmentation[END_REF] We can provide the following equivalent formulation of the CCMF problem (2.3) using an

The problem (2.6) may also be written as min

F -c ⊤ F + i C (F ) + i D (F ). ( 15
)
In order to employ the parallel proximal algorithm of Section 1.

C k = {F ∈ R m | ∀i ∈ S k |A ⊤ | i F 2 ≤ g 2 i , } (16)
C 1 = {F ∈ R m | (F 2 1 + F 2 2 + F 2 4 ) 1 2 ≤ g 2 2 (F 2 3 + F 2 6 ) 1 2 ≤ g 2 4 (F 2 5 + F 2 7 ) 1 2 ≤ g 2 6 } C 2 = {F ∈ R m | (F 2 1 + F 2 3 ) 1 2 ≤ g 2 1 (F 2 2 + F 2 5 )
where s denotes the source node, and t the sink node. Similarly to the case of C, the projections on convex sets D i are given by simple projections onto hyperplanes.

s. t. λ ≥ 0.

The derivatives with respect to λ and ν are given by

∇ λ (E) = - 1 4 |A| ⊤ (c + Aν) • /(|A|λ) 2 + g 2 , ∇ ν (E) = 1 2 A ⊤ (c + Aν) • /(|A|λ)
A conjugate gradient descent has been tested for the optimization of [START_REF] Beucher | Micromorph reference manual, applications and solutions[END_REF]. The constraint of positivity for λ was enforced by setting to neglectable but positive values the negative values of λ at each step. A second strategy was to employ a method alternating between minimizations in function of λ and minimization in function of ν. Both approaches were converging toward the optimal solution until a certain point but the precision required was rarely reached.

Further possible tests would include implementations of Beck and Teboule's FISTA algorithm, and projected gradient method (i.e [START_REF] Bertsekas | Projected newton methods for optimization problems with simple constraints[END_REF]).

A.2 Power watershed

A.2.1 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.4) is a Voronoi diagram of the seeds.

A proof showing that solving the minimization problem (4.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.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.4)

q e ij ∈E w p ij |x i -x j | q = q e ij ∈E w ij p q |x i -x j | q (20)
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

The minimization problem can be written as

) 22
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).

A.2.2 Illustration of Theorem 4.3.1

In theorem 4.3.1, the condition for seeds to be the maxima of the weight function is necessary as illustrated in the following Figure 3. 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 q-cut (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. 4.3.1. Note that the threshold of lim p→∞ x is an MSF cut, as stated later in Prop. 5.

A.2.3 Proof of Property 4

If q is a real number such that 1 < q < ∞, then the solution x to problem (4.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 non-decreasing. The function f : x → w ⊤ 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 p⊤ |Ax| q . ( 23
)
As h is a non-decreasing convex function, and g is convex, h • g is convex. As f is a nondecreasing 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.

A.2.4 Proof of Theorem 4.4.1

The energy of the solution computed by the power watershed algorithm converges to the energy which minimizes E p,q when p → ∞.

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 ≤ δ. (24)
Proof.

E p,q (x)

w M p = e M |x i -x j | q + e ij =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 + e ij =e M w ij w M p | xi -xj | q - e ij =e M w ij w M p |x * i -x * j | q . ( 25
)
The first part of ( 25) is bounded by 0,

e M | xi -xj | q - e M |x * i -x * j | q ≤ 0, (26)
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 7 are equal, leading to a sum equal to 0. Else (if there are plateaus with different labels), e M |x i -x j | q subject to the boundary constraints is minimized on the plateaus, so the solution is optimal.

The last part of ( 25) is also negative, -

e ij =e M w ij w M p |x * i -x * j | q ≤ 0. ( 27
)
x ⋆ to (4.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) [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in N-D images[END_REF] in which p(x) models how the parameters of the Bernoulli variables vary spatially. The spatial smoothness prior is parameterized by

v i ∈V p(s i |x i ) v i ∈V p(I i |s i ),
p(x) ∝ exp   -λ e ij ∈E w p ij |x i -x j | q   , (33)
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

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   - v i ∈V w p i0 |x i -0| q - v i ∈V w p i1 |x i -1| q   , (35)
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 (34) is equal to our energy minimization problem from (4.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 [8,[START_REF] Kaufhold | Energy Formulations of Medical Image Segmentations[END_REF][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] Krajsek | Diffusion filtering without parameter tuning: Models and inference tools[END_REF][START_REF] Michael | Image statistics and anisotropic diffusion[END_REF] and MRFs requiring continuous-valued estimations have appeared in both early work on computer vision [START_REF] Bouman | A generalized Gaussian image model for edge-preserving MAP estimation[END_REF][START_REF] Geman | Constrained restoration and the discovery of discontinuities[END_REF][START_REF] Geman | Statistical methods for tomographic image reconstruction[END_REF][START_REF] Levitan | A maximum a posteriori probability expectation maximization algorithm for image reconstruction in emission tomography[END_REF] 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].

A.2.6 Using mathematical morphology for an efficient preprocessing step

One difficulty in Algorithm 7 is dealing with the set of merged nodes. More precisely, when solving (4.7) 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) ≥ λ}.

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 8. 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 7.

The components of w, ordered by the inclusion relation, form a tree called the maxtree [START_REF] Salembier | Anti-extensive connected operators for image and sequence processing[END_REF] or the component tree [START_REF] Breen | Attribute openings, thinnings, and granulometries[END_REF][START_REF] Guillataud | Contribution l'analyse dendroniques des images[END_REF][START_REF] Hanusse | Sémantique des images par analyse dendronique[END_REF]. Several efficient algorithms exist to compute the List of Figures

  vertices of a graph. E : Set of edges of a graph. m : Number of edges in a graph. n : Number of vertices in a graph. A : Incidence node-edge matrix of a graph. g : Vector of nodes weights. w : Vector of edge weights. s : Index of a source node. t : Index of a sink node. c : Vector indicator of the sink-source edge in a transport graph. |A| : Matrix composed of the component-wise absolute values of the elements of A. |A| i : The i th row of |A|: sum operator through neighbors of node i. Let u, v be real vectors of R N . u • v : Hadamar product (component-wise) of u by v. u p : Component-wise power operator of u. u ⊤ v : Scalar product of u by v. u • /v : Component-wise division of u by v. Let C be a nonempty closed set, C ⊂ R N . ı C (u) : Indicator function of C. σ C (u) : Support function of C. P C u : Projection of u on C.

Figure 1 . 1 :

 11 Figure 1.1: Clustering and restoration problems, in images and arbitrary graphs.

= p 1 p 2 p 3 p 4 e 1 1 Figure 1 . 2 :

 1112 Figure 1.2: A graph and its incidence matrix A ∈ R m×n .

  of data fidelity term in an image denoising context could beD(x) = Ω (x(z) -f (z)) 2 dz,appropriate in the case where f is an image corrupted with Gaussian noise.

Figure 1 . 3 :

 13 Figure 1.3: Image denoising example by TV minimization on a non-local graph.

Figure

  Figure from X. Bresson [44].

FigureFigure 1 . 4 :

 14 Figure from C. Rother[START_REF] Rother | GrabCut" -Interactive foreground extraction using iterated graph cuts[END_REF] V. Kolmogorov and R.Zabih[START_REF] Kolmogorov | Computing visual correspondence with occlusions via graph cuts[END_REF]

Figure 1 . 5 :

 15 Figure 1.5: Analogies and illustration of results obtained optimizing the combinatorial Dirichlet problem (a) A graph with two constraints on x: one value set to '0', and another to '1'. Given those boundary constraints, the labeling x optimizing (1.5) is shown on the nodes. (b) Analogy with an electrical network resistance: electric voltage of all nodes are given by the solution to (1.5) in (a). (c) Analogy with random walks: in (a), each node value x corresponds to the probability for a random walker to reach the node labeled '1' before the node labeled '0'. (d) Example of constraints (foreground and background seeds) for seeded image segmentation. (e) Obtained result x optimizing (1.5).

Figure 1 . 6 :

 16 Figure 1.6: Example of shortest path computations for image segmentation of thin, tubular objects. (a,b) Figure from F. Benmansour [17], (c) Example of shortest path forest for general purpose image segmentation.

Figure 1 . 7 :

 17 Figure 1.7: Illustration of two Maximum Spanning Forest (MSF) algorithms behavior, and segmentation result on a real image.

2 .

 2 Compute the primal-dual search direction ∆ y such that M ∆ y = r. 3. Determine a step length s > 0 and set y = y + s∆ y . (y = [x, λ, ν] ⊤ ) Until ||r p || 2 ≤ ǫ, ||r d || 2 ≤ ǫ, and η ≤ ǫ.

. 10 )

 10 For example, if ϕ = ι C , prox ι C = P C , where ι C denotes the indicator function over a closed convex set C and P C denotes the projection onto C. More formally, given x ∈ R N , and a nonempty closed convex set C ⊂ R N ,

Figure 1 . 8 :

 18 Figure 1.8: Figure illustrating the behavior of the Parallel Proximal Algorithm for minimizing a sum of two convex functions.

Figure 2 . 1 :

 21 Figure 2.1: Example of metrication (blockiness) artifacts in cells segmentation obtained with classic max flow, and solution obtained optimising the continuous max flow problem.

Figure 2 . 2 :

 22 Figure 2.2: Different discretization of metric and different type of cuts obtained with classical max-flow and with the combinatorial continuous max-flow (CCMF) formulations.

  ν (d) Threshold ν at .5

Figure 2 . 3 :

 23 Figure 2.3: Solution to the dual of the CCMF problem in a seeded segmentation example.

Algorithm 4 :

 4 Primal dual interior point algorithm optimizing CCMF Data:

Figure 2 . 4 :

 24 Figure 2.4: Comparison of brain segmentation results using GC, AT-CMF, and CCMF.(a) Original image with foreground and background seeds. (b,c,d) Segmentation obtained with (b) graph cuts (GC), (c) AT-CMF, threshold of P (obtained after 10000 iterations), (d) CCMF, threshold of ν (15 iterations).

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

Figure 2 . 5 :

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

Figure 2 . 6 :

 26 Figure 2.6: Comparison of perimeters computed analytically and with Graph cuts and CCMF. To estimate the boundary length with CCMF, segmentations have been produced using the input images in (b) with the foreground and background seeds that appears in (a). (c) : plot showing the linear relation between the true perimeters and the energies obtained with CCMF. The symbols used to represent the dots correspond to the objects segmented in the images. A similar relation of proportionality is obtained with CP-TV. (d) : value of the cut obtained for a diagonal and a vertical line of approximately 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.

Figure 2 . 7 :

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

Figure 2 . 8 :

 28 Figure 2.8: GC and CCMF segmentation results on a image from the GrabCut database.

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

Figure 2 . 9 :

 29 Figure 2.9: Example of graph construction to incorporate unary terms for unsupervised segmentation.

  2.10). Another strategy would be to compute the solution in the dual graph, that is in the graph where nodes are replaced by edges and (weighted) edges are replaced by (weighted) nodes. We have shown in an example the ability of CCMF for classification, a task that can not be treated with finite element-based methods such as AT-CMF. Evaluating the performances of CCMF for classification may be a topic for future research.

  Figure 2.10: Zachary's karate club network and CCMF classification result.

Figure 2 . 11 :

 211 Figure 2.11: Liver segmentation in 3D by CCMF.

• 2 =

 2 a • a. The constraint given in the convex set of (3.11) also appears in the Combinatorial Continuous Maximum Flow (CCMF) problem of Chapter 2,[START_REF] Couprie | Combinatorial Continuous Maximum Flow[END_REF]. The problem studied in the present work may be seen as an extension of the CCMF problem (applied to clustering problems in graphs, such as image segmentation)

Algorithm 5 :

 5 Algorithm 5, this requires computing in parallel projections onto each set C r with r ∈ {1, . . . , s}, which are defined ∀F ∈ R m as P Cr (F ) = arg min Φ∈Cr Φ -F . Note that the convergence of the sequence (F k ) k generated by this algorithm to a solution F of (3.17) is guaranteed, which allows us to deduce a solution to(3.16) by using Relation(3.18). Parallel proximal algorithm solving(3.19) Fix γ > 0 and ν ∈]0, 2[. Set k = 0. Choose y 1,0 = y 2,0 = . . . = y s+1,0 ∈ R m and F 0 . repeat for r = 1, . . . , s + 1 do in parallel

Figure 3 . 2 :

 32 Figure 3.2: Coloring a graph using the DSAT algorithm [43]. (a) Sort nodes by decreasing degree. Initial labeling DSAT equal to node degrees. (b) Choose a node with maximum DSAT. Assign color number 1. Update DSAT: For each non-colored node v, if at least one neighbor of v is colored, DSAT(v) equals the number of different colors used in the neighborhood of v. (c) Choose a non-colored node of maximum DSAT. In case of conflict, choose a node of maximum degree. Assign the smallest possible color. Update DSAT. (d,e,f) Repeat the procedure in (c) until all nodes are colored.

  ation of the algorithm. M+SFBF belongs to the rich class of primal-dual proximal methods which also include the methods by Chen and Teboulle [58], Chambolle and Pock [54] and Esser et al. [90]. Property 3. Problem (3.16) admits the equivalent formulation

Figure 3 . 3 :

 33 Figure 3.3: Denoising an image with constant levels of additive Gaussian noise.(a) Noisy image, (b) TV [50] λ = 0.06, (c) Weighted TV on a 4-connected graph [89] , weights given by the gradient of the image, (d) PPXA norm ℓ ∞ , (e) PPXA norm ℓ 2 (λ = 0.22, ǫ = 0.3).

Figure 3 . 4 :

 34 Figure 3.4: Denoising the 'man' image (zoom) corrupted with Gaussian noise of variance σ 2 = 10.

  with different values of noise variance. Result examples are shown in Figure 3.4. Those experiments show that DCTV leads to improved results when the variance of the noise is lower than 50. Visually, DCTV results are sharper and feature better contrast than the weighted TV results. This may explain the slight degradation of performance in presence of heavy noise (variance ranging from 50 to 100). Geometrically, the improvement of DCTV over the weighted TV can be interpreted in the following way: the convex set we use for projection adapts itself to the local neighborhood, and this avoids oversmoothing as a result. We show, in the experiment of Figure 3.5, an example of the contrast preservation properties of DCTV. In terms of computation time, DCTV is competitive with the most efficient weighted TV algorithm. Denoising the 512 × 512 Lena image corrupted with Gaussian noise (σ 2 = 15) requires 0.38 seconds for split Bregman, versus 0.7 seconds for PPXA to converge on an Intel

Figure 3 . 5 :

 35 Figure 3.5: Synthetic denoising example demonstrating improved contrast preservation using DCTV.

Figure 3 .

 3 6 presents a non local DCTV result obtained using this strategy compared with a non local TV result. About 15 iterations are necessary for M+SFBF to produce this result, taking about two seconds using a C parallel implementation on the same 2.5GHz 8-core Xeon system used previously. The projections on the different (9 in the example of Figure 3.6) convex sets are performed in parallel. In comparison, the split Bregman algorithm used for optimizing non local TV only takes four iterations to converge in about one second.

  (a) Original image (b) Noisy SNR=22.6dB (c) Non-local graph (Figure P. Coupé, [68]) (d) Unweighted TV SNR=25.4dB (e) Nonlocal TV SNR=26.1dB (f) Nonlocal DCTV SNR=26.1dB

Figure 3 . 6 :

 36 Figure 3.6: Local and non-local denoising.

Figure 3 . 7 :

 37 Figure 3.7: Comparison of convergence speed of the PPXA (blue) and M+SFBF (red) algorithms for denoising an image.

  (a) using a regular, 4-connected lattice, (b) using a non-local graph.

Figure 3 . 8 :

 38 Figure 3.8: Denoising and deblurring an MRI image corrupted with synthetic Gaussian 5 × 5 blur and Gaussian noise (σ 2 = 10).

  1 , . . . , λ N ) ∈]0, +∞[ N . We present in Figure 3.9 an example of fusion between a blurry image and a noisy image in a 4-connected lattice (N = 2 and H 2 = I).

Figure 3 . 9 :

 39 Figure 3.9: Image fusion from an image corrupted with Gaussian noise with variance σ 2 2 = 20 (b) and a blurry image (c) (uniform blur kernel of size 5 × 5, with additive Gaussian noise with variance σ 2 1 = 1). (d) DCTV result with λ 1 = 0.17 and λ 2 = 0.24.

Figure 3 . 10 :

 310 Figure 3.10: Example of mesh denoising using DCTV on the spacial coordinates of the nodes.

  ) where b X , b Y , b Z represent vectors of additive noise with unit magnitude variance.

Figure 3 .

 3 Figure 3.11, we used the Graph Analysis Toolbox [116] that contains implementations of several filtering techniques in arbitrary graphs. From a graph representing the spacial resolution of the Bottlenosed dolphin visual system, the image (c) represents image resampled with the dolphin's mesh, corresponding to the input f of our method. It is then straightforward to run the M+SFBF algorithm, using the incidence matrix of the graph represented in (b).

Figure 3 . 11 :

 311 Figure 3.11: Filtering image data on a biologically sampled image [116]. Noise of variance σ 2 = 10 was added to the resampled values of the image (c) to produce (d).

Definition 4 . 3 . 1 .

 431 (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 segmentation result defined in Eq 4.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

Figure 4 . 1 :

 41 Figure 4.1: Illustration of different steps in the proof of Thm. 1 for q = 2.

Fig. 4 . 1 ,

 41 where F and B are displayed in (a), and two possible extensions in bold in (b) and (c), with their corresponding labels.

Theorem 4 . 3 . 1 .

 431 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 ,

  Figures 4.1(c) and 4.1(d)).

Figure 4 . 2 :

 42 Figure 4.2: Example of behavior 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.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.

Figure 4 . 3 :

 43 Figure 4.3: Illustration of progressive convergence of the random walker result to the power watershed result as p → ∞, using q = 2. Top row: Segmentation results obtained by random walker with weights at the power p = 1, p = 8, p = 25 and finally by the power watershed algorithm. Bottom row: Corresponding potentials for p = 1, p = 8, p = 25 and the input seeds.

. 10 )

 10 with boundary conditions given by already labeled nodes. It is exactly what is performed by Algorithm 7 in the "If" part.

Figure 4 . 4 :

 44 Figure 4.4: Unseeded segmentation using unary terms.(a) Example of the unseeded segmentation of a 3 × 3 image computed with a maximum spanning forest (watershed). (b) An image of blood cells, (b) graph cuts result(c) Maximum Spanning Forest (watershed) result.

Figure 4 . 5 :

 45 Figure 4.5: Example segmentations with more than two labels.

11 )

 11 An example of result is shown on Fig.4.5.

Figure 4 . 6 :

 46 Figure 4.6: Example segmentations using the provided (left images) and skeletonized (right 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).

Figure 4 . 7 :

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

Figure 4 . 8 :

 48 Figure 4.8: 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.

Figure 4 . 9 :

 49 Figure 4.9: Segmentation of an artificial image.(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.

Figure 4 . 10 :

 410 Figure 4.10: 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.

function in ( 5 . 3 )

 53 is given by the Welsch function σ(z) = 1 -exp -αz 2 .Black et al. explored the effect of different robust estimator forms of σ(z). However, a common characteristic of all robust estimators is that they exhibit linear or quadratic growth near z = 0, but provide a nearly constant output as z → ∞. The transition between quadratic growth and constant growth is controlled by the parameter α. If α → ∞, then the robust estimator ultimately achieves its constant growth phase for any z = 0. In this way, the σ(z) function may be viewed as a differentiable approximation to the ||•|| 0 norm in which ||0|| 0 = 0 and ||z|| 0 = 1 for any z = 0. Consequently, the anisotropic diffusion filtering algorithm may be viewed as gradient descent on an energy functional which is an approximation to the || • || 0 norm.Previous filtering algorithms that explicitly formulated filtering from the standpoint of optimizing the || • || 0 norm took the formE(x) = σ(Ax) + λh(x, f ),(5.6)in which f represents the intensities of the input unfiltered image, h(x, f ) represents a loss function enforcing data fidelity, and λ is a free parameter. When h(x, f) = ||x -f || 2 2, the gradient of (5.6) becomes dE dx = A ⊤ σ ′ (Ax)Ax + λ (f -x), which reaches a stable point whenA ⊤ σ ′ (Ax)A + λI x = λf.(5.7)Since (5.7) may be viewed as a backward Euler solution for the anisotropic diffusion equation (5.2) when λ = 1 dt , then the anisotropic diffusion algorithm may be seen as the optimizing a robust estimator of the image gradient balanced against a loss function of the form ||x -f || 2 2

Figure 5 . 1 :Algorithm 8 :

 518 Figure 5.1: Graph necessary to build for applying the power watershed diffusion algorithm to a 3 × 3 image.The unary weights in blue and pairwise weights in red are updated at each step.

  Therefore, we suggest to employ the filtering algorithm given in Algo. 8. This technique can be initialized by the image f , or, in the case of very noisy images by a smoother version of the image, for instance by the application of a Gaussian or a median filter. (a) Noisy image, PSNR = 24.24dB (b) PM, PSNR = 34.03dB (c) PM, PSNR = 30.46dB (d) PW, PSNR = 31.54dB

Figure 5 . 2 :

 52 Figure 5.2: Comparison of Perona-Malik(PM), and power watershed(PW) algorithms for denoising a synthetic image. (b) PM used with 80 iterations α = 0.0015, leading to a good PSNR but with a few remaining isolated noisy pixels. (c) PM, best compromise found for this image to remove the isolated pixels with 50 iterations and α = 0.0005. (d) Any median filtered image as initialization and λ = 0.975 results in a better PSNR while removing isolated noisy pixels.

Figure 5 . 3 :

 53 Figure 5.3: Example of piecewise smooth image denoising using the power watershed algorithm. Perona-Malik PW Fig. 5.2, 104 × 100 Nb iter. 50 80 5 Time (s) 0.19 0.30 0.17

Fig. 5 . 5 Table 5 . 1 :

 5551 Fig. 5.5, 250 × 300 Nb iter. 50 80 6 Time (s) 1.38 2.14 1.78 Fig. 5.4, 299 × 364 Nb iter. 50 80 6 Time (s) 1.95 3.08 2.43

Figure 5 .

 5 2(b) and Figure 5.2(c)) between complete denoising and good contrast edge restoration.Examples on real images in Figure5.4 and Figure5.5 show that the power watershed algorithm may be useful as a filtering step before segmentation. In terms of computation time, each iteration of the power watershed algorithm for anisotropic diffusion operates in quasi-linear time in practice. We noticed empirically that the convergence is very fast (fewer than 10 iterations). The result obtained at the first iteration is generally close to the final

Figure 5 . 4 :

 54 Figure 5.4: Filtering of a liver image by power watershed.

Figure 5 . 5 :

 55 Figure 5.5: Segmentation of an image of concrete filtered using PW. The segmentation in three classes (bubbles in black, stones in light grey) is useful for the study of the material's mechanical properties. The segmentation is obtained by two thresholds of the filtered image by PW.

Figure 5 . 6 :

 56 Figure 5.6: Illustration of the relationship between disparity between two views and scene depth estimation.

Figure

  Figure from the IMEC's website (http://www.imec.be)

  ) p |∇u(z)| q dz subject to u(z) = 0 ∀z ∈ Ω in ,and u(z) = 1 ∀z ∈ Ω out ,(5.11)

Figure 5 . 9 :

 59 Figure 5.9: Illustration of surface extraction from dots using PW on a small graph.

Figure 5 . 10 :

 510 Figure 5.10: Comparison of surface reconstruction from a set of points in 2D, using the total variation method, Graph cuts, a maximum spanning forest algorithm and finally the power watershed algorithm.

Figure 5 . 11 :

 511 Figure 5.11: Power watershed results obtained from points clouds.(a) Total size of scans : 362272 points, Grid size : 234 × 297 × 301. (b) Number of scans used : 341072 points, Grid size : 275 × 276 × 668 (c) Total size of scans : 2748318 points, Grid size 382 × 270 × 171.

Figure 5 . 13 :

 513 Figure 5.13: Robustness to seed quantity. This figure shows slices of solutions obtained for reconstructing the bunny surface, obtained on a grid of size 251 × 248 × 195, using different seeding strategies. The interior seed is colored in red, and the background seed in blue, and resulting isolines in green. (a,b,c) : Results obtained with TV minimization, and (d) Result obtained with Power Watershed.

Figure 5 . 14 :

 514 Figure 5.14: Comparison of TV, GC, and PW resulting surfaces from the Stanford bunny dataset. Grid size : 234 × 297 × 301. (a) Total variation minimization result, (b) Graph cuts result, (c) power watershed result. Isosurfaces at 0.5 have been extracted on all results, and were downsampled by 2 to render the surfaces.

From a theoretical view

  point, 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.Dual Constrained TV-based regularization for image restoration -Inspired from the CCMF formulation, the energy formulation has been adapted and extended to multi-label problems. Constraining the dual projection variable in the primal dual total variation problem provides a more flexible way of penalization of variations. Fast parallel proximal optimization tools have been tested and adapted for the optimization of this problem. Experimental results show an improvement of the approach in term of restoration quality and contrast preservation.

  Power watershed: A unifying graph-based optimization framework -We have introduced a framework that generalizes several state-of-the-art graph-based segmentation algorithms, namely graph cuts, random walker, shortest paths, and watershed. This generalization allowed us to exhibit a new family of watershed algorithms using different exponents,

3 . 3 . 1 ,

 331 we need to compute the proximal operator of ı C and ı D , corresponding to P C and P D by definition. Projecting directly on C and D is analytically intractable so we split the two convex sets into several parts. We decompose ı C (F ) into the sum of several indicator function on convex sets C k , k = 1...p. Given a decomposition of the set of nodes V = ∩ p k=1 S k , we define

Figure 1 : 2 For a 4 - 3 D

 1243 Figure 1: Decomposition of C = C 1 ∩ C 2

Figure 2 :

 2 Figure 2: Segmentation of a 256 by 256 image using PPXA.(c) Saturated vertices after 130 seconds (4000 iterations) (d) Saturated vertices after about one hour (120000 iterations). In this case, the obtained map is precise enough to obtain, by propagation of the marker, the segmentation in (e).

  |x i -x j | q = lim q→∞ max e ij ∈E w ij p q |x i -x j |.

q

  |x i -x j |, s.t. x(F ) = 1, x(B) = 0.(

Figure 3 :

 3 Figure 3: Example of seeded graph where the q-cut does not correspond to an MSF cut.

x

  p(x)p(I|x) = arg max p(x) s p(I|s)p(s|x).

1. 1

 1 Clustering and restoration problems, in images and arbitrary graphs. . . . . . 1.2 A graph and its incidence matrix A ∈ R m×n . . . . . . . . . . . . . . . . . . . 1.3 Image denoising example by TV minimization on a non-local graph. . . . . . 1.4 A transport graph with its resulting max-flow/min-cut solution, and illustration of applications in image segmentation and stereovision. . . . . . . . . . . 1.5 Analogies and illustration of results obtained optimizing the combinatorial Dirichlet problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Example of shortest path computations for image segmentation of thin, tubular objects. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Illustration of two Maximum Spanning Forest (MSF) algorithms behavior, and segmentation result on a real image. . . . . . . . . . . . . . . . . . . . . . . . 1.8 Figure illustrating the behavior of the Parallel Proximal Algorithm for minimizing a sum of two convex functions. . . . . . . . . . . . . . . . . . . . . . . 2.1 Example of metrication (blockiness) artifacts in cells segmentation obtained with classic max flow, and solution obtained optimising the continuous max flow problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Different discretization of metric and different type of cuts obtained with classical max-flow and with the combinatorial continuous max-flow (CCMF) formulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .2.3 Solution to the dual of the CCMF problem in a seeded segmentation example. 2.4 Comparison of brain segmentation results using GC, AT-CMF, and CCMF. .

2. 5

 5 Comparison of AT-CMF and CCMF results on the catenoid test problem. . . 2.6 Comparison of perimeters computed analytically and with Graph cuts and CCMF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Segmentation of an artificial image with AT-CMF (top row) and CCMF (bottom row). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 GC and CCMF segmentation results on a image from the GrabCut database.

  The appendix includes some additional technical points refereed in the chapters, and an index is available at the end of the manuscript. Finally, the conclusion contains a summary of the work, some perspectives, and a table presenting the studied and proposed energy formulations on a unified way.

	This generalization allows us to exhibit a new case, for which we developed a globally op-
	timal optimization method. In Chapter 5, we further generalize and extend the framework
	to applications beyond image segmentation, specifically to image filtering, stereovision and
	surface reconstruction.

In the next two chapters, two new flow-based methods are introduced. Chapter 2 deals with image segmentation, and Chapter 3 adresses image restoration problems. The optimization tools introduced in Sections 1.3.2 and 1.3.3 are employed for solving the considered problems, exploiting primal-dual formulations. In Chapter 4, we introduce a framework that unifies and generalizes several state-of-the-art graph-based segmentation algorithms, namely graph cuts, random walker, shortest paths, and watershed.

  tunately, this algorithm has no stopping criteria and has not been proved to converge. Works by Pock et al[START_REF] Pock | A convex formulation of continuous multi-label problems[END_REF], Zach et al[START_REF] Zach | Fast global labeling for real-time stereo using multiple plane sweeps[END_REF] and Chambolle et al[START_REF] Chambolle | A convex approach for computing minimal partitions[END_REF] present different algorithms for

optimizing comparable energies for solving multilabel problems, but again, those algorithms are not proved to converge. Some other works have been presenting provably converging algorithms, but with relatively slow convergence speed. For example, Pock and Chambolle introduce in [54] a general saddle point algorithm that may be used in various applications, but needs half an hour to segment a 350 × 450 image on a CPU and still more than a dozen seconds on a GPU (Details are given in the experiments section).

  1.1. The vectors λ ∈ R n and ν ∈ R n correspond to two Lagrange multipliers.

	The expression of the CCMF dual may be written in summation form as
			weighted cut			smoothness term
	min λ,ν	vi∈V	λ i g 2 i	+	1 4	eij ∈E\{s,t}

  Compute the primal-dual search direction ∆ y such that M ∆ y = r. 3. Determine a step length s > 0 and set y = y + s∆ y . (y = [F, λ, ν] ⊤ ) until ||r p || 2 ≤ ǫ, ||r d || 2 ≤ ǫ, and η ≤ ǫ.

	ǫ > 0.
	Result: F solution to the CCMF problem (2.6) such that F st is maximized under the
	divergence free and capacity constraints, and ν, λ solution to the CCMF dual
	problem (2.8).
	repeat
	1. Compute the surrogate duality gap η = -f (F ) ⊤ λ and set t = µn/ η. 2.

Table 2 .

 2 2: Dice coefficient (percentage) computed between the segmentation mask and the ground truth image (provided by GrabCut database).

	The con-

  .26) This interpretation provides more insight in the choice for C. Another interest of the formulation in (3.21) is that it allows us to propose an alternative parallel algorithm to solve the problem. More precisely, we propose to employ Algorithm 6 which is derived from the M+S Forward-backward-forward algorithm in [46][Prop. 4.4], recalled in Section 1.3.3.2.

Table 3

 3

	σ 2	5	10	15	20	25	50	100
	SNR values obtained by optimizing weighted TV with Split Bregman
	house	32.5 28.8 27.0 25.7	25.0	21.6	17.3
	man	28.4 24.8 23.0 21.7	20.6	17.2	12.8
	lena	31.5 27.7 25.9 25.0	24.2	21.2	16.9
	barbara 27.1 22.8 20.2 18.5	17.3	15.2	12.8
	mean 29.9 26.0 24.0 22.7 21.8 18.8	14.9
			SNR values obtained by optimizing DCTV
	house	32.6 29.2 27.3 25.7	25.1	21.6	16.9
	man	28.9 25.5 23.4 21.9	20.8	17.0	12.6
	lena	31.9 28.5 26.5 25.1	24.3	21.1	16.7
	barbara 27.1 23.0 20.7 19.0	17.7	15.1	12.5
	mean 30.1 26.6 24.5 22.9 22.0 18.7	14.7

.1 reports SNR values for DCTV and weighted TV results obtained on each image corrupted

Table 3 .

 3 1: Quantitative denoising experiment on standard images corrupted with additive Gaussian noise with variance σ 2 . The edges weights were set to The regularization parameter

λ was chosen empirically to give the best results for both methods.

Table 3 .

 3 .8. Furthermore, the results are obtained twice faster, for the same number of iterations, using a Matlab implementations for both methods.

	(σ 2 , b)	(5,5)	(5,7)	(10,5)	(10,7)
	SNR values obtained using Wiener deconvolution
	house	19.2	19.8	18.8	19.2
	lena	18.2	18.9	18.1	18.8
	barbara	14.7	14.5	14.2	14.0
	man	15.0	15.6	15.2	15.9
	mean	16.8	17.2	16.5	17.2
	SNR values obtained using Hybrid TV [66]
	house	26.2	25.3	24.9	23.8
	lena	24.9	23.8	23.9	22.9
	barbara	16.4	15.7	15.8	15.4
	man	21.6	20.6	20.6	19.6
	mean	22.3	21.4	21.3	20.4
		SNR values obtained by optimizing DCTV
	house	26.3	24.7	24.6	23.3
	lena	25.5	23.9	24.0	22.8
	barbara	16.9	15.8	16.0	15.5
	man	22.1	20.6	20.8	19.6
	mean	22.7	21.3	21.4	20.3

2: Quantitative deconvolution experiment on standard images corrupted with additive white Gaussian noise with variance σ 2 , and a Gaussian blur of kernel size b × b. The DCTV weights have been set according to

Table 6 .

 6

		x∈{0,1} n	w ⊤ |Ax| + D(x)	See Sect. 1.2.2
	Combinatorial Con-tinuous Max Flow	λ>0, δ∈R + arg min ν∈{δ,1+δ}	λ ⊤ g 2 + (	1 4	• /(|A|λ)) ⊤ (Aν) 2 + D(ν, λ) See Sect. 2.2.3
	Combinatorial Total Variations	arg min x	g ⊤ (|A| ⊤ (Ax) 2 )	1 2 + D(x)	See Sect. 2.3.3
	Dual Total Variations Constrained	arg min x	max F ∈C	F ⊤ (Ax) + D(x)	See Sect. 3.2
	Combinatorial Dirichlet	arg min x	(Ax) ⊤ diag(w)(Ax) + D(x)	See Sect. 1.2.3
	Shortest Path Forest (geodesics)	lim
	Power watershed	lim				See Sect. 4.3.2
	ℓ 0 -norm anisotropic diffusion	lim α→∞	arg min x	σ α (Ax) + D(x)	See Sect. 5.1.2.1

q→∞ arg min

x w ⊤ q |Ax| q + D(x)

See

[START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF]

p→∞ arg min x w ⊤ p |Ax| q + D(x) 1: The different combinatorial regularization terms studied and proposed in this thesis. The solution x is a labeling defined on a graph with incidence matrix A, edge weights w and node weights g that are assumed to be positive. The data fidelity term is noted D. More details about each approach are given in the designated sections and references.

This chapter was published in part in Proc. of ICASSP

[START_REF] Couprie | Dual constrained TV-based regularization[END_REF]

Published in Proc. of ICIP

[START_REF] Couprie | Anisotropic diffusion using power watersheds[END_REF] 2 Published in Proc. of ISMM

[69]

Saddle point problem, 49 Saturated edge, 10 Seeded image segmentation, 8 Seeds, 8 Shortest Path forests, 10 Shortest Paths, 10 Shortest paths (geodesics), 71, 74 Shrinking bias, 87 Sink, 9 Source, 9 Spanning tree, 11 Spliting method, 50 Stereovision, 100 Strictly convex function, 77 Surface reconstruction, 102 Total variation, 6, 31, 47, 106 Touch expand, 108 Transport graph, 9, 25 Tree, 11 Unary terms, 40, 82 Watershed, 70, 75, 78 Weak duality, 23, 32 Weight function, 26, 51, 73

Acknowledgments / Remerciements

I would like to acknowledge, in no particular order, the three persons who super-

J'ai bénéficié au long de ma thèse d

Algorithm 7: 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.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. phantom seeds v F and v B , i.e.,

We prove in the next section that the labeling x obtained by Algorithm 7 optimizes (4.4).

An illustration for this section is given in Figure 4.3. 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.4) regardless of the value of p.

An implementation of Algorithm 7 when q = 2 can be downloaded from sourceforge [START_REF] Couprie | Power watersheds: a unifying graph based optimization framework[END_REF].

Some implementation details are given in Appendix A.2.6.

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.4) subject to the boundary constraints. We note the labeling obtained by Algorithm 7 by x.

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

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 independent of β. We used the color gradient given by max 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, as it is with this set of seeds the worst performer. Graph cuts

Formulation

Black et al. showed that anisotropic diffusion could be viewed as the minimization of a robust estimator. We employ this same robust estimator formulation of anisotropic diffusion, but

show that the power watershed may be used to iteratively optimize the formulation.

Given the classical continuous form

the anisotropic diffusion algorithm is defined using the discrete calculus notations as

where x is the image intensities of a filtered image, and x at time 0 equals the input, unfiltered image, and A is, as usual, the incidence matrix of a graph. The function g(x) is used to prevent blurring over edges and may be of any decreasing form. One possibility for g(x) suggested in [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] is

where α is a free parameter.

The typical method for solving the anisotropic diffusion equation in (5.2) is via a forward Euler method in which the iteration

is applied, using a time step dt which is designed to obey the CFL conditions and ensure stability of the solution. In practice, the number of iterations applied to produce the filtered image is a free parameter which is set to define the desired level of smoothing.

Black et al. [START_REF] Black | Robust anisotropic diffusion[END_REF] showed that the anisotropic diffusion equation given by (5.2) could be viewed as the gradient of the energy weights have been chosen equal to a constant.

The result obtained by the power watershed algorithm is given in Figure 5.8.

Figure 5.8: Two images and the stereo reconstruction using a maximum spanning forest algorithm. This preliminary result is obtained using a different graph than the one described in Fig. 5.7, and is composed of different levels, each level corresponding to a different disparity.

This prospective result shows that power watershed is suitable to dense stereo estimation.

Future work will determine if the global energy optimized with power watershed brings benefits to stereo applications. Much work remains to be done, exploiting the multilabel character of the power watershed. In a transport graph G with m edges, n nodes, we define a vector c of R m , composed of zeros except for the element corresponding to the source/sink edge which is -1. Let λ and ν be two vectors of R n . The CCMF problem max

has for dual min

equivalently written in (2.10), and the optimal solution (F * , λ * , ν * ) verifies max

and the n following equalities

Proof. The Lagrangian of (1), whose objective is equivalent to min F c ⊤ F , is

with λ ∈ R n and ν ∈ R n two Lagrange multipliers. We have to find F in the dual function is such as ∇ F L(F, λ, ν) = 0.

Substituting F in the Lagrangian function, we obtain

The Lagrangian may be simplified as:

The dual of (1) is also min

Furthermore, the KKT optimality conditions give

Using the expression of F * from equation ( 5), and substituting it in (7), we obtain

incidence matrix A of the transport graph deprived of the sink-source edge:

As F verifies the divergence free constraint

and as the equation u s -u t = 1 may be written equivalently a ⊤ u = 1, we conclude that

and by weak duality (Property 6),

A.

and C the capacity convex set

Appendix 123

We obtain the following problem min

optimized by the parallel proximal algorithm given in Algo 9.

Algorithm 9: Parallel proximal algorithm optimizing CCMF ( 15)

From the optimal flow vector F , one can obtain a map of the saturated nodes (See Def. 2.2.1), as displayed in Figure 2, but the convergence is very slow, and subject to large oscillations, making it difficult to find a stopping criterium. Furthermore, the value of F is not sufficient to obtain the optimal value of ν and λ that can not be recovered analytically.

With some patience, after about one hour (on CPU) for the image of Fig. 2, a segmentation can be obtained, by propagation of the seeds into the saturation map.

The most important issue is that this algorithm only gives a value of the flow, and the final deduced segmentation is not very robust when the optimal value of flow is not precise enough.

A.1.3.2 First order method for solving the dual problem

In order to obtain directly the optimal labeling, we perform a direct optimization of the CCMF dual problem. Our goal is to optimize It only remains to bound the middle part of ( 25)

with M 2 the number of edges of weight inferior to w M , and w M 2 the second maximum weight.

Thus we have

A.2.5 Interpretation as a Markov Random Field

An optimum x ⋆ of Eq. 4.5 may be interpreted as a probability for each pixel to belong to the object (as opposed to background). More rigorously, as shown in [START_REF] Singaraju | Continuous valued MRFs for image segmentation[END_REF], 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 [START_REF] Angulo | Stochastic watershed segmentation[END_REF], who use random placements of seeds to define the most probable locations of the watershed lines).

In this section, we follow the development of [START_REF] Singaraju | Continuous valued MRFs for image segmentation[END_REF], 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

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

However, the generalized mean value theorem in [START_REF] Singaraju | Continuous valued MRFs for image segmentation[END_REF] guarantees that the optimal solution component tree, some quasi-linear [START_REF] Najman | Building the component tree in quasi-linear 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] Matas | Parallel Algorithm for Concurrent Computation of Connected Component Tree[END_REF][START_REF] Wilkinson | Concurrent computation of attribute filters on shared memory parallel machines[END_REF]. From Prop. 8, it is easy to see how to use this tree in Algorithm 7. 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 has been used as a preprocessing step for watershed segmentation with markers. It is called geodesic reconstruction from the markers [START_REF] Beucher | The morphological approach to segmentation: The watershed transformation[END_REF][START_REF] Meyer | Morphological segmentation[END_REF], 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 9. 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. 9 also suggests that geodesic reconstruction can be used as a preprocessing in Alg. 7.

Note that there exist some very efficient and easy to implement algorithms to compute a geodesic reconstruction [START_REF] Géraud | Algorithms for mathematical morphology[END_REF][START_REF] Najman | Geodesic Saliency of Watershed Contours and Hierarchical Segmentation[END_REF][START_REF] Vincent | Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms[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. 9 also suggests links between our framework and the classical watershed-based segmentation framework [START_REF] Beucher | The morphological approach to segmentation: The watershed transformation[END_REF][START_REF] Meyer | Morphological segmentation[END_REF][START_REF] Meyer | Segmentation, minimum spanning tree and hierarchies[END_REF]. The framework of watershed cuts [START_REF] Cousty | Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle[END_REF][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 10. 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.