Xavier Desquesnes
email: xavier.desquesnes@unicaen.fr

Abderrahim Elmoataz

Olivier Lézoray
email: olivier.lezoray@unicaen.fr

process for local and non-local image and data processing

Keywords: Eikonal equation, weighted graph, non-local image processing, active contour, PdE, fast marching, high dimensional data

HAL is a

Introduction

The main goal of this paper is to propose an adaptation of the Eikonal equation on weighted graphs, using the framework of Partial difference Equations [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF], with the motivation of extending this equation's applications, to any discrete data that can be represented by graphs. This adaptation leads to a finite difference equation whose coefficients are data geometry dependent, and that leads to an efficient algorithm, as an extension of the Fast Marching algorithm, to propagate multiple fronts without restriction in the direction of their Many applications involve data defined on topologically complex domains. These data can be defined on manifolds (e.g., a sphere) or irregularly shaped domains, defined on network-like structures (e.g., network communities), or defined as high dimensional point clouds such as collections of features vectors. Such unorganized data can be conveniently represented as graphs, where the vertices represent initial data and the edges represent interactions between them. Moreover, the use of a graph representation for usual images also enables to take into account local and non-local interactions and leads to very powerful tools for non-local image processing [START_REF] Buades | Nonlocal image and movie denoising[END_REF][START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF].

Processing and analyzing such structured types of data is a major challenge for both the image and machine learning communities. Hence, it is very important to transfer many tools which were initially developed on usual Euclidean spaces and proven to be efficient for many problems dealing with usual image and signal domains, to graphs and networks.

Classical approaches for graph processing mainly come from graph theory and one can quote two main categories for these methods. Methods of the first category are based on the minimization of an energy with applications in semisupervised segmentation. One can quote graph cuts [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in n-d images[END_REF], random walks [START_REF] Grady | Random walks for image segmentation[END_REF] or recently the power watershed [START_REF] Couprie | Power Watersheds: A Unifying Graph-Based Optimization Framework IEEE Trans[END_REF]. A second category groups techniques based on spectral graph theory [START_REF] Chung | Spectral graph theory[END_REF]. They have been successfully used for image filtering [START_REF] Hein | Manifold denoising[END_REF], image segmentation [START_REF] Shi | Normalized cuts and image segmentation[END_REF], data clustering [START_REF] Von Luxburg | A tutorial on spectral clustering[END_REF], or network communities extraction in complex networks [START_REF] Kurucz | Spectral clustering in social networks[END_REF], and so on.

There has been also recently much interest in transposing signal processing tools used in image and signal process-ing on graphs. One can quote the generalization of wavelets approach to graphs, with the work of Coifman et al. on diffusion wavelets [START_REF] Coifman | Diffusion wavelets[END_REF], Jansen on multiscale methods [START_REF] Jansen | Second Generation Wavelets and Applications[END_REF], or recently Hammond et al. on spectral transform [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF].

Similarly, there are recent works that aim to transpose Partial Differential Equations (PDEs) on graphs. These works exploit discrete calculus to perform such transcription [START_REF] Hirani | Discrete exterior calculus[END_REF]. Discrete Calculus has been used in recent years to produce a reformulation of continuous problems onto a graph is such a manner that the solution behaves analogously to the continuous formulation. See [START_REF] Grady | Discrete Calculus Applied Anal[END_REF] and references therein for a complete overview on that subject with applications in image processing and machine learning.

To transpose PDEs on graphs, one approach consists in exploiting Partial difference Equations (PdEs) on graphs. Conceptually, PdEs are methods that mimic PDEs on the graphs general domain, by replacing differential operators by difference operator on graphs. Historically, it was first introduced in the seminal paper of Courant, Friedrichs and Lewy [START_REF] Courant | On the partial difference equations of mathematical physics[END_REF]. Then, the study of PdEs has appeared to be a subject on its own interest, dealing with existence and qualitative behavior problems [START_REF] Bensoussan | Difference equations on weighted graphs[END_REF][START_REF] Neuberger | Nonlinear elliptic partial difference equations on graphs[END_REF][START_REF] Park | Positive solutions for discrete boundary value problems involving the p-laplacian with potential terms[END_REF]. Introduction of such methods for image processing started with the work of Chan et al. [START_REF] Chan | The digital tv filter and nonlinear denoising[END_REF] who introduced the TV digital filter for image denoising, which is the discrete analogue of total variation on unweighted graphs. Zhou has also used TV on graphs for semi-supervised classification [START_REF] Zhou | Discrete Regularization[END_REF].

Following the line of research we developed in previous works, we base the contributions of this paper on difference operators on graphs [START_REF] Hein | Manifold denoising[END_REF]. The motivation is that these operators allows to simply adapt continuous formulations to graphs by replacing continuous operators by their discrete adaptation. In particular, it allows most techniques based on the p-Laplacian and gradients to be handled with such operators on graphs in a very straightforward, simple and similar manner [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF][START_REF] Ta | Nonlocal pdes-based morphology on weighted graphs for image and data processing[END_REF]. Such an approach enables an adaptation to graphs that is not necessarily consistent with the continuous formulation (see in [START_REF] Grady | Discrete Calculus Applied Anal[END_REF]). This point is however not a problem for the paper and will be investigated in future works.

In previous works, using the PdEs formalism, we have introduced non-local difference operators on graphs, and used the framework of PdEs to transcribe PDEs on graphs [START_REF] Bougleux | Discrete regularization on weighted graphs for image and mesh filtering[END_REF]. In particular, in [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF], we have introduced a non-local discrete regularization on graphs of arbitrary topologies as a framework for image and data filtering and clustering. With the same ideas, we have proposed PdE morphological processes on graphs that are a transcription of continuous morphological PDEs [START_REF] Ta | Partial difference equations over graphs: Morphological processing of arbitrary discrete data[END_REF]. Recently, we have also adapted a time -dependent version of the Eikonal equation with PdEs morphological processes on graphs [START_REF] Ta | Adaptation of eikonal equation over weighted graphs[END_REF][START_REF] Ta | Nonlocal pdes-based morphology on weighted graphs for image and data processing[END_REF].

Eikonal equation background

The Eikonal equation is a special case of the following general continuous static Hamilton-Jacobi equations:

H(x, f , ∇ f ) = 0 x∈Ω ⊂ IR m f (x) = ψ(x) x∈Γ ⊂ Ω , (1)
where ψ is a positive function defined on a domain Ω and f (x) is the traveling time or distance from source Γ . Then, the Eikonal equation can be expressed by using the following Hamiltonian:

H(x, f , ∇ f ) = ∇ f (x) -P(x), ( 2
)
where P is a given potential function. This equation can be linked to the level-set formulation for advancing fronts introduced by Sethian [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics[END_REF] ∂ φ (x,t) ∂t = F (x) ∇φ (x,t) ,

where φ is the level-set representation of Γ , and F = 1/P. The relation between such formulation and the Eikonal equation stems from the following change of variable : φ (x,t) = tf (x), under the condition that F is positive on the whole domain Ω .

Solutions of static equation ( 2) are usually based on a discretization of the Hamiltonian where the approximations are performed by the Godunov methods [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF] or with the Lax-Friedrich schemes [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics[END_REF]. Numerous numerical schemes have been proposed and investigated for solving the nonlinear system described by [START_REF] Bai | Geodesic matting: A framework for fast interactive image and video segmentation and matting[END_REF]. Among the existing ones, we can quote the following schemes. (i) An iterative scheme has been proposed by [START_REF] Rouy | A viscosity solutions approach to shapefrom-shading[END_REF] based on fixed point methods that solve a quadratic equation. (ii) The fast sweeping methods [START_REF] Zhao | Fast sweeping method for eikonal equations[END_REF] that use Gauss-Seidel type of iterations to update the distance function field. The key point of fast sweeping is to update the points in a certain order. (iii) Tsitsiklis [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF] was the first to develop a Dijkstra like method and proposed an optimal algorithm for solving the Eikonal equation. Based on the same idea, Sethian [START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF][START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics[END_REF] proposed the fast marching methods.

Another approach to solve equation ( 2) is to consider a time dependent version of the equation:

     ∂ f (x,t) ∂t = -∇ f (x) + P(x) x∈Ω ⊂ IR m f (x,t) = ψ(x) x∈Γ ⊂ IR m f (x, 0) = ψ 0 (x) x∈Ω . ( 4
)
At steady state, the solution of the system (4) satisfies the Eikonal equation [START_REF] Bai | Geodesic matting: A framework for fast interactive image and video segmentation and matting[END_REF]. Recently, we proposed in a conference paper [START_REF] Ta | Adaptation of eikonal equation over weighted graphs[END_REF] an adaptation of the time dependent formulation of the Eikonal equation over weighted graphs. Based on PdE, the analogue of (4) on a weighted graph

G = (V, E, w) is      ∂ f (u,t) ∂t = -(∇ - w f )(u) p + P(u) u∈V f (u,t)=ψ(u) u∈V 0 ⊂ V f (u, 0)=ψ 0 (u) u∈V (5)
where V corresponds to the set of vertices of the graph and V 0 is the initial set of seed vertices. Operator ∇ - w corresponds to the weighted internal morphological gradient on graphs (detailed in Sect. 2.2) and . p denotes the L p -norm. One can see that formulation [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in n-d images[END_REF] needs numerous iterations due to finite propagation speed and CFL conditions to converge to the solution of the Eikonal equation.

Contributions

This work generalizes and extends significantly our previous works on the Eikonal equation. First, we propose an adaptation of the stationary version of the Eikonal equation over arbitrary weighted graphs. Based on PdEs, our adaptation of ( 2) is given by this finite difference equation

F (u) (∇ - w f )(u) p = 1 u ∈ V f (u) = 0 u ∈ V 0 . (6)
Explicits solutions of this equation are given for particular values of p ∈ {1, 2, ∞}. An efficient algorithm to obtain such solutions, using the Fast Marching's updating scheme, is proposed, and proofs of existence and uniqueness are also provided. This formulation generalizes front propagation methods on weighted graphs and recovers well known schemes as Osher-Sethian discretization or Dijkstra shortest path, for specific graphs and values of p. Such an adaptation of the Eikonal equation on graphs enables the transcriptions of efficient algorithms from the field of image processing to a huge variety of discrete data that can be represented by a weighted graph. Then, using this adaptation, we also propose a new fast algorithm for propagation and tracking of many concurrent fronts on a weighted graph, the complexity of which is independent of the number of these fronts. Such an algorithm leads to several applications on weighted graphs such as semi-supervised image segmentation or data clustering.

In these two previous contributions, we only considered the case where the front is evolving in the outward normal direction (with the speed F defined non-negative everywhere), but it is also interesting to consider both inward and outward directions and both positive and negative speeds. In particular, in the case of semi-supervised clustering, such inward and outward evolutions enables to minimize or overcome errors due to wrong initialization. Finally, we generalize the previous algorithm to the case where the speeds of the different fronts can be either positive or negative, which leads to a graph-based active contour model with many contours. This generalization provides a complete tool for multiple fronts propagation on arbitrary graphs, and offers a novel extension to the classical Fast Marching, when applied to regular grid graphs.

Paper organization

The rest of this paper is organized as follows. In Section 2, we provide definitions and notations used in this work. In Section 3, we present our new finite difference equation on weighted graphs with proofs for existence and uniqueness, and explicit solutions for values of p ∈ {1, 2, ∞} are also given. Section 4 introduces two efficient algorithms for labels propagation, using the previous equation and an adaptation of the Fast Marching algorithm. Then, Section 5 presents several experiments which illustrate the behavior and efficiency of the proposed formulations and algorithms, as geodesic distance, semi-supervised image segmentation (with non-local configuration), active contours or data clustering. Finally, Section 6 concludes this paper.

Operators on graphs

As the core structure of our approach, in this section we provide notations and basics on weighted graphs. We recall our formulations of differences, morphological differences and gradients on weighted graphs [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] Ta | Partial difference equations over graphs: Morphological processing of arbitrary discrete data[END_REF][START_REF] Ta | Adaptation of eikonal equation over weighted graphs[END_REF]. The latter formulations constitute the basis of our proposed numerical scheme for solving the Eikonal equation on weighted graphs.

Notations

We consider the general situation where any discrete domain can be viewed as a weighted graph. Let G=(V, E, w) be a weighted graph composed of two finite sets: vertices V and weighted edges E ⊆ V ×V . An edge (u, v) ∈ E connects two adjacent (neighbor) vertices u and v . The neighborhood of a vertex u is noted N(u) = {v ∈ V \ {u} : (u, v) ∈ E}. The weight w(u, v) of an edge (u, v) can be defined with a function w : V ×V → IR + if (u, v) ∈ E, and w(u, v) = 0 otherwise. For the sake of simplicity, w(u, v) will be denoted by w uv . Graphs are assumed to be simple, connected and undirected implying that function w is symmetric.

Let f : V → IR be a real-valued function that assigns a real value f (u) to each vertex u ∈ V . We denote by H (V ) the Hilbert space of such functions and similarly by H (E), the Hilbert space of functions that assign a real value to each edge of E. These two spaces are endowed with the following inner products:

f , h H (V ) = ∑ u∈V f (u)h(u), (7)
with f , h ∈ H (V ), and

F, H H (E) = ∑ u∈V ∑ v∈V F(u, v)H(u, v), ( 8
)
where F, G ∈ H (E).

Given a function f ∈ H (V ), the integral of f is defined as

V f = ∑ u∈V f (u) (9)
and it's L p norm is given by

f p = ∑ u∈V | f (u)| p 1/p , 1 p < ∞ f ∞ = max u∈V | f (u)| , p = ∞ (10)
Let A be a set of connected vertices with A ⊂ V such that for all u ∈ A , there exists a vertex v ∈ A with (u, v) ∈ E. We denote by ∂ + A and ∂ -A : the external and internal boundary sets of A , respectively

∂ + A = {u ∈ A c : ∃v ∈ A with (u, v) ∈ E} and ∂ -A = {u ∈ A : ∃v ∈ A c with (u, v) ∈ E} ( 11
)
where

A c = V \ A is the complement of A .

Operators and gradients on weighted graphs

The weighted gradient operator or weighted difference oper-

ator of a function f ∈ H (V ) , noted → ∇ w : H (V ) → H (E), respectively d w , is defined on an edge (u, v) ∈ E by → ∇ w f (u, v) de f . = d w f (u, v) de f . = w(u, v) f (v) -f (u) . ( 12
)
The adjoint of the weighted gradient operator, noted

→ ∇ * w f : H (E) → H (V ), is defined by : → ∇ w f , H de f . = f , → ∇ * w H , (13)
with f ∈ H (V ) and H ∈ H (E), and can be expressed as

→ ∇ * w H (u) de f . = ∑ v∼u w(u, v) H(v, u) -H(u, v) . ( 14
)
where v ∼ u means that v is adjacent to u. This adjoint is linear and measures the flow of a function in H (E) at each vertex of the graph. Similarly to the continuous case, the divergence of a function

F ∈ H (E) is defined by div w F = - → ∇ * w F.
These two definitions of the weighted gradient operator and it's adjoint allow to define a family of first and second order operators on graphs, as the p-Laplace operator. But in this paper we only focus on the first order weighted gradient operator.

Based on the weighted gradient operator definition, two weighted directional gradient operators are defined. The weighted directional external and internal gradient operators are defined as

→ ∇ ± w : H (V ) → H (E), with → ∇ + w f (u, v) de f . = w(u, v) f (v) -f (u) + → ∇ - w f (u, v) de f . = w(u, v) f (v) -f (u) -, (15)
with the following notations: (x) + = max(x, 0) and (x) -=min(x, 0).

The weighted gradient of a function f ∈H (V ) at vertex u is defined as the vector of all weighted gradient with respect to the set of edges (u, v) ∈ E

(∇ w f )(u) de f . = → ∇ w f (u, v) v∈V . (16)
In the sequel, weighted gradient will refer to this gradient defined on vertices. Similarly, the weighted morphological internal and external gradients at a vertex u are expressed as

(∇ + w f )(u) = → ∇ + w f (u, v) v∈V and (∇ - w f )(u) = → ∇ - w f (u, v) v∈V . (17)
The corresponding L p -norms of gradients ( 17) and ( 16) for a vertex u are

(∇ + w f )(u) p = ∑ v∼u w p/2 uv |(D f (u)) + | p 1/p , (∇ - w f )(u) p = ∑ v∼u w p/2 uv |(D f (u)) -| p 1/p and (∇ w f )(u) p = ∑ v∼u w p/2 uv |D f (u)| p 1/p ( 18
)
with 0 < p < ∞ and where D f (u) = ( f (v)f (u)). The relation between these directional gradients norms was given in [START_REF] Ta | Adaptation of eikonal equation over weighted graphs[END_REF] as

(∇ w f )(u) p p = (∇ + w f )(u) p p + (∇ - w f )(u) p p ( 19
)
For the L ∞ -norm, we have

(∇ + w f )(u) ∞ = max v∼u √ w uv |(D f (u)) + | , (∇ - w f )(u) ∞ = max v∼u √ w uv |(D f (u)) -| and (∇ w f )(u) ∞ = max v∼u √ w uv |D f (u)| ,
with the following property :

(∇ w f )(u) ∞ = max( (∇ + w f )(u) ∞ , (∇ - w f )(u) ∞ ). (20)
Properties of these gradients can be found in [START_REF] Ta | Partial difference equations over graphs: Morphological processing of arbitrary discrete data[END_REF][START_REF] Ta | Adaptation of eikonal equation over weighted graphs[END_REF]. One can note that the general definitions presented in this section are defined on graphs of arbitrary topology. They can be used to process any discrete regular or irregular data sets that can be represented by a weighted graph. Moreover, local and non-local settings are directly handled in these definitions and both are expressed by the graph topology in terms of neighborhood connectivity [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF].

Morphological evolution equations

Time dependent Hamilton-Jacobi equation formulated in ( 4) is linked with mathematical morphology processes and can be viewed as morphological evolution equations. We have shown in [START_REF] Ta | Partial difference equations on graphs for mathematical morphology operators overs images and manifolds[END_REF][START_REF] Ta | Partial difference equations over graphs: Morphological processing of arbitrary discrete data[END_REF]

∂ t f = + ∇ f p and ∂ t f = -∇ f p ( 21
)
where f is a modified version of f 0 and with the initial condition ∂ t=0 f = f 0 . With different values of p, one obtains different structuring elements: a rhombus for p = ∞, a disc with p = 2 and a square with p = 1.

We have proposed in [START_REF] Ta | Partial difference equations over graphs: Morphological processing of arbitrary discrete data[END_REF] a graph-based versions of (21) by using gradients defined in [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF]. Given a weighted graph G = (V, E,w) and a function f ∈ H (V ), the analogues of [START_REF] Grady | The piecewise smooth mumford-shah functional on an arbitrary graph[END_REF] on G are

∂ t f (u) = + (∇ + w f )(u) p and ∂ t f (u) = -(∇ - w f )(u) p . (22)
Intuitively, given a set of vertices A ⊂ V and using external and internal graph boundaries [START_REF] Chan | Active contours without edges[END_REF], equation of dilation over A can be interpreted as a growth process that adds vertices from ∂ + A to A . By duality, erosion over A can be interpreted as a contraction process that removes vertices from

∂ -A to A .
As mentioned in the introduction, one approach to solve the Eikonal equation ( 1) is to consider a time dependent version (4) of this equation. The time dependent equation can be viewed as an erosion process regarding the minus sign and a null potential function P. Then, using the corresponding internal gradient ∇ - w involved in discrete PdEs based erosion process, one can directly obtain the time-dependent adaptation of the eikonal equation on graphs. Given a graph G=(V, E, w) and a function f ∈H (V ), we have

     ∂ f (u,t) ∂t = -(∇ - w f )(u) p + P(u) u∈V f (u,t)=ψ(u) u∈V 0 ⊂ V f (u, 0)=ψ 0 (u) u∈V, ( 23
)
where V 0 corresponds to the set of initial seed vertices. This equation can be solved by steepest gradient descent method, using the following iterative numerical scheme for all u ∈ V , with f n (u) ≈ f (u, n∆t):

f n+1 (u) = f n (u) -∆t ∇ - w f n (u) p -P(u) (24)
A complete and precise definition with numerical schemes for p ∈ {1, 2, ∞} can be found in [START_REF] Ta | Adaptation of eikonal equation over weighted graphs[END_REF].

Adaptation of the Eikonal equation over weighted graphs

As mentioned in the introduction, the time dependent version of the Eikonal equation on weighted graphs needs numerous iterations to converge to the solution and our approach, detailed in Section 2.3 is not adapted to large graphs. Therefore, we propose a new adaptation of the Eikonal equation and a new algorithm to solve such finite difference equation ( 6), that overcomes these limitations.

We consider the Eikonal equation that describes the evolution of a propagation front Γ ,

F (x) ∇ f (x) = 1 f (x) = 0, x ∈ Γ , (25)
where f (x) is the arrival time of the front at x and F 0. The associated level-set function φ is defined as

∂ φ (x,t) ∂t = F (x) ∇φ (x,t) , (26)
and we know that φ (x,t) = tf (x), then

∂ (t -f (x)) ∂t = F (x) ∇(t -f (x)) = F (x) ∇ f (x) = 1. ( 27
)
Transposed on a weighted graph G = (V, E, w) with functions f ∈ H (V ) and F ∈ H (V ) and using morphological equations defined in ( 22), the level-set formulation ( 26) can be rewritten as

∂ φ (u,t) ∂t = F (u) ∇ w φ (u,t) p , (28)
and can be expressed as a morphological process with the following sum of dilation and erosion

∂ φ (u,t) ∂t = F (u) + ∇ + w φ (u,t) p p + F (u) -∇ - w φ (u,t) p p 1/p . ( 29
)
Only considering the case F 0, and with φ (u,t

) = t - f (u), one has ∂ φ (u,t) ∂t = F (u) ∇ + w (t -f ) (u) p = F (u) ∇ - w f (u) p = 1. (30)
Finally, with P = 1/F we obtain a discrete adaptation of the Eikonal equation on weighted graph, which describes a morphological erosion process, and defined as

∇ - w f (u) p = P(u). ∀u ∈ V f (u) = 0 ∀u ∈ V 0 (31)

Existence and uniqueness of the solution

The proof of existence and uniqueness of equation ( 31) can be established considering the equation

∇ - w f (u) p = P(u) u ∈ A ⊂ V f (u) = 0 u ∈ ∂ A ( 32
)
which can be expressed as

S u u, f , f (v) v∼u = 0 with S u u, f , [ f (v)] v∼u = ∇ - w f (u) p -P(u) = 0 u ∈ A S u u, f , [ f (v)] v∼u = 0 u ∈ ∂ A. ( 33
)
It can be easily shown that this scheme satisfies the following properties for all f .

-

H 1 . ∂ S u [ f ]/∂ f (v) 0, ∀u = v -H 2 . S u [ f + M] = S u [ f ], ∀M ∈ IR∀u -H 3 . S u [λ f ] = λ S u [ f ] + (λ -1)P(u), ∀λ 0, ∀u -H 4 . lim f (u)→+∞ S u (u, f , f (v) v∼u ) = -P(u) -H 5 . lim f (u)→-∞ S u (u, f , f (v) v∼u ) = +∞
Uniqueness of the solution Equation ( 31) has a unique solution.

Proof Let f and g be two distinct solutions of (31) such that max

A ( f (u) -g(u)) > 0. ( 34
)
Then for a given λ 1, we have

M H = max A ( f (u) -λ g(u)) > 0. ( 35
)
We denote u 0 the u for which

M H = f (u 0 ) -λ g(u 0 ). Then we can deduce that, with h(u) = M H + λ g(u) : f (u) h(u) ∀u ∈ A f (u 0 ) = h(u 0 ) (36)
This implies that

0 = S u 0 [ f ] S u 0 [h] (due to H 1 ) = S u 0 [λ g] (due to H 2 ) = λ S u 0 [g] + (λ -1)P(u 0 ) > 0 (due to H 3 ). (37)
There is a contradiction unless f (u) = g(u).

Existence of solution The proof of existence for equation [START_REF] Jansen | Second Generation Wavelets and Applications[END_REF]'s solution can be easily shown. We know that equation

S u (u, f , f (v) v∼u ) is continuous. Due to properties H 4 and H 5 , we can deduce that ∃ -∞ < f (u) < +∞ such that S u (u, f , f (v) v∼u ) = 0.

Numerical schemes and algorithms

From (31) and using norms defined in ( 18) and ( 20) with the property min(x, 0) =max(-x, 0) , we obtain the following equations for the L p and L ∞ norms.

-Case p ∈ {1, 2}:

∑ v∼u w p/2 uv max 0, f (u) -f (v) p 1/p = P(u). ( 38
)
With a simple transformation of variables and some conventional notations, equation ( 38) can be rewritten as

n ∑ i=1 (x -a i ) + h i p = C p , (39)
where

x = f (u), h i = 1/w uv , n = card(N(u)), a = { f (v i )|v i ∈ N(u) with i = 1, ..., n}, C = P(u). -Case p = ∞ : max v∼u √ w uv max 0, f (u) -f (v) = P(u). ( 40
)
Using the same transformation of variables we obtain

max i (x -a i ) + h i = C . (41)
Local solutions (i.e., solution for a vertex, assuming the others are held fixed) of equation ( 31) over a weighted graph are given by ( 39) and [START_REF] Shamir | A survey on Mesh Segmentation Techniques[END_REF]. Both equations are clearly independent of the graph formulation, and can be applied to weighted graphs of arbitrary topology.

For the case p = {1, 2}, local solution x at a particular vertex can be easily obtained with an iterative algorithm as described in Algo. 1. The algorithm is based on the knowledge that there exists a k with 1 k n such that a k x a k+1 , and x is the unique solution of the equation. Then, the algorithm consists in sorting increasingly the values a i and computing temporary solutions xm with the following equations [START_REF] Shi | Normalized cuts and image segmentation[END_REF] and [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF] until the condition xm a m+1 is satisfied.

For p = 1, the temporary local solution is given by :

xm = m ∑ i=1 a i h i + C m ∑ i=1 1 h i . ( 42
) Algorithm 1 x computation (Local solution) Sort increasingly the a i , i = 1, ..., n a n+1 ← ∞ m ← 1 x ← ∞ while xm a m+1 and m n -1 do xm ← solution of ∑ m i=1 x-a i h i p = C p with p = 1 or 2 m ← m + 1 end while x ← xm For p = 2, one has : xm = m ∑ i=1 a i h 2 i + m ∑ i=1 C 2 h 2 i -∑ j>i (a i -a j ) 2 h 2 i h 2 j m ∑ i=1 1 h 2 i ( 43
)
The numerical scheme for the L ∞ norm is much simpler.

Considering equation ( 41), the local solution x for a particular node is directly given by

x = n min i=1 (a i + h i C ), (44)
where n corresponds to the number of neighbors. One can remark that this equation is equivalent to Dijkstra like local update equation.

Many solvers can be adapted to compute global solution (i.e., on the whole graph) as, for instance, the Fast Iterative Method introduced by Jeong and Whitaker [START_REF] Won-Ki | A fast iterative method for eikonal equations[END_REF]. In this paper, we prefer the Fast Marching method which has the advantage to be monotonic. On an arbitrary graph, the Fast Marching consists in an active list (A) of vertices for which the solution is already known and fixed, and in a narrow band (NB) of vertices which are not yet fixed and have at least one neighbor in the active list. Vertices which are neither active or in the narrow band are said far away (FA). The narrow band is built as a sorted heap, and at each iteration, the first vertex is removed from the narrow band, added to the active list, and it's neighbors are updated if they are not yet fixed and added to the narrow band if they are far away. Each neighbors v are updated simultaneously by computing new local solutions x v from the new value of f (u) (and using Algo.1) and updating f (v) ← min f (v), x v under the condition that the new local solution x v is inferior to previous local solution. This is iterated until the narrow band is empty. More details on the initial Fast Marching algorithm and a proof that the initial algorithm constructs a viable solution are given in [START_REF] Sethian | A fast marching level set methods for monotonically advancing fronts[END_REF].

Accuracy and Complexity of the Algorithm Let G = (V, E, w) be a weighted graph. The costs C u,p to update a vertex u (i.e., compute a new value of f (u)), according to Algo. 1 and equations ( 41), ( 42) and [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF], are given by

C u,1 (u) = D u (u) 2 case p = 1, C u,2 (u) = D u (u) 3 case p = 2, C u,∞ (u) = D u (u) case p = ∞, (45)
where the quantity D u (u) corresponds to the number of neighbors of u which are not far away.

In order to generalize the cost functions C u,p , we can consider that the previous quantity D u (u) is constant and corresponds to the maximum degree of a vertex of V . Such constant is defined as D = max u∈V (card(N(u))). Thus, in the sequel we will use the following constants to denote the update costs.

C u,1 = D 2 C u,1 (u)∀u ∈ V, C u,2 = D 3 C u,2 (u)∀u ∈ V, C u,∞ = D C u,∞ (u)∀u ∈ V. ( 46
)
Based on the previous constants, the activation costs C a,p of a vertex (i.e., the cost to change the state of a vertex from far away to active, what implies to update every of it's not yet activated neighbors) can be defined as

C a,p = DC u,p , p ∈ {1, 2, ∞}. (47)
while we know that N(u) D ∀u ∈ V . And we have

C a,1 = D 3 case p = 1, C a,2 = D 4 case p = 2, C a,∞ = D 2 case p = ∞. (48)
Finally, the total cost of the algorithm can be summarized as the following complexity

O(C a,p N log(N)) ( 49
)
where N is the number of vertices and the log(N) factor corresponds to the managing cost of the heapsort. In practice, processing such an algorithm on strongly connected graphs is very rare, and in the most cases we have C a,p ≪ N.

Grid graph example

In this section, we show that with an adapted graph topology and an appropriated weight function, the proposed formulation recovers the well known Osher-Sethian scheme that has been proposed in literature to solve the Eikonal equation.

Let G = (V, E, w) be a weighted graph. In the case where p = 2, we have the following scheme

∑ v∼u w uv max 0, f (u) -f (v) 2 = P(u). ( 50
)
If the graph represents a m-dimensional grid in IR m , this numerical scheme recovers the Osher-Sethian discretization models on a m-dimensional grid.

Let u be a given vertex of V , that defines a vector of mdimensional spatial such that u = (i 1 h 1 , ..., i m h m ) where the h j represent the grid spacing and i j ∈ IN. The neighborhood of u can then be defined as N(u) = {v|v = u ± h j e j with j ∈ 1, ..., m} where e j = (q k ) k=1,...,m is a vector of IR m such as q k = 1 if j = k and q k = 0 otherwise. Finally, with the notations

D + j f (u) = ( f (u + h j e j ) -f (u)) /h j , D - j f (u) = ( f (u) -f (u -h j e j )) /h j ( 51
)
and with w uv j = 1/h j for all v j ∈ N(u), [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF] can be rewritten as 2 . This equation ( 52) corresponds to the Osher-Sethian Hamiltonian discretization scheme on a m-dimensional grid, and can be solved by the initial Fast Marching algorithm [START_REF] Sethian | A fast marching level set methods for monotonically advancing fronts[END_REF].

m ∑ j=1 max 0, D - j f (u) 2 + min 0, D + j f (u) 2 = P(u) (52) since min(0, a -b) 2 = max(0, b -a)

A new class of fast algorithms for semi-supervised graph clustering

In Section 3, we have introduced a new finite difference equation [START_REF] Jansen | Second Generation Wavelets and Applications[END_REF] which is an adaptation of the Eikonal equation over weighted graph, and proposed an efficient algorithm based on the Fast Marching method to solve such an equation. In this section, we propose to extend the previous algorithm so as to define a new class of fast algorithms for diffusion processes, that enable the simultaneous propagation of many concurrent labels on a weighted graph. Efficient algorithms for concurrent label propagation on a weighted graph have a great interest for graph partitioning and have led to numerous applications in semi-supervised image segmentation and data clustering. We will first present the algorithm when the speed is always non-negative ( F : V → [0, +∞]), then its extension to the case where the speed can be either positive or negative (F : V → [-∞, +∞]).

Label propagation with non-negative speed

Given a graph G = (V, E, w), let L = {l 1 , .., l n } be a set of labels and S 0 = S 0 1 ∪ ... ∪ S 0 n be the set of vertices initially marked by a label, where S 0 i is the set of vertices initially marked by l i . We call S 0 i the seed of l i . The objective of label propagation for graph clustering is to mark each vertex u of V with a label l i , under the condition that u is closer to at least one vertex of S 0 i than other vertices of S 0 (according to the topology of G, the weight function w and a speed function F l i : V → IR, l i ∈ L), so as to obtain a partition of the graph in n consistent clusters. The propagation of each label l i is driven by a front Γ i , initialized as the boundary of S 0 i , and the final partition S i is given by the set of vertices reached by Γ i until it is stopped by another front or the boundary of the domain.

lab(u) = Initial label of u. f (u) = 0 ∀ u ∈ S 0 ; f (u) = +∞ ∀ u ∈ V \ S 0 s(u) = +∞ ∀ u ∈ V \ S 0 A = S 0 ; NB = {u | ∃ v ∈ A and v ∈ N(u)} FA = V \ A ∪ NB 2. Process: while FA = / 0 do u ← first element of NB remove u from NB and add u in A. for all v ∈ N(u) ∩ A do compute local solution t ← f (v) if t < f (v) then f (v) = t if v ∈ FA then remove v from FA and add v in NB. else update position of v in NB end if if f (u)/w uv < s(v) then s(v) = f (u)/w uv lab(v) = lab(u) end if end if end for end while
In practice, the propagation is performed in the same time than the computation of the arrival-time (given by resolution of (31)): each time a vertex u is reached by a front, the label of the front is propagated to u. The proximity condition is respected while in the case of several front propagation, the front that arrives at a vertex is necessarily the front coming from the nearest source of that vertex (according to weight and speed functions). This is a consequence of the Fast Marching algorithm which activates vertices from the smallest to the greatest distance (or arrival-time). Indeed, a vertex u is considered reached by a front when it is activated by the algorithm, and the propagated label comes necessarily from it's neighbor v which is already activated and such that

s(u) = f (v) w uv = min z∼u f (z) w uz . ( 53
)
The weight 1/w is used to penalize neighbors which are close to a source but have very weak link with u. The whole algorithm for label propagation on weighted graphs is given in Algo. 2.

Remark 1 In the case where several neighbors have the same contribution and where these neighbors have different labels, several labels could be affected to node u. But, as this situation is very rare and happens only where fronts collapse, we choose, in this work, to label the vertex with only one arbitrary label from the list. This choice is application dependent and does not act as a rule.

Remark 2 Due to the multitude of fronts, we have to define one speed per front and the speed of a front Γ i at a vertex u is given by F l i (u). Consequently, each time a node u is activated, it's neighbors v such that v / ∈ A are updated by :

(∇ - w f )(v) p = P l i (v), (54)
where l i is the label of u and

P l i (v) = 1/F l i (v).
One can remark that the algorithm complexity does not depend on the number of labels (the complexity is the same than for the classical Fast Marching algorithm described in Section 3). Then, this algorithm opens the way to a new class of fast algorithms for semi-supervised graph clustering, which allows to cluster data or segment images in many partitions, without loss of efficiency.

As noticed in Section 3, the label propagation is performed with a speed F which is always non-negative. In the following section, we propose an adaptive version of the algorithm that can consider a front evolution (and then label propagation) with positive or negative speed.

Label propagation without restrictions on speed sign

The case where the front's speeds can be either positive or negative is very interesting because it enables the labels to be propagated when the speed is positive and removed when the speed is negative, This also enables the tracking of multiple fronts subjected to changing sign speeds.

In the continuous case, Carlini has proposed an approach for tracking an hypersurface with non always positive velocity, called Generalized Fast Marching Method (GFMM) [START_REF] Carlini | Convergence of a generalized fast-marching method for an eikonal equation with a velocity-changing sign[END_REF].

Let Γ be a front evolving on a continuous domain noted Ω , an represented by a characteristic discontinuous function θ (x,t) defined as

θ t (x,t) = F (x,t)|∇(θ (x,t)| θ (x, 0) = 1 Ω 0 -1 Ω c 0 (55)
where F (x,t) is a changing-sign speed and 1 Ω 0 -1 Ω c 0 is equal to 1 on Ω 0 and -1 on its complementary. The position of the front Γ is then given by the discontinuities of Fig. 1 Label propagation with changing-sign speed. The first graph is a 4-adjacency grid graph, and the second a non-Cartesian graph. Vertices inside the purple line represent S 0 i . The speed is positive on the left side and negative on the right side. Green points and red points represent vertices in initial configuration of positive front Γ + i , respectively vertices in initial configuration of negative front Γ - i . Blue points represent vertices in the narrow band from the front Γ + i (surrounded in green) and from the front Γ - i (surrounded in red). (Color online).

function θ and the time evolution is obtained by solving the stationary Eikonal equation [START_REF] Hirani | Discrete exterior calculus[END_REF]. Then, the main idea of the Generalized Fast Marching Method is to split the front in two fronts: Γ t + = ∂ {x|θ (x,t) > 0} and Γ t -= ∂ {x|θ (x,t) < 0} and perform two Fast Marching along two narrow bands (NB t + = Γ t + ∩ {x, F (x,t) < 0} and NB t -= Γ t -∩ {x, F (x,t) > 0}), with a modified version of (55). Interested readers can refer to [START_REF] Carlini | Convergence of a generalized fast-marching method for an eikonal equation with a velocity-changing sign[END_REF] for a detailed presentation.

Then, we propose an extension of such an approach for front tracking on weighted graphs, along with a generalization of the previous algorithm (Algo. 2) with positive and negative speeds. This can be done by representing each front Γ i by two fronts Γ + i and Γ - i that describe the part of front Γ i subjected to a positive speed and respectively negative speed. Then, the propagation of a label l i is performed through two evolving fronts: Γ + i which describes the growth of the set S i (where the speed F l i is positive) and Γ - i which describes the decay of the set S i (where the speed F l i is negative). We have

Γ i = Γ + i ∪ Γ - i .
On the first hand, the front Γ + i is initialized by the set {u|u ∈ ∂ -S 0 i and F l i (u) > 0}. In other words, this corresponds to the set of vertices of S 0 i which lies in the inner narrow band of S 0 i and where the speed F l i of the front is positive. Neighbors of the front are added to the narrow band under an additional condition: as the speed of the front has to be positive at their position. This additional condition can be easily justified since the front can't grow where its speed is negative.

On the other hand, the front Γ - i is initialized by the set {u|u ∈ ∂ + S 0 i and F l i (u) < 0}. Neighbors are added to the narrow band under the condition than the speed F l i is negative. Because this front describes the decay of S i , each time a vertex is reached by Γ - i , the vertex is unmarked instead of being marked by l i .

Finally, in order to keep the distance (or arrival time) always positive, the signed version of the speed is replaced by it's modulus and the potential is given by

P l i (u) = 1/|F l i (u)|.
Such an approach of label propagation with changingsign speed for two fronts is illustrated in Fig. 1 and Fig. 2. The entire process is summarized in Algo. Fig. 2 Illustration of label propagation with changing-sign speed, with a toy example. The set S i is represented in blue. The inner narrow band which represents the position of front Γ + i is superimposed in green, and the outer narrow band which represents the front Γ - i is superimposed in red. The six images show different steps of the evolution of the set S i and both two fronts, with the speed F l i which is positive on the gray background and negative on the white background. (Color online).

Remark 3 In the case where p = ∞, and according to equation [START_REF] Ta | Partial difference equations on graphs for mathematical morphology operators overs images and manifolds[END_REF], the arrival-time of a vertex is computed with the contribution of exactly one edge. Geometrically, this means

A = {u | u ∈ ∂ + S 0 i and F l i (u) < 0} ∪ {u | u ∈ ∂ -S 0 i and F l i (u) > 0} NB = {u | ∃ v ∈ A ∩ N(u) and F l i (v) × F l i (u) > 0} FA = V \ (A ∪ NB) f (u) = 0; lab(u) = l i ∀ u ∈ A f (u) = min(1/w uv ), v ∈ A ∩ N(u) lab(u) = l i ; s(u) = f (u) ∀ u ∈ NB f (u) = +∞; lab(u) = 0; s(u) = +∞ ∀ u ∈ FA 2. Process: while FA = / 0 do u ← first element of NB remove u from NB and add u in A. for all v ∈ N(u) ∩ A and F l i (v) × F l i (u) > 0 do compute local solution t ← f (v) if t < f (v) then f (v) = t if v ∈ FA then remove v from FA and add v in NB. else update position of v in NB end if if f (u)/w uv < s(v) then s(v) = f (u)/w uv lab(v) = lab(u) end if end if end for end while
that the fronts propagation follows the edges (one can see edges as pipes in which the fronts evolve), and the front that reaches a particular vertex comes from exactly one of its neighbors. The label assignment is then automatically deduced.

On the contrary, in the case where p = {1, 2} and according to equations ( 42), [START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF], the arrival-time of a vertex is computed with the contribution of several edges. Thus, only considering the graph structure, the front propagation is interpreted as coming from several edges. In practice, we consider that the front comes from the edge with the highest contribution, and labels are propagated in this way. If we consider a sufficiently higher-dimensional domain in which the graph is embedded (as it is the case for images, meshes or manifolds), the front is coming between all the involved edges.

Experiments and Applications

Graph construction

There exists several popular methods that transform a set of unorganized discrete data into a neighborhood graph. Considering a set of vertices V such that functions of H (V ) represent the data, the construction of such a graph consists in modeling the neighborhood relationships between the data through the definition of a set of edges E, using a pairwise distance measure µ : V × V → IR + . In the general case (where data are unorganized) one can quote:

ε-neighborhood graphs where two data v i , v j ∈ V are connected by an edge of E if µ(v i , v j ) ε, with ε > 0.

k-nearest neighbor graphs (k-NNG) where each vertex v i is connected with it's k-nearest neighbors according to µ. Such construction implies to build a directed graph, as the neighborhood relationship is not symmetric. Nevertheless, an undirected graph can be obtained while adding an edge between two nodes v i and v j if v i is among the k-nearest neighbors of v j or if v j is among the knearest neighbors of v i .

In the case of structured data (i.e., images or meshes):

-Grid graphs which are the natural structure to describe an image with a graph. Each pixel is connected by an edge to its adjacent pixels. -Region adjacency graphs (RAG), initially designed for images where vertices correspond to image regions and the set of edges is obtained by considering an adjacency distance. A RAG can be built for any structured data represented by a graph, where a region R i is defined as a set of connected vertices such that R i = V and

R i = / 0. Two regions R i , R j are adjacent if ∃ v i ∈ R i and v j ∈ R j | v i ∼ v j .
In both cases, weights are computed according to a measure of similarity g : E → IR + , which satisfies:

w(u, v) = g(u, v) if (u, v) ∈ E 0 otherwise ( 56
)
The similarity is usually based on a pairwise distance measure between data features, where each vertex u ∈ V is represented by a feature vector F u ⊂ IR m . For a given edge (u, v) ∈ E and a distance measure µ : V ×V → IR + , we can have

g 0 (u, v) = 1, g 1 (u, v) = (µ(F u , F v ) + ε) -1 with ε > 0, ε → 0, g 2 (u, v) = exp(-µ(F u , F v )/σ 2 ) with σ > 0, ( 57
)
where σ depends on the variation of the function µ over the graph and controls the similarity scale. Several choices can be considered for the expression of the feature vectors, depending on the nature of the features to be used for the graph processing. In the context of image processing, one can quote the simplest grayscale or color feature vector F u , or the patch feature vector F τ u = v∈W τ (u) F v (i.e, the set of values F v where v is in a square window W τ (u) of size (2τ + 1) × (2τ + 1) centered at a vertex pixel u), in order to incorporate non-local features. In this section, we show the adaptivity of our framework and illustrate the behavior of proposed numerical schemes to compute weighted geodesic distances. Figure 3 presents the application of the schemes formulated by ( 39) and [START_REF] Shamir | A survey on Mesh Segmentation Techniques[END_REF]. Results are obtained with different graph topologies and weight functions. They are all illustrated with color distance maps where iso-levels are superimposed in white. The potential function P is constant and equal to 1 and the propagation is performed for a unique label which initial seed is located either at the top left corner of the original image or at the center. The original image is a 256 × 256 grayscale image. First rows show results obtained with a weighted 4-adjacency grid graph and different weight functions w 1 , w 2 and w 3 . The first weight function (w 1 = 1) is a constant weight function. Using such weight function and in the case where p = 2 we obtain the same result as using the initial Fast Marching algorithm [START_REF] Sethian | A fast marching level set methods for monotonically advancing fronts[END_REF], what verifies that in such configuration our formulation recovers the Osher-Sethian scheme. Similarly, in the case where p = ∞, our algorithm is equivalent to a Dijkstra algorithm. The second weight function is designed to favor horizontal diffusion to vertical diffusion by penalizing the weights of vertical edges and given by

Weighted geodesic distances

w 2 (u, v) = 1 if v x = u x 0.25 if v y = u y . ( 58
)
In that latter case, the seed is placed at the center of the image. The third weight function (w 3 = exp(-d 2 /σ 2 )) is designed to hold the similarity between each connected pixel (based on pixel intensity). Shape information is naturally captured by the weights that stop the front evolution at the boundaries. Clearly, in this case, one can obtain a segmentation of the different shapes simply by thresholding the obtained distance map. We provide results for different values of p = 1, p = 2 and p = ∞ at respectively first, second and third column. The last row illustrates distance computation on a Region Adjacency Graph (RAG), built from the original image. The first column shows the partition. Second one shows the associated graph where each node represents a region of the RAG and two nodes are linked if their regions are adjacent. Finally, the third column presents the distance map obtained from the top-left region. In that case, the weight function holds the similarity between two adjacent regions in the sense of mean intensity of each region, and the distance was computed with p = ∞.

Figure 4 presents illustrations of geodesic distances computed on an irregular 3d mesh and an irregular 2d graph with constant speed/potential. The two graphs are said irregular, because each vertex has a variable number of edges. In these cases, our results does not necessarily coincide with solutions of the Eikonal equation via Fast Marching algorithms. Indeed, the finite difference equation ( 31) used to compute these geodesics is not the Eikonal equation (but is adapted from) and only consider values of functions defined on graph vertices. There exists some algorithms, using geometric informations and local solvers designed for these particular cases, which better approximate the solution of the Eikonal equation on triangulated domains, as the Fast Marching algorithm for triangulated domains proposed by Kimmel [START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF], or Fast Sweeping on triangulated meshes by Qian [START_REF] Qian | Fast sweeping methods for eikonal equations on triangular meshes[END_REF].

In the following section, we will present the accuracy of the proposed fast algorithm to compute graph partitions.

Seeds

Distance map

-4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • c(u, v
) is a sequence of vertices (u 1 , ..., u j ) such as u = u 1 and v = u j with (u i , u i+1 ) ∈ E and i = 1, .., j -1. Let ρ : E → IR be a metric defined as

ρ(u, v) = min ρ∈C G (u,v) j-1 ∑ i=1 f (u i+1 ) -f (u i ) . ( 59
)
Given a set of K labels L = {l i } i=1,...,K and a set of K seeds S 0 = {s 0 i } i=1,...,K ⊆ V (where s 0 i is the seed of l i ), the energy ρ S : V → IR m induced by the metric ρ for all seed of S 0 is

ρ S (u) = min s i ∈S ρ(s i , u)∀u ∈ V. ( 60
)
The region R i of a label l i , is the set of vertices which are closer to s 0 i than to any other seeds with respect to the metric ρ. R i can be defined as

R i = {u ∈ V : ρ(s 0 i , u) ρ(s 0 j , u), ∀ j = 1, ..., K.}, (61)
A partition of G, for a given set of labels (with seeds S 0 ) and a metric ρ, is finally the set R(L, ρ) = {R i , ∀l i ∈ L}. For a graph G, to find it's partition corresponds to seek a minimal cost path over G, and can be easily computed with a simple Dijkstra like algorithm.

This general formulation of graph partitioning according to a metric can be easily linked with our definition of label propagation for graph clustering (see Sec. 4.1). Indeed, such a partition can be processed using the new finite difference equation [START_REF] Jansen | Second Generation Wavelets and Applications[END_REF] and the new class of algorithms proposed in Algo. 2 or Algo. 3, since we know that the label propagated at each vertex u ∈ V is the label which seed is the nearest of u (see Sec. 4 for more details). Given a vertex u ∈ V with a label l i , we have the following relation

f (u) = ρ S (u) = ρ(u, s 0 i ) = min s 0 j ∈S 0 ρ(u, s 0 j ) . (62)
Moreover, graph partitioning is intrinsic to our algorithm, because of the label propagation according to the distance measure (or arrival time), which naturally cluster the graph. The resulting partition corresponds to a generalized Voronoï diagram. In the next paragraph, we present a methodology to compute graph partition using our equations and algorithms, we call supervertices.

From Super Pixels to Super Vertices. Initially developed by Ren and Malik [START_REF] Ren | Learning a classification model for segmentation[END_REF], superpixels are an efficient way to reduce image complexity by grouping pixels in a region map while preserving contours. With TurboPixels [START_REF] Levinshtein | Turbopixels: Fast superpixels using geometric flows[END_REF], Levinshtein et al. have proposed an implementation of superpixels in which an image partition is obtained by dilating a regular grid of seeds so as to adapt to local image structure, where the dilation is performed using a level-set approach.

In this paper, we propose to extend such approach of image Seeds Distance Supervertices simplification to weighted graphs, using our fast algorithm for multilabel propagation on graphs. Because this adaptation groups vertices instead of pixels, we name it supervertices. As TurboPixels, our approach uses a set of labels whose seeds are placed on initial data, but seeds dilation becomes a label propagation which is controlled by our fast algorithm instead of iterative evolving equations. Figure 5 presents two examples of supervertices on 3D meshes. A review of some other mesh partitioning methods can be found in [START_REF] Shamir | A survey on Mesh Segmentation Techniques[END_REF] and references therein. In the case of image processing, the process is very similar to superpixels: the simplification is performed from a regular grid of distinct labels whose seeds position can be perturbed in the direction of the descending image gradient to avoid placing seeds close to string boundaries. This is illustrated in Fig. 6 on two images from the Berkeley Database1 . Such partition can be easily transformed in a Region Adjacency Graph (RAG), which can be used as a simplified version of the initial image that preserves texture information and strong boundaries. Then, it becomes very interesting to perform any graph-based algorithm on such a reduced version of the image. This is illustrated in Section 5.4.4 with active contour and semisupervised image segmentation using RAG.

Supervertices partition of images

Fig. 6 Image partition with supervertices. Both images were partitioned using a regular grid of seeds (perturbed according to the image gradient). See text for details.

Semi-Supervised Image segmentation and Data Clustering

In this Section, we present the behavior and efficiency of our algorithm for semi-supervised image segmentation and data clustering, using graph representation of these data. First we give qualitative and quantitative comparison with a recent and efficient graph-based approach for semi-supervised image segmentation. Then, we illustrate the algorithm behavior with local or non-local configurations and with non-negative speeds or either positive or negative speeds, through several examples involving different kinds of data.

In the case of image semi-supervised segmentation, graph-based approaches have became very popular in recent years. Many graph-based algorithms for image segmentation have been proposed, such as graph cuts [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in n-d images[END_REF], random walk [START_REF] Grady | Random walks for image segmentation[END_REF], shortest-path, watershed or frameworks that unify some of the previous methods (as powerwatershed) [START_REF] Couprie | Power Watersheds: A Unifying Graph-Based Optimization Framework IEEE Trans[END_REF][START_REF] Sinop | A seeded image segmentation framework unifying graph cuts and random walker which yields a new algorithm[END_REF]. Recently, these algorithms were all placed into a common framework [START_REF] Couprie | Power Watersheds: A Unifying Graph-Based Optimization Framework IEEE Trans[END_REF] that allows them to be seen as special cases of a single general semi-supervised algorithms, and leads to the following general formulation: Given a graph G = (V, E, w), the general algorithm consists in finding the func-tion f which minimizes the following energy function:

E p,q ( f ) = ∑ u∼v∈E w p uv | f (u) -f (v)| q (63)
where f represents the target configuration. Using such a formulation in the case of two classes A and B such that V A ∪ V B = V and V A ∩ V B = / 0, the graph clustering can be performed as follow:

1. The first step consists in finding an optimum f for (63), which is obtained by solving the equation

     f (A) = 1 f (B) = 0 f = arg min f E p,q ( f ) (64)
Such formulation can be interpreted using PdEs as arg min ∑ (∇ w f )(u) with w depending on p and q, and A and B represents the internal Dirichlet boundary condition. 2. Finally, the segmentation is performed via a threshold and the labeling function L : V → {1, 0} is given by

L(u) = 0 if f (u) < 1/2 1 if f (u) 1/2 (65)
In the case of more than two classes, the method needs to find as many optimum f i as distinct labels, where the shape is given by the set of vertices of the label and the background by the set of vertices of all other labels. Then, the labeling function is given by:

L(u) = arg max i f i (u) (66)
Another approach [START_REF] Bai | Geodesic matting: A framework for fast interactive image and video segmentation and matting[END_REF][START_REF] Falcão | The image foresting transform: Theory, algorithms, and applications[END_REF] to perform a graph clustering consists in computing a graph partition from the set of user's seeds and a metric, as it is described in the previous Section.

In this paper, we consider this approach which can be easily performed using Algo. 2 and Algo. 3.

Comparison with power watershed

In this paragraph, we propose qualitative and quantitative comparisons between our method and the recent power watershed [START_REF] Couprie | Power Watersheds: A Unifying Graph-Based Optimization Framework IEEE Trans[END_REF] which is one of the most efficient methods proposed in the previously cited framework. Results for the power watershed are obtained using the source code from Couprie's website2 Figure 7 presents a qualitative comparison between the two algorithms for semi-supervised image segmentation with two labels (shape and background). The comparison is performed on images from the Microsoft Grabcut database 3 .

The initial set of seeds is replaced by an eroded version in order to penalize both algorithms (as all propagation algorithms are very efficient when desired boundaries and seeds are very close). Results are presented for values of p ∈ {1, 2, ∞} and the power watershed (given by (63) with p = 2, q = 2). The mean computation time is given for each method (although it is quite difficult to compare runtime from two different implementations) and one can remark that the runtime difference between each p-norm is negligible, while the case p = 2 outperform the three other methods (p = 1, ∞ and power watershed), with smoother results.

On the contrary to others methods, the proposed algorithm complexity and computation time does not depend of the number of labels (neither seeds). Indeed, whatever the number of labels, each vertex is activated exactly once by the algorithm (and we recall that the arrival time and the label of each vertex are fixed when the vertex is activated, see Sec.4.1). This is illustrated on Fig. 8, which presents the runtime variation of both methods, in function of an increasing number of labels. The labels propagation is performed on a cytological slide photography (size 1280 × 2960 ) on which we consider a label for the background and a label per cell. For each experiment, an increasing number of label's seeds are placed on the image, and the non-seeded cells are considered as background. Figure 8 presents the whole set of labels and the slide segmentation with all labels for both two methods. One can remark that where the power watershed runtime grows as the number of labels, the proposed method runtime is relatively constant (the little increase is implementation dependent and is not a consequence of the algorithm).

Comparison in function of the number of distinct labels Elapsed time (seconds) Number of labels Fig. 9 Runtime variation in function of the number of labels. Variation is given for the power watershed (in red) and the proposed approach (in blue). The graph illustrates the runtime variation while the number of labels increase from 20 to 240. See text for a comment. (Color online).

Non-local image segmentation

Our approach using graphs has the advantage to naturally enable local and non-local configurations in the same formulation. In this paragraph, we show the benefits of non-local schemes as compared to local ones for semi-supervised image segmentation, especially noisy images or images that contain fine and repetitive structures. In order to hold nonlocal configurations, the graph is construct as a k-grid graph, (i.e a pixel u is linked to every pixels in a k × k window centered on u), and each pixel is associated with a patch feature vector (where each vector can be seen as a patch of texture). Figure 10 presents several examples of textured image segmentation. First column presents initial image and seeds. Second column shows local results obtained with a 4-adjacency grid graph using color feature vectors. Finally, the last column shows non-local results with a 11 × 11 window and patchs of 3 × 3 . All results were computed with a potential function P = 1 and p = 2. These results demonstrate the benefits of non-local configurations especially for textured images, where classical methods fail to found correctly the desired object. In non-local configuration cases, the graph weights better capture fine and repetitive information contained in the image.

Active contour

In this paragraph, we illustrate the behavior of our algorithm to perform active contour on graphs using an adaptation of the Chan and Vese model [START_REF] Chan | Active contours without edges[END_REF]. Although many recent approaches outperform this model [START_REF] Boykov | An integral solution to surface evolution pdes via geo-cuts[END_REF][START_REF] El Zehiry | Graph cut optimization for the mumford-shah model[END_REF][START_REF] Grady | The piecewise smooth mumford-shah functional on an arbitrary graph[END_REF], our choice lies in that we are not interested in present a new active contour model on graph, but only illustrate the adaptivity of our finite difference equation and algorithms. Given a front Γ represented by the level-set function φ , the Chan and Vese active contours model (without regularization term) can be summarized by the following PDE

∂ φ (x,t) ∂t = ((F x -µ 1 (t)) 2 -((F x -µ 2 (t)) 2 ∇φ (x,t) , (67)
where F x is the feature vector of x, µ 1 (t) and µ 2 (t) are the mean outside the front Γ at time t, respectively the mean inside. With F (x,t) = ((F x -µ 1 (t)) 2 -((F x -µ 2 (t)) 2 such an equation ( 67) is equivalent to the level-set function defined in [START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF], which is the continuous version of our finite difference equation. Then, transposed on weighted graphs, the active contours model for a front Γ i can be written as

(∇ - w f )(u) p = P l i (u,t), (68)
with P l i (u,t) = 1/F l i (u,t) and

F l i (u,t) = ((F u -µ 1 (t, l i )) 2 -(F u -µ 2 (t, l i )) 2 .
Finally, the front evolution is performed using Algo. 3 with the advantages given by our formulation as the speed of the Fast Marching and the ability to perform several active contours evolution simultaneously. It is illustrated by Fig. 13, where different steps of the contours evolutions are presented. The graph is a 4-adjacency grid graph, and the 4 initial sets of labels (or in other words the initial position of contours) were deliberately roughly placed, to illustrate the behavior of the algorithm in such configuration. Figure 11 presents another example which involves 9 distinct fronts/contours. Here again, the seeds are deliberately roughly placed to illustrate the robustness of the method to imprecise initialization Remark 4 In Algo.3, in order to transcript the time dependence of speed, the speeds F l i are periodically updated.

Region Adjacency Graph

Another advantage of our graph-based formulation is that the proposed algorithm can be applied to any graph, and therefore any graph representing images. To illustrate such Initial labels Final Segmentation mentation algorithm is performed. Figure 14 illustrates nonlocal semi-supervised image segmentation using a RAG. In order to let the labels grow beyond local neighborhood, each vertex neighborhood is extended by a k-nearest neighborhood based on mean color value. The final graph is a RAG ∪ k-NNG. One can remark that every cells are extracted even those without initial marker. Another example, presented in Fig. 12, illustrates active contour on a RAG. Finally, using region based graphs has the following advantages.

-Fast computation. By reducing the number of graph vertices (approximately 98% of reduction in term of vertices as compared to the pixel based graph), the computing time is reduced due to the reduced number of data to consider. -Non-local segmentation. The first example shows nonlocal object segmentation, where only one object is initially marked and the others are found by our method even if the objects are not spatially close or connected -Minimal number of initial seeds. With an appropriate graph structure, the segmentation requires a minimal number of initial seeds.

Data clustering

Finally, this last experiment illustrates the adaptivity and behavior of the proposed algorithm for high-dimensional unorganized data clustering. The clustering is performed similarly to the non-local image segmentation case, using a set of seeds and a weighted graph. The data consists in a set of

Conclusion

In this paper, we have considered a new finite difference equation which is an adaptation of the Eikonal equation and defined on weighted graphs. We have shown that this equation has a unique solution and proposed numerical schemes for its resolution with different L p norms, with p ∈ {1, 2, ∞}.

Based on such an equation, we have proposed efficient algorithms for multiple labels propagation on graphs of arbitrary topology which enables numerous applications as graph clustering, geodesic distance computation or front tracking. This can be used on many domains such as images, meshes, data or any structure that can be represented by a weighted graph.

Algorithm 2

 2 Labels propagation with non-negative speed 0. List of variables S 0 : the set of seed vertices. A : the set of active vertices. NB : the set of vertices in the narrow band FA : the set of vertices said as far away lab : the label indicator function 1. Initialization:

3

 3

Algorithm 3

 3 Labels propagation with changing-sign speed 0. List of variables: S 0 : the set of seed vertices A : the set of active vertices. NB : the set of vertices in the narrow band FA : the set of vertices said as far away lab : the label indicator function ∂ + S 0 i : outer boundary of the set of seeds S i ∂ -S 0 i : inner boundary of the set of seeds S i F l i : speed function associated to label l i . 1. Initialization:

Fig. 3

 3 Fig.3Adaptive weighted distance computation with different p values, weight functions w and graph topologies G. Results represent color distance maps with iso-level sets. At top: original image. First, second and third columns: results with p=2, 1 and ∞. First row, the weight w 1 is constant. Second row, the weight w 2 is designed to penalize vertical diffusion. Third row, the weight w 3 holds original image informations. First rows, the graph is built on the whole image. On the last, the distance is computed on a Region Adjacency Graph (RAG) from a partition of the original image.

Fig. 5

 5 Fig. 5 Supervertices on meshes. Graphs used in this example are faces graphs, where each face of meshes is represented by a vertex and linked by an edge to any adjacent face. In both cases, the weight function is given by f = 1 and P = 1. First column presents initial meshes with superimposed random seeds. Second column shows the distance maps computed where red means near a seed and purple means far from a seed. Finally, last column shows the resulting partition in super vertices.

Fig. 7 2 Fig. 8

 728 Fig.7Qualitative comparison between power watershed and the proposed approach, using two images from the grabcut database. See text for more details.SeedsPower watershed Eikonal p = 2

Fig. 11

 11 Fig. 11 Another example of active contours on a natural image. The evolution of the 9 contours is performed simultaneously without lost of efficiency.

Fig. 12 Fig. 13 Fig. 14 Fig. 15

 12131415 Fig. 12 Active contours on a weighted RAG from one contour. (a) The image partition the RAG is built on. (b,c,d) Different steps of the contour evolution. Points represent vertices and white lines represent edges. The contour is given by the red points and blue points represent the vertices which lie inside the contour. (e) The final contour transposed on initial image. (Color online).

Figure 15 presents

 15 Figure 15 presents the graph with 4 initial sets of labels, and the final clustering.

  : IR m → IR m and erosions ε : IR m → IR m of a function f 0 : IR m → IR by structuring sets B={x : x p ≤1}, with the general PDEs

	that morphological gradients constitute
	numerical schemes for solving time dependent morphologi-
	cal dilation and erosion processes over graphs and therefore
	solve the time dependent Eikonal equation.
	Continuous scale morphology (see for instance in [7, 38]
	and references therein) defines flat dilations δ

http://www.eecs.berkeley.edu/Research/Projects/CS/vision/bsds/

http://www.esiee.fr/ coupriec/

http://research.microsoft.com/en-us/um/cambridge/projects/ visionimagevideoediting/segmentation/grabcut.htm

This work was supported under a doctoral grant of the Conseil Régional de Basse-Normandie and the Coeur et Cancer association.

Fig. 4 First row: geodesic distance on a 3d mesh from two distinct seeds with, at left the seeds, at right the rainbow distance map where red means near and blue means far. Second row: geodesic distance on an irregular graph in IR 2 with one seed.

Graph partitioning

In this Section, we present the accuracy of the proposed algorithm (Algo. 2) for graph partitioning using a distance function (i.e. a metric), with some examples of data simplification. The use of a distance function for graph partitioning is not new, and we recommend interested readers to refer to [START_REF] Arbelaez | A Metric Approach to Vector-Valued Image Segmentation[END_REF][START_REF] Bai | Geodesic matting: A framework for fast interactive image and video segmentation and matting[END_REF][START_REF] Falcão | The image foresting transform: Theory, algorithms, and applications[END_REF] an references therein for a more complete review on this topic.

First, we recall a definition of metric-based graph partition. Let G = (V, E, w) be a weighted graph and C G (u, v) be the set of paths connecting two vertices u, v ∈ V . A path