Manuel Mera

Andrey Polyakov

Wilfrid Perruquetti

Gang Zheng

Finite-time Attractive Ellipsoid Method Using Implicit Lyapunov Functions

A finite-time version, based on Implicit Lyapunov Functions (ILF), for the Attractive Ellipsoid Method (AEM) is developed. Based on this, a robust control scheme is presented to ensure finite-time convergence of the solutions of a chain of integrators with bounded output perturbations to a minimal ellipsoidal set. The control parameters are obtained by solving a minimization problem of the "size" of the ellipsoid subject to a set of Linear Matrix Inequalities (LMI's) constraints, and by applying the implicit function theorem. A numerical example is presented to support the implementability of these theoretical results.

I. INTRODUCTION

The modeling and control of dynamical systems subject to bounded perturbations is one of the most relevant problems in control theory, this is because it is unrealistic to expect actual dynamical systems to be free of noise, or that mathematical models represent perfectly, without uncertainties, real systems. For this reason the research done in the robust control area has been extensive.

The minimization of the effects due to external perturbations and uncertainties, on the performances of dynamical systems is a well-known area where several results have been obtained, and a considerable amount of different approaches have been developed, such as sliding-mode control [START_REF] Utkin | Sliding Modes in Control and Optimization[END_REF], which in general only works for matched perturbations; H ∞ Control [START_REF] Orlov | Advanced H ∞ Control[END_REF], which usually asks for a vanishing condition on the perturbations; and neural networks [START_REF] Haykin | Neural networks and learning machines[END_REF], which has some implementability issues. Among these a remarkable one is the Attractive Ellipsoid Method (AEM).

The main advantages of the AEM reside in the fact that it is applicable to systems with unmatched and non-vanishing perturbations, also the control gain is usually obtained by the solution of a linear minimization problem, this is a minimization problem with a linear objective function and linear constraints, this suggests an advantage on the computational aspect of the controller design and its implementability.

The history of this approach can be track down to the late 60's and early 70's, when the characterization of uncertain dynamics by ellipsoidal sets was first introduced in [START_REF] Schweppe | Recursive state estimation: Unknown but bounded errors and system inputs. Automatic Control[END_REF] and [START_REF] Bertsekas | Recursive state estimation for a set-membership description of uncertainty. Automatic Control[END_REF]. Then, the use of ellipsoids as estimations of sets guaranteed to contain a significant variable was further developed in [START_REF] Kurzhanskii | Control and observation under uncertainty[END_REF], [START_REF] Chernousko | State estimation for dynamic systems[END_REF], [START_REF] Polyak | Ellipsoidal parameter or state estimation under model uncertainty[END_REF]. The concept of the asymptotically attractive (invariant) ellipsoid as used in this paper was formalized in [START_REF] Usoro | Ellipsoidal set-theoretic control synthesis[END_REF], [START_REF] Polyak | Suppression of bounded exogenous disturbances: Output feedback[END_REF] for linear systems and later extended to nonlinear systems in [START_REF] Poznyak | Practical output feedback stabilisation for a class of continuous-time dynamic systems under sample-data outputs[END_REF], [START_REF] Mera | Robust control for a class of continuous-time dynamical systems with sample-data outputs[END_REF] and [START_REF] Poznyak | Attractive Ellipsoids in Robust Control[END_REF].

The AEM is based on [START_REF] Gonzalez-Garcia | Output linear controller for a class of nonlinear systems using the invariant ellipsoid technique[END_REF] the Lyapunov analysis, so it is natural to use the Implicit Lyapunov Function (ILF) method, originally presented in [START_REF] Korobov | A general approach to the solution of the bounded control synthesis problem in a controllability problem[END_REF], and later revisited in [START_REF] Adamy | Soft variable-structure controls: a survey[END_REF], [START_REF] Polyakov | Finitetime stabilization using implicit Lyapunov function technique[END_REF], to extend this approach and obtain additional features such as finite-time convergence to the ellipsoidal region. This is a very desirable property considering that many applications require control algorithms fast enough to guarantee the convergence of system dynamics to a desirable value or set of values in finite-time. The finite-time stabilization [START_REF] Moulay | Finite time stability and stabilization of a class of continuous systems[END_REF] is crucial when the transient dynamics are required to end as fast as possible or in a specific time interval, as in robotics and aeronautics. Some examples of this can be found in [START_REF] Bhat | Finite time stability of continuous autonomous systems[END_REF], [START_REF] Haimo | Finite time controllers[END_REF] and [START_REF] Orlov | Finite time stability and robust control synthesis of uncertain switched systems[END_REF].

The ILF, as presented in [START_REF] Polyakov | Finitetime stabilization using implicit Lyapunov function technique[END_REF], [START_REF] Polyakov | Homogeneous differentiator design using implicit Lyapunov function method[END_REF], [START_REF] Polyakov | Finite-time and fixed-time stabilization: Implicit Lyapunov function approach[END_REF], is used to construct a robust control strategy by using a Lyapunov function defined implicitly in an algebraic equation. The stability analysis then does not require the explicit solution of this equation, instead the stability conditions can be revised directly using the implicit function theorem (see [START_REF] Courant | Introduction to Calculus and Analysis[END_REF]).

To simplify this analysis an homogeneous function is proposed. This is because homogeneity is a very useful property in the study of finite-time stability (see [START_REF] Zubov | Methods of A.M. Lyapunov and Their Applications[END_REF], [START_REF] Levant | Homogeneity approach to high-order sliding mode design[END_REF], [START_REF] Perruquetti | Finitetime observers: application to secure communication[END_REF]). Specifically, if an asymptotically stable system is homogeneous of negative degree, then it is finite-time stable.

The AEM with the ILF derives in a minimization problem of the size (trace) of the ellipsoid characterizing matrix subject to a set of constraints obtained from the finite-time stability conditions, which can be expressed as a set of Linear Matrix Inequalities (LMI's). The implementable control strategy is obtained from the solution of this optimization problem, and the application of a numerical procedure.

II. PROBLEM STATEMENT Considering the system

ẋ(t) = Ax(t) + bu(t), y(t) = x(t) + ξ (t), (1)
where x ∈ R n is the state vector, u ∈ R is the control input, y ∈ R n is the measurable output, ξ ∈ R n is an unknown but bounded perturbation,

ξ (t) T Qξ (t) ≤ 1, ∀t ≥ 0, (2)
with Q = Q T > 0, no additional restriction is imposed on ξ , however it is not considered to be random,

A =      0 1 • • • 0 . . . . . . . . . . . . 0 0 • • • 1 0 0 • • • 0      , and b =      0 . . . 0 1     
The objective of this paper is to design a robust feedback control strategy for the system (1) that ensures finite-time convergence of the system solutions to an ellipsoidal region despite having noisy measurements.

The ILF approach is used to develop an extended version of the AEM which includes finite-time convergence to the ellipsoid. The main tool for these two approaches is the Lyapunov function method, which in this case is implicitly defined as a solution of an algebraic equation. The analysis of the ellipsoidal set attractivity does not require the solution of this equation, because it is possible to use the implicit function theorem to check all the required conditions directly from the implicit formulation.

III. PRELIMINARIES

A. Finite-Time Attractive Ellipsoid

Considering the system of the form

ẋ(t) = f (x(t), ξ (t)), x(0) = x 0 , ξ (0) = ξ 0 ; (3)
where

• x(t) ∈ R n is the state vector, • ξ (t) ∈ R n is an unknown but bounded perturbation vector ξ (t) ≤ l 0 , ∀ t ≥ 0, l 0 ∈ R + ,
• f : R n → R n is a nonlinear continuous vector field.

Assume that the following ellipsoid

E x := x ∈ R n : x T P-1 x ≤ 1 , P = PT > 0. (4)
is an invariant set of (3). Additionally we consider the notation, for any

θ ∈ R n θ E x := d(θ , E x ) = inf η∈E x d(θ , η),
for distance from a point θ to a set E x . Definition 1: [START_REF] Lasalle | The Stability of Dynamical Systems[END_REF] The set E x , is an asymptotically attractive set for the system (3

) if x(t, x 0 ) E x → 0, as t → ∞, for any x 0 ∈ R n .
The asymptotic attractivity of E x guarantees convergence of the system solutions, with any initial condition x 0 , to the ellipsoid. However, it is always desirable to have some estimation for the time of convergence to the set. Addressing this, the definition of finite-time attractive ellipsoid is presented.

Definition 2:

The set E x is a finite-time attractive ellipsoid of system (3) if it is an asymptotically attractive ellipsoid, and there exists a function T (called the settling time function): R n \E x → R + such that for any x 0 / ∈ E x (i.e.

x 0 ∈ R n : x T 0 P-1 x 0 > 1) lim t→T (x 0 ) x(t, x 0 ) E x = 0,
in terms of the matrix P this means lim t→T (x 0 )

x T (t, x 0 ) P-1 x(t, x 0 ) ≤ 1.

Note that the Definitions 1 and 2 do not ask for stability of the set, only consider the attractivity property.

B. Implicit Lyapunov Method for Attractive Ellipsoid

For the very specific case when ξ (t) = 0, ∀t, we have: Theorem 1 ( [START_REF] Polyakov | Homogeneous differentiator design using implicit Lyapunov function method[END_REF]): If there exists a continuous function

G : R + × R n → R (V, x) → G(V, x)
satisfying the conditions C1) G is continuously differentiable outside the origin for all positive V ∈ R + and for all x ∈ R n \{0}; C2) for any x ∈ R n \{0} there exists

V ∈ R + such that G(V, x) = 0; C3) let Φ = {(V, x) ∈ R + × R n \{0} : G(V, x) = 0}, then, ∀(V, x) ∈ Φ lim x →0 V = 0 + , lim V →0 x = 0, lim x →∞ V = +∞; C4) the inequality ∂ G(V,x) ∂V < 0 holds for all V ∈ R + and for all x ∈ R n \ {0}; C5) ∂ G(V,x) ∂ x f (x, 0) < 0 holds for all (V, x) ∈ Φ;
then the origin of ( 3) is globally uniformly asymptotically stable.

Theorem 2 ([21]):

If there exists a continuous function G : R + × R n → R that satisfies the conditions C1)-C4) of Theorem 1 and the condition C5bis) there exist c > 0 and 0 < µ < 1 such that

∂ G(V,x) ∂ x f (x, 0) ≤ cV 1-µ ∂ G(V,x)
∂V for (V, x) ∈ Φ, then the origin of the system (3) is globally uniformly finite-time stable and T (x 0 ) ≤ V µ 0 cµ is the settling time function, where G(V 0 , x 0 ) = 0. Now considering the more general case when ξ (t) ≤ l 0 ∀ t in (3), it is unrealistic to expect that the system trajectories converge to the origin. Still, it is possible for the system trajectories to converge to a certain set depending on l 0 , the bound of the perturbations. In order to define this set in terms of the conditions given by the ILF we present the next Corollary.

Corollary 1: Let the conditions C1)-C4) from Theorem 1 hold. If the inequality ∂ G(V,x) ∂ x f (x, ξ ) < 0 holds for all (V, x) ∈ Φ and all ξ ≤ l 0 such that V (x) > 1, then the set

W := {x ∈ R n : (V, x) ∈ Φ, 0 < V ≤ 1} ∪ {0},
is a globally asymptotically stable set of (3).

These conditions repeat in an implicit form the requirements for the Lyapunov-like (storage) function of the AEM, this is a positive function which is not necessarily monotonically non-increasing in all R n , but only outside of a bounded set. Conditions C1, C2, C4 and the implicit function theorem [START_REF] Courant | Introduction to Calculus and Analysis[END_REF], denote that the equation G(V, x) = 0 defines a unique function V : R n \{0} → R + for all x ∈ R n \{0}. The function V is continuously differentiable outside the origin and ∂V

∂ x = -∂ G(V,x) ∂V -1 ∂ G(V,x) ∂ x
for G(V, x) = 0, ∀x = 0. The implicit function theorem and conditions C3, C4 and Corollary 1 prove that V ≤ 0, while V > 1.

The next corollary specifies the finite-time attractive set.

Corollary 2: If there exists a continuous function G : R + × R n that satisfies the conditions C1-C4, the one in Corollary 1, and additionally it fulfills that there exists c > 0 and µ ∈ (0, 1) such that ∂ G(V,x)

∂ x f (x, ξ ) ≤ -cV 1-µ ∂ G(V,x)
∂V for (V, x) ∈ Φ and V (x) > 1, then W is a globally finite-time attractive set for system (3), and

T (x 0 ) ≤ V µ 0 -1 cµ is the settling time estimation with G(V 0 , x 0 ) = 0.
One of the main features of the control design using the AEM is the possibility to minimize in certain sense the ellipsoid size, insuring that the system solutions converge to a minimal set in R n regardless of perturbations or uncertainties [START_REF] Poznyak | Attractive Ellipsoids in Robust Control[END_REF]. The usual characterization of the size of the ellipsoid in the AEM is conveniently described by the trace of the matrix P. The "natural" restrictions of the system and its parameters are written in the form of LMI's. Hence, the control parameters that ensure convergence to this minimal set are obtained through the solution of a linear minimization problem which objective function is tr( P), and the corresponding constraints are a set of LMI's. The same approach is used in the next sections to obtain a characterization of the minimal finitetime version of the Attractive Ellipsoid.

IV. CONTROL DESIGN

A. Finite-time attractive set

Let us consider the following ILF [START_REF] Polyakov | Finitetime stabilization using implicit Lyapunov function technique[END_REF]

G(V, x) := x T D r (V -1 )PD r (V -1 )x -1, (5)
where

V ∈ R + , x ∈ R n , P = P T ∈ R n×n , P > 0, D r (λ ) := diag{λ r i },
with

r i = 1 + (n -i)µ, i ∈ {1, . . . , n}, 0 < µ < 1.
Denoting

H µ := diag{r i },
and defining V y := V (y).

Note that when µ = 0 equation ( 5) coincides with the storage function in the AEM.

Theorem 3: Let µ ∈ (0, 1), and some tunning parameters δ and β such that δ > 1, β ∈ R + , the matrix

Ω =       AX + bY + XA T +Y T b T +H µ X + XH µ bY bY X Dm Y T b T -γ -2+2µ δ X 0 0 Y T b T 0 -β X 0 Dm X 0 0 -r 1 β I n      
, where X = X T ∈ R n×n , Y ∈ R 1×n and Dm ∈ R n×n is a diagonal matrix which its i-th entry m i is given by

m i := max{|γ -(n-i+1)µ -1|, |γ -(n-i+1)µ -1|},
for i ∈ {1, . . . , n}, and

γ = r 1 r 2 1 - r 1 r 3 , γ = r 2 r 1 1 + r 1 r 3 ,
with r 1 := 1/λ max (X), r 2 := 1/λ min (X), and r 3 := λ min (Q). If the LMI system

         0 < Q -1 ≤ X 0 < δ 2 X ≤ XH µ + H µ X 0 < QH µ + H µ Q Ω ≤ 0 (6)
is feasible for some X and Y ; the control input

u(V y , y) = V 1-µ y kD r (V -1 y )y, (7)
where k := Y X -1 , V y ∈ R + is such that G(V y , y) = 0, with G defined as in [START_REF] Courant | Introduction to Calculus and Analysis[END_REF].

Then the ellipsoid (4), with P = D r (γ)XD r (γ) is a finitetime attractive ellipsoid for the closed-loop system (1) with the settling time function estimation given by T (x 0 ) ≤ 1 1-δ -1 γ µ (V µ (y 0 ) -1).

B. Minimization of the Attractive Ellipsoid

Let us consider the minimization of the ellipsoidal set described in Theorem 3. In order to do so we minimize the trace of the matrix X subject to the constraints defined in the LMI's system, as usual. Note that the minimization of X is equivalent to the minimization of the ellipsoidal set characterized by D r (γ -1 )PD r (γ -1 ). Then, the minimization problem statement is the following tr(X) → min X,Y [START_REF] Haykin | Neural networks and learning machines[END_REF] subject to (6)

C. On the boundedness of the control

If it is necessary to restrict the magnitude of the control signal inside the ellipsoid, it is sufficient to include an additional constraint σ -2 k T k ≤ P to the optimization problem [START_REF] Haykin | Neural networks and learning machines[END_REF]. This, in terms of the matrices X and Y and using the Schur's complement can be represented in LMI form as X Y T Y σ 2 ≥ 0, this inequality implies that for any trajectory inside the ellipsoid (V y ≤ 1) we will have

u 2 = V 2-2µ y y T D r (V -1 y )k T kD r (V -1 y )y ≤ σ 2 V 2-2µ y y T D r (V -1 y )PD r (V -1 y )y ≤ σ 2 V 2-2µ y ≤ σ 2 .
V. NUMERICAL RESULTS Considering a chain of integrators as in (1) with A ∈ R 4×4 , with the initial conditions x 0 = 1, 1, 1. 1 T the selected perturbation considered was ξ (t) = 0.5 + 0.25 sin(100t) + 0.25 sin(t) 0.0002, 0.002, 0.02 0.2 T , and the matrix Q is defined as [START_REF] Polyakov | Finitetime stabilization using implicit Lyapunov function technique[END_REF], which is based on the bisection numerical method, to find the zeros of (5). The simulations have been done using the explicit Euler method with a sampling period with 0.001 seconds. The obtained approximation for the settling time was T (x 0 ) ≤ 18.5037 seconds. Figure 1 and figure 2 show the states convergence to the ellipsoidal set in the x 1 vs. x 2 and x 3 vs x 4 projections respectively. Figure 3 presents the evolution of the states in time, while figure 4 shows the control output in the same way.

Q =     24

VI. CONCLUSIONS

A theoretical extension of the well known Attractive Ellipsoid Method was developed in this paper using the Implicit Lyapunov Function approach. Finite-time convergence despite noisy measurements to an ellipsoidal set is proved. This is achieved using a control input obtained from an ILF and the solution of a linear, respect to the matrix variables, minimization problem. The implementable numerical algorithm, based on the above extension is validated through a numerical example.

  Selecting µ = 0.4, δ = 1.5, σ 2 = 125 and β = 1, the ellipsoidal matrix and the matrix gain obtained from the solution of the LMI problem stated in theorem 3-0.0492 0.0342 -0.0004 -0.0492 0.0680 -0.0658 0.0179 0.0342 -0.0658 0.1005 -0.0834 -0.0004 0.0179 -0.0834 0.2803     , k = -111.6465 -121.0438 -52.2740 -10.3568 , and γ = 14.9111. The function V is calculated by the algorithm presented in

Fig. 1 .Fig. 2 . 4 Fig. 3 .

 1243 Fig.1. System trajectories projection on the subspace x 1 vs. x 2 .

Fig. 4 .

 4 Fig. 4. Control Vs. Time.


Vaïa Machairas
email: vaia.machairas@mines-paristech.fr

Matthieu Faessel
email: matthieu.faessel@mines-paristech.fr

David Cárdenas-Peña

Théodore Chabardes
email: theodore.chabardes@mines-paristech.fr

Thomas Walter
email: thomas.walter@mines-paristech.fr

Etienne Decencière
email: etienne.decenciere@mines-paristech.fr

Dr Yonggang V Shi

M Machairas

T Faessel

E Chabardes

Decencière

IEEE TRANSACTIONS ON IMAGE PROCESSING 1 Waterpixels

Keywords: Superpixels, watershed, segmentation. 17 I AQ:2 AQ:3 27

Many approaches for image segmentation rely on a 1 first low-level segmentation step, where an image is partitioned 2 into homogeneous regions with enforced regularity and adherence 3 to object boundaries. Methods to generate these superpixels have 4 gained substantial interest in the last few years, but only a few 5 have made it into applications in practice, in particular because 6 the requirements on the processing time are essential but are not 7 met by most of them. Here, we propose waterpixels as a general 8 strategy for generating superpixels which relies on the marker 9 controlled watershed transformation. We introduce a spatially 10 regularized gradient to achieve a tunable tradeoff between the 11 superpixel regularity and the adherence to object boundaries. 12 The complexity of the resulting methods is linear with respect 13 to the number of image pixels. We quantitatively evaluate our 14 approach on the Berkeley segmentation database and compare 15 it against the state-of-the-art.

16

2) connected partition: each SP is made of a single connected component and the SPs constitute a partition of the image;

3) adherence to object boundaries: object boundaries should be included in SP boundaries;

4) regularity: SPs should form a regular pattern on the image. This property is often desirable as it makes the SP more convenient to use for subsequent analysis steps. The requirements on regularity and boundary adherence are to a certain extent oppositional, and a good solution typically aims at finding a compromise between these two requirements.

In addition to these requirements on superpixel quality, computational efficiency is an absolutely essential aspect, as the partition into superpixels is typically only the first step of an often complex and potentially time consuming workflow. Methods of linear complexity are consequently of particular interest.

We therefore hypothesized that the Watershed transformation [1], [2] should be an interesting candidate for superpixel generation, as it has been shown to achieve state-of-the-art performance in many segmentation problems, it is non-parametric, and there exist linear-complexity algorithms to compute it, as well as efficient implementations [3], [4]. The only often cited drawback, oversegmentation, does not seem to be problematic for superpixel generation, as long as we can control the degree of oversegmentation (number of superpixels), and the regularity of the resulting partition.

Given these considerations, we propose a strategy for applying the watershed transform to superpixel generation, where we use a spatially regularized gradient to achieve a tunable trade-off between superpixel regularity and adherence to object boundaries. We quantitatively evaluate our method on the Berkeley segmentation database and show that we outperform the best linear-time state-of-the art method: Simple Linear Iterative Clustering (SLIC) [5]. We call the resulting superpixels "waterpixels." is spent on the computation of the SP. Type B methods may spend more time on finding appropriate seeds, but can therefore afford not to iterate the SP generation.

B. Building Superpixels From Seeds

In the second step, the partition into superpixels is built from the seeds. Among the methods with linear complexity, there are two main strategies for this: Shortest Path Methods (Type 1) [11], [13]: these methods are based on region growing: they start from a set of seeds (points or regions) and successively extend them by incorporating pixels in their neighborhood according to a usually image dependent cost function until every pixel of the image plane has been assigned to exactly one superpixel. This process may or may not be iterated.

Shortest Distance Methods (Type 2) [5], [12]: these are iterative procedures inspired by the field of unsupervised learning, where at each iteration step, seeds (such as centroids) are calculated from the previous partition and pixels are then re-assigned to the closest seed (like for example the k-means approach).

Even though methods inspired by general clustering methods (type 2) seem appealing at first sight, in particular when they globally optimize a cost function, this class of methods does not guarantee connectivity of the superpixels for arbitrary choices of the pixel-seed distance (see [5], [12]). For instance, the distance metric proposed in [5] (a combination of Euclidean and grey level distance), leads to non-connected superpixels, which is undesirable. To solve this issue, a postprocessing step is necessary, consisting either in relabeling the image so that every connected component has its own label (see [12]), leading to a more irregular distribution of SP sizes and shapes, or in reassigning isolated regions to the closest and large enough Superpixel, as in [5], leading to non-optimality of the solution and an unpredictable number of superpixels. In addition, such postprocessing increases the computational cost and can turn out to be the most time-consuming step when the image contains numerous small objects/details compared to the size of the Superpixel.

On the contrary, methods based on region growing (type 1) inherently implement a "path-type" distance, where the distance between two pixels does not only depend on value and position of the pixels themselves, but on values and positions along the path connecting them. Type 1 methods imply connected superpixel regions, for which the number of superpixels is exactly the number of seeds.

C. Other Properties

It is generally accepted that a good superpixel-generation method should provide to the user total control over the number of resulting Superpixels. While this property is achieved by [11]- [14], some only reach approximatively this number because of post-processing (either by splitting too big superpixels, or removing small isolated superpixels as in [5]). Another parameter is the control on superpixels regularity in the trade-off between regularity and adherence to contours. Only [5] and [12] enable the user to weight the importance Let f : D → V be an image, where D is a rectangular subset of Z 2 , and V a set of values, typically {0, . . . , 255} when f is a grey level image, or {0, . . . , 255} 3 for color images.

The first step consists in computing the gradient image g of the image f . The choice of the gradient operator depends on the image type, e.g. for grey level images we might choose a morphological gradient. This gradient will be used to choose the seeds (section III-B) and to build the regularised gradient (III-C).

For the definition of cells, we first choose a set of N points {o i } 1≤i≤N in D, called cell centers, so that they are placed on the vertices of a regular grid (a square or hexagonal one for example). Given a distance d on D, we denote by σ the grid step, i.e. the distance between closest grid points.

A Voronoi tesselation allows to associate to each o i a Voronoi cell . For each such cell, a homothety centered on o i with factor ρ (0 < ρ ≤ 1) leads to the computation of the final cell C i . This last step allows for the creation of a margin between neighbouring cells, in order to avoid the selection of markers too close from each other.

B. Selection of the Markers

As each cell is meant to correspond to the generation of a unique waterpixel, our method, through the choice of one marker per cell, offers total control over the number of SP, with a strong impact on their size and shape if desired.

First, we compute the minima of the gradient g. Each minimum is a connected component, composed of one or more pixels. These minima are truncated along the grid, i.e. pixels which fall on the margins between cells are removed.

Second, every cell of the grid serves to define a region of interest in the gradient image. The content of g in this very region is then analyzed to select a unique marker, as explained in the next paragraph.

For each cell, the corresponding marker is chosen among the minima of g which are present in this very cell.

If several minima are present, then the one with the highest surface extinction value [19] is used. We have found surface extinction values to give the best performances compared with volume and dynamic extinction values (data not shown).

It may happen that there is no minimum in a cell. This is an uncommon situation in natural images. In such cases, we must add a marker for the cell which is not a minimum of g, in order to keep regularity. One solution could be to simply choose the center of the cell; however, if this point falls on a local maximum of the gradient g, the resulting SP may coincide with the maximum region and therefore be small in size (leading to a larger variability in size of the SP). We propose instead to take, as marker, the flat zone with minimum value of the gradient inside this very cell.

In both cases (i.e. either there exists at least one minimum in the cell or there is not), the selected marker has to be composed of a unique connected component to ensure regularity and connectivity of the resulting superpixel. However, it might not be the case, respectively if more than one minimum

288 ∀ p ∈ D, d Q ( p) = 2 σ min i∈[1,N] d( p, q i ) ( 1
)
where σ is the grid step defined in the previous section. The normalization by σ is introduced to make the regularization independent from the chosen SP size.

We have studied two possible choices of the q i . The first one is to choose them equal to the markers: q i = M i . Resulting waterpixels are called m-waterpixels. The second one consists in setting them at the cell centers: q i = o i , which leads to c-waterpixels. We have found that the first gives the best adherence to object boundaries, while the second produces more regular superpixels.

The spatially regularized gradient g reg is defined as follows:

g reg = g + kd Q (2)
where g is the gradient of the image f , d Q is the distance function defined above and k is the spatial regularization parameter, which takes its values within + . The choice of k is application dependent: when k equals zero, no regularization of the gradient is applied; when k → ∞, we approach the Voronoi tessellation of the set {q i } 1≥i≥N in the spatial domain.

In the final step, we apply the watershed transformation on the spatially regularized gradient g reg , starting the flooding from the markers {M i } 1≤i≤N , so that an image partition {s i } 1≤i≤N is obtained. The s i are the resulting waterpixels.

IV. EXPERIMENTS

In order to evaluate waterpixels, the proposed method has been applied on the Berkeley segmentation database [20] and benchmarked against the state-of-the-art. This database is divided into three subsets, "train", "test" and "val", containing respectively 200, 200 and 100 images of sizes 321 × 481 or 481 × 321 pixels. Approximately 6 human-annotated ground-truth segmentations are given for each image. These ground-truth images correspond to manually drawn contours.

A. Implementation

We have found that it is beneficial to pre-process the images from the database using an area opening followed by an area closing, both of size σ 2 /16 (where σ is the chosen step size of the regular grid). This operation efficiently removes details which are clearly smaller than the expected waterpixel area and which should therefore not give rise to a superpixel contour.

The Lab-gradient is adopted here in order to best reflect our visual perception of color differences and hence the pertinence of detected objects. The margin parameter ρ, described in III-A, is set to 2 3 .

The cell centers correspond to the vertices of a square or an hexagonal grid of step σ . The grid is computed in one pass over the image, by first calculating analytically the coordinates of the set of pixels belonging to each cell and then assigning to them the label of their corresponding cell. We will display the results for the hexagonal grid, as hexagons are more isotropic than squares. Interestingly, they also lead to a better quantitative performance, which was intuitively expected.

The implementation of the waterpixels was done using the Simple Morphological Image Library (SMIL) [21]. SMIL is a Mathematical Morphology library that aims to be fast, lightweight and portable. It brings most classical morphological operators re-designed in order to take advantage of recent computer features (SIMD, parallel processing, …) to allow handling of very large images and real time processing. As a gradient-based approach, the quality of the watershed is dependant on the borders contrast. If we look at the contours of objects missed by waterpixels, we see that it is due to the 368 weakness of the gradient, as illustrated in Figure 4.

B. Qualitative Analysis

369

C. Evaluation Criteria

370

SP methods produce an image partition {s i } 1≤i≤N . In order 371 to compute the SP borders, we use a morphological gradient 372 with a 4 neighborhood. Note that the resulting contours are 373 two pixels wide. To this set S c , we add the one pixel wide 374 image borders S b . The final set is denoted C. The ground 375 truth image corresponding to the contours of the objects to be 376 segmented, provided in the Berkeley segmentation database, 377 is called GT .

378

In superpixel generation, we look for an image decomposi-379 tion into regular regions that adhere well to object boundaries. 380 We propose to use three measures to evaluate this trade-381 off, namely boundary-recall, contour density and average 382 mismatch factor, as well as computation time.

383

There are two levels of regularity: (1) the number of pixels 384 required to describe the SP contours, which can be seen as a 385 measure of complexity of individual SP, and (2) the similarity 386 in size and shape between SP.

387

The first property is evaluated by the Contour Density, 388 which is defined as the number of SP contour pixels divided 389 by the total number of pixels in the image:

390 C D = 1 2 |S c | + |S b | |D| (3) 391
Note that |S c | is divided by 2 since contours are 392 two-pixel-wide.

393

The second property, i.e. similarity in size and shape, is 394 evaluated by an adapted version of the mismatch factor [22]. 395 The mismatch factor measures the shape and size dissimilarity 396 between two regions. Given two sets, A and B, the mismatch 397 factor m f between them is defined as:

398 m f (A, B) : = |A ∪ B \ A ∩ B| |A ∪ B| 399 = 1 - |A ∩ B| |A ∪ B| (4) 400
The mismatch factor and the Jaccard index thus sum to 401 one. Aiming to measure the superpixel regularity, we adapted 402 the mismatch factor to estimate the spread of size and shape 403 distribution. Hence, the average mismatch factor M F is 404 proposed as:

405 M F = 1 N N i=1 m f (s * i , s * ) (5) 406
where s * i is the centered version of superpixel s i , and s * is 407 the average centered shape of all superpixels. The complete 408 definition of the average mismatch factor is given in Appendix. 409

Note that although compactness is sometimes used in super-410 pixels evaluation (see [23]), it is a poor measurement for 411 region regularity. For example, perfectly-rectangular regions 412 are regular but not compact (because they are different from 413 discs). Waterpixels can in principle tend towards differently 414 shaped superpixels (rectangles, hexagons or other), depending 415 on the grid and the regularization function used. Since the 416 average mismatch factor compares each superpixel against an 417

424 B R = |{ p ∈ GT, d( p, C) < 3}| |GT | ( 6
)
425
where d is the L 1 (or Manhattan) distance.

426 While precision cannot be directly used in the context of over-segmentations, boundary-recall has to be, in this particular case of superpixels, interpreted with caution. Indeed, as noted also by Kalinin and Sirota [24], very tortuous contours systematically lead to better performances: because of their higher number, SP contour pixels have a higher chance of matching a true contour, increasing artificially the boundaryrecall. Hence, we propose to always consider the trade-off between boundary-recall and contour density to properly evaluate the adherence to object boundaries, penalizing at the same time the cost in pixels to describe SP contours.

D. Quantitative Analysis and Comparison With State-of-the-Art

In this paragraph, we will use m-waterpixels and denote them directly as "waterpixels" for the sake of simplicity.

During the design of the algorithm, we used intermediate results from the train and test subsets of the Berkeley database. Therefore, we report the results obtained for the validation subset ("val"), which contains 100 images. Results for boundary-recall, average mismatch factor and contour density are averaged for this subset and shown in Figure 5. than for SLIC (k = 15). This means that for the same number 474 of superpixels, waterpixels contours are shorter than SLIC 475 contours, which is partly explained by less tortuous contours. 476 We are currently working on a new implementation of minima computation/evaluation which would be less dependent on the number of superpixels.

To conclude this section, waterpixels are generally faster to compute than SLIC superpixels, and they are at least as performant in the trade-off between adherence to object boundaries and regularity in shape and size, while using much less pixels to describe their contours.

V. DISCUSSION AND PERSPECTIVES

We have shown that waterpixels produce competitive results with respect to the state-of-the-art. These advantages are valuable in the classification/detection/segmentation pipeline, where superpixels play the part of primitives. Moreover, there is one major difference in the construction of the algorithm: the SLIC approach does not impose any connectivity constraint. The resulting superpixels are therefore not necessarily connected, which requires some ad hoc postprocessing step. In contrast, waterpixels are connected by definition, and the connectivity constraint is actually implemented in the distance used.

The proposed approach is gradient-based. Standard methods can be used to compute this gradient, or a specific gradient computation method can be designed for a given application. In any case, this offers flexibility to waterpixels. One limitation though is the quality of the signal in such a gradient image. As seen in 4, alteration by noise or insufficiently contrasted contours may lead to the prevalence of regularity over adherence to object boundaries. If filtering steps are usually enough to deal with noise and remove non pertinent small details, parameter values have to be optimized for each database. Future work will aim at overcoming this limitation by adding a learning step of optimal filtering values for specific databases.

The general design of waterpixels offers many prospects. Among them, one promising field of improvement resides in the placement of markers, as they constitute the main degree of freedom of the method. We are currently investigating the possibility to select the markers in an optimal manner, for example by formulating the marker placement as a p-dispersion problem (see [25]) in an augmented space.

The speed of waterpixels contributes to expanding their s * i of s i is obtained by translating s i so that its barycenter is 586 the origin of the coordinates system.

587

The average shape s * of the {s i } is computed as follows.

588

Let first define function S:

589 S : D -→ N x p -→ N i=1 1 i (x p ) (7

595

The average centered shape s * is then the set S t 0 , where t 0 596 is the maximal threshold value which enables s * to have an 597 area greater than or equal to μ A :

598 t 0 = max{t |S t | ≥ μ A } ( 8
) 599 s * = S t 0 (9) 600
Finally, the mean mismatch factor of superpixels {s i } 1≤i≤N is: V. Machairas, M. Faessel, T. Chabardes, and E. Decencière are with MINES ParisTech, Paris 75006, France, also with PSL Research University, Paris 75005, France, and also with the Center for Mathematical Morphology, Fontainebleau 77305, France (e-mail: vaia.machairas@mines-paristech.fr; matthieu.faessel@mines-paristech.fr; theodore.chabardes@mines-paristech.fr; etienne.decenciere@mines-paristech.fr). D. Cárdenas-Peña is with the Signal Processing and Recognition Group, Universidad Nacional de Colombia, Manizales 170001-17, Colombia (e-mail: dcardenasp@unal.edu.co).

M F = 1 N N i=1 m f (s * i , s * ). ( 10
T. Walter is with MINES ParisTech, Paris 75006, France, also with PSL Research University, Paris 75005, France, also with the Centre for Computational Biology, Fontainebleau 77305, France, also with the Institut Curie, Paris 75248, France, and also with Inserm, Paris 75248, France (e-mail: thomas.walter@mines-paristech.fr).

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TIP.2015.2451011 2) connected partition: each SP is made of a single connected component and the SPs constitute a partition of the image;

3) adherence to object boundaries: object boundaries should be included in SP boundaries;

4) regularity: SPs should form a regular pattern on the image. This property is often desirable as it makes the SP more convenient to use for subsequent analysis steps. The requirements on regularity and boundary adherence are to a certain extent oppositional, and a good solution typically aims at finding a compromise between these two requirements.

In addition to these requirements on superpixel quality, computational efficiency is an absolutely essential aspect, as the partition into superpixels is typically only the first step of an often complex and potentially time consuming workflow. Methods of linear complexity are consequently of particular interest.

We therefore hypothesized that the Watershed transformation [1], [2] should be an interesting candidate for superpixel generation, as it has been shown to achieve state-of-the-art performance in many segmentation problems, it is non-parametric, and there exist linear-complexity algorithms to compute it, as well as efficient implementations [3], [4]. The only often cited drawback, oversegmentation, does not seem to be problematic for superpixel generation, as long as we can control the degree of oversegmentation (number of superpixels), and the regularity of the resulting partition.

Given these considerations, we propose a strategy for applying the watershed transform to superpixel generation, where we use a spatially regularized gradient to achieve a tunable trade-off between superpixel regularity and adherence to object boundaries. We quantitatively evaluate our method on the Berkeley segmentation database and show that we outperform the best linear-time state-of-the art method: Simple Linear Iterative Clustering (SLIC) [5]. We call the resulting superpixels "waterpixels."

B. Building Superpixels From Seeds

In the second step, the partition into superpixels is built from the seeds. Among the methods with linear complexity, there are two main strategies for this: Shortest Path Methods (Type 1) [11], [13]: these methods are based on region growing: they start from a set of seeds (points or regions) and successively extend them by incorporating pixels in their neighborhood according to a usually image dependent cost function until every pixel of the image plane has been assigned to exactly one superpixel. This process may or may not be iterated.

Shortest Distance Methods (Type 2) [5], [12]: these are iterative procedures inspired by the field of unsupervised learning, where at each iteration step, seeds (such as centroids) are calculated from the previous partition and pixels are then re-assigned to the closest seed (like for example the k-means approach).

Even though methods inspired by general clustering methods (type 2) seem appealing at first sight, in particular when they globally optimize a cost function, this class of methods does not guarantee connectivity of the superpixels for arbitrary choices of the pixel-seed distance (see [5], [12]). For instance, the distance metric proposed in [5] (a combination of Euclidean and grey level distance), leads to non-connected superpixels, which is undesirable. To solve this issue, a postprocessing step is necessary, consisting either in relabeling the image so that every connected component has its own label (see [12]), leading to a more irregular distribution of SP sizes and shapes, or in reassigning isolated regions to the closest and large enough Superpixel, as in [5], leading to non-optimality of the solution and an unpredictable number of superpixels. In addition, such postprocessing increases the computational cost and can turn out to be the most time-consuming step when the image contains numerous small objects/details compared to the size of the Superpixel.

On the contrary, methods based on region growing (type 1) inherently implement a "path-type" distance, where the distance between two pixels does not only depend on value and position of the pixels themselves, but on values and positions along the path connecting them. Type 1 methods imply connected superpixel regions, for which the number of superpixels is exactly the number of seeds.

C. Other Properties

It is generally accepted that a good superpixel-generation method should provide to the user total control over the number of resulting Superpixels. While this property is achieved by [11]- [14], some only reach approximatively this number because of post-processing (either by splitting too big superpixels, or removing small isolated superpixels as in [5]). Another parameter is the control on superpixels regularity in the trade-off between regularity and adherence to contours. Only [5] and [12] enable the user to weight the importance of regularity compared to boundary adherence, so it can be 160 adapted to the application.

161

As far as performance is concerned, one of the main 162 criteria is undoubtedly the complexity that the method 163 requires. Indeed, for Superpixels to be used as primitives for 164 further analysis such as classification, their computation should 165 neither take too long nor too much memory. This is the reason 166 why we focus on linear complexity methods. Among them, 167 SLIC appears to offer the best performance with regards to the 168 trade-off between adherence to boundaries and regularity [5].

169 Moreover, since its recent inception, this method has become 170 very popular in the computer vision community. We will 171 therefore use it as reference for the quantitative evaluation of 172 our method.

173

D. Superpixels and Watershed 174

In principle, the watershed transformation (see [15] for a 175 review) is well suited for SP generation: 3D applications [16], [17], in particular when shape regularity 187 of the elementary regions was not required.

188

Previous publications claimed that the watershed transfor-189 mation does not allow for the generation of spatially regular 190 SP [5], [11]. Recently, we and others [6], [18] have shown

A. Gradient and Cells Definition

Let f : D → V be an image, where D is a rectangular subset of Z 2 , and V a set of values, typically {0, . . . , 255} when f is a grey level image, or {0, . . . , 255} 3 for color images.

The first step consists in computing the gradient image g of the image f . The choice of the gradient operator depends on the image type, e.g. for grey level images we might choose a morphological gradient. This gradient will be used to choose the seeds (section III-B) and to build the regularised gradient (III-C).

For the definition of cells, we first choose a set of N points {o i } 1≤i≤N in D, called cell centers, so that they are placed on the vertices of a regular grid (a square or hexagonal one for example). Given a distance d on D, we denote by σ the grid step, i.e. the distance between closest grid points.

A Voronoi tesselation allows to associate to each o i a Voronoi cell . For each such cell, a homothety centered on o i with factor ρ (0 < ρ ≤ 1) leads to the computation of the final cell C i . This last step allows for the creation of a margin between neighbouring cells, in order to avoid the selection of markers too close from each other.

B. Selection of the Markers

As each cell is meant to correspond to the generation of a unique waterpixel, our method, through the choice of one marker per cell, offers total control over the number of SP, with a strong impact on their size and shape if desired.

First, we compute the minima of the gradient g. Each minimum is a connected component, composed of one or more pixels. These minima are truncated along the grid, i.e. pixels which fall on the margins between cells are removed.

Second, every cell of the grid serves to define a region of interest in the gradient image. The content of g in this very region is then analyzed to select a unique marker, as explained in the next paragraph.

For each cell, the corresponding marker is chosen among the minima of g which are present in this very cell.

If several minima are present, then the one with the highest surface extinction value [19] is used. We have found surface extinction values to give the best performances compared with volume and dynamic extinction values (data not shown).

It may happen that there is no minimum in a cell. This is an uncommon situation in natural images. In such cases, we must add a marker for the cell which is not a minimum of g, in order to keep regularity. One solution could be to simply choose the center of the cell; however, if this point falls on a local maximum of the gradient g, the resulting SP may coincide with the maximum region and therefore be small in size (leading to a larger variability in size of the SP). We propose instead to take, as marker, the flat zone with minimum value of the gradient inside this very cell.

In both cases (i.e. either there exists at least one minimum in the cell or there is not), the selected marker has to be composed of a unique connected component to ensure regularity and connectivity of the resulting superpixel. However, it might not be the case, respectively if more than one minimum

∀ p ∈ D, d Q ( p) = 2 σ min i∈[1,N] d( p, q i ) ( 1
)
where σ is the grid step defined in the previous section. The normalization by σ is introduced to make the regularization independent from the chosen SP size.

We have studied two possible choices of the q i . The first one is to choose them equal to the markers: q i = M i . Resulting waterpixels are called m-waterpixels. The second one consists in setting them at the cell centers: q i = o i , which leads to c-waterpixels. We have found that the first gives the best adherence to object boundaries, while the second produces more regular superpixels.

The spatially regularized gradient g reg is defined as follows:

g reg = g + kd Q (2)
where g is the gradient of the image f , d Q is the distance function defined above and k is the spatial regularization parameter, which takes its values within + . The choice of k is application dependent: when k equals zero, no regularization of the gradient is applied; when k → ∞, we approach the Voronoi tessellation of the set {q i } 1≥i≥N in the spatial domain.

In the final step, we apply the watershed transformation on the spatially regularized gradient g reg , starting the flooding from the markers {M i } 1≤i≤N , so that an image partition {s i } 1≤i≤N is obtained. The s i are the resulting waterpixels.

IV. EXPERIMENTS

In order to evaluate waterpixels, the proposed method has been applied on the Berkeley segmentation database [20] and benchmarked against the state-of-the-art. This database is divided into three subsets, "train", "test" and "val", containing respectively 200, 200 and 100 images of sizes 321 × 481 or 481 × 321 pixels. Approximately 6 human-annotated ground-truth segmentations are given for each image. These ground-truth images correspond to manually drawn contours.

A. Implementation

We have found that it is beneficial to pre-process the images from the database using an area opening followed by an area closing, both of size σ 2 /16 (where σ is the chosen step size of the regular grid). This operation efficiently removes details which are clearly smaller than the expected waterpixel area and which should therefore not give rise to a superpixel contour.

The Lab-gradient is adopted here in order to best reflect our visual perception of color differences and hence the pertinence of detected objects. The margin parameter ρ, described in III-A, is set to 2 3 .

The cell centers correspond to the vertices of a square or an hexagonal grid of step σ . The grid is computed in one pass over the image, by first calculating analytically the coordinates of the set of pixels belonging to each cell and then assigning to them the label of their corresponding cell. We will display the results for the hexagonal grid, as hexagons are more isotropic than squares. Interestingly, they also lead to a better quantitative performance, which was intuitively expected.

The implementation of the waterpixels was done using the Simple Morphological Image Library (SMIL) [21]. SMIL is a Mathematical Morphology library that aims to be fast, lightweight and portable. It brings most classical morphological operators re-designed in order to take advantage of recent computer features (SIMD, parallel processing, …) to allow handling of very large images and real time processing. As a gradient-based approach, the quality of the watershed is dependant on the borders contrast. If we look at the contours of objects missed by waterpixels, we see that it is due to the 368 weakness of the gradient, as illustrated in Figure 4.

B. Qualitative Analysis

369

C. Evaluation Criteria

370

SP methods produce an image partition {s i } 1≤i≤N . In order 371 to compute the SP borders, we use a morphological gradient 372 with a 4 neighborhood. Note that the resulting contours are 373 two pixels wide. To this set S c , we add the one pixel wide 374 image borders S b . The final set is denoted C. The ground 375 truth image corresponding to the contours of the objects to be 376 segmented, provided in the Berkeley segmentation database, 377 is called GT .

378

In superpixel generation, we look for an image decomposi-379 tion into regular regions that adhere well to object boundaries. 380 We propose to use three measures to evaluate this trade-381 off, namely boundary-recall, contour density and average 382 mismatch factor, as well as computation time.

383

There are two levels of regularity: (1) the number of pixels 384 required to describe the SP contours, which can be seen as a 385 measure of complexity of individual SP, and (2) the similarity 386 in size and shape between SP.

387

The first property is evaluated by the Contour Density, 388 which is defined as the number of SP contour pixels divided 389 by the total number of pixels in the image:

390 C D = 1 2 |S c | + |S b | |D| (3) 391
Note that |S c | is divided by 2 since contours are 392 two-pixel-wide.

393

The second property, i.e. similarity in size and shape, is 394 evaluated by an adapted version of the mismatch factor [22]. 395 The mismatch factor measures the shape and size dissimilarity 396 between two regions. Given two sets, A and B, the mismatch 397 factor m f between them is defined as:

398 m f (A, B) : = |A ∪ B \ A ∩ B| |A ∪ B| 399 = 1 - |A ∩ B| |A ∪ B| (4) 400
The mismatch factor and the Jaccard index thus sum to 401 one. Aiming to measure the superpixel regularity, we adapted 402 the mismatch factor to estimate the spread of size and shape 403 distribution. Hence, the average mismatch factor M F is 404 proposed as:

405 M F = 1 N N i=1 m f (s * i , s * ) (5) 406
where s * i is the centered version of superpixel s i , and s * is 407 the average centered shape of all superpixels. The complete 408 definition of the average mismatch factor is given in Appendix. 409

Note that although compactness is sometimes used in super-410 pixels evaluation (see [23]), it is a poor measurement for 411 region regularity. For example, perfectly-rectangular regions 412 are regular but not compact (because they are different from 413 discs). Waterpixels can in principle tend towards differently 414 shaped superpixels (rectangles, hexagons or other), depending 415 on the grid and the regularization function used. Since the 416 average mismatch factor compares each superpixel against an 417 426 While precision cannot be directly used in the context of over-segmentations, boundary-recall has to be, in this particular case of superpixels, interpreted with caution. Indeed, as noted also by Kalinin and Sirota [24], very tortuous contours systematically lead to better performances: because of their higher number, SP contour pixels have a higher chance of matching a true contour, increasing artificially the boundaryrecall. Hence, we propose to always consider the trade-off between boundary-recall and contour density to properly evaluate the adherence to object boundaries, penalizing at the same time the cost in pixels to describe SP contours.

D. Quantitative Analysis and Comparison With State-of-the-Art

In this paragraph, we will use m-waterpixels and denote them directly as "waterpixels" for the sake of simplicity.

During the design of the algorithm, we used intermediate results from the train and test subsets of the Berkeley database. Therefore, we report the results obtained for the validation subset ("val"), which contains 100 images. Results for boundary-recall, average mismatch factor and contour density are averaged for this subset and shown in Figure 5. than for SLIC (k = 15). This means that for the same number 474 of superpixels, waterpixels contours are shorter than SLIC 475 contours, which is partly explained by less tortuous contours. 476 We are currently working on a new implementation of minima computation/evaluation which would be less dependent on the number of superpixels.

To conclude this section, waterpixels are generally faster to compute than SLIC superpixels, and they are at least as performant in the trade-off between adherence to object boundaries and regularity in shape and size, while using much less pixels to describe their contours.

V. DISCUSSION AND PERSPECTIVES

We have shown that waterpixels produce competitive results with respect to the state-of-the-art. These advantages are valuable in the classification/detection/segmentation pipeline, where superpixels play the part of primitives. Moreover, there is one major difference in the construction of the algorithm: the SLIC approach does not impose any connectivity constraint. The resulting superpixels are therefore not necessarily connected, which requires some ad hoc postprocessing step. In contrast, waterpixels are connected by definition, and the connectivity constraint is actually implemented in the distance used.

The proposed approach is gradient-based. Standard methods can be used to compute this gradient, or a specific gradient computation method can be designed for a given application. In any case, this offers flexibility to waterpixels. One limitation though is the quality of the signal in such a gradient image. As seen in 4, alteration by noise or insufficiently contrasted contours may lead to the prevalence of regularity over adherence to object boundaries. If filtering steps are usually enough to deal with noise and remove non pertinent small details, parameter values have to be optimized for each database. Future work will aim at overcoming this limitation by adding a learning step of optimal filtering values for specific databases.

The general design of waterpixels offers many prospects. Among them, one promising field of improvement resides in the placement of markers, as they constitute the main degree of freedom of the method. We are currently investigating the possibility to select the markers in an optimal manner, for example by formulating the marker placement as a p-dispersion problem (see [25]) in an augmented space.

The speed of waterpixels contributes to expanding their 551 possible applications. For example, it could be interest-

Fig. 1 .

 1 Fig. 1. Superpixels illustration. The original image comes from the Berkeley segmentation database. (a) Original image. (b) Waterpixels.

Fig. 2 .

 2 Fig. 2. Illustration of waterpixels generation: (a): original image; (b) corresponding Lab gradient; (c): selected markers within the regular grid of hexagonal cells (step σ = 40 pixels); (d): distance function to markers; (g): distance function to cell centers; (e) and (h): spatially regularized gradient respectively with distance functions to selected markers (d) and to cell centers (g); (f) and (i): Resulting waterpixels obtained by respectively applying the watershed transformation to (e) and (h), with markers (c).

Figure 3

 3 Figure 3 shows various images from the Berkeley segmentation database and their corresponding waterpixels (m-waterpixels and c-waterpixels, hexagonal and square grids, different steps). Figures 3.b and 3.c (zooms of original image presented in 3.a for m-waterpixels and c-waterpixels respectively) show the influence of the regularization parameter k (0, 4, 8, 16) for an homogeneous (blue sky) and a textured (orange rock) regions. As expected, when k → ∞, m-waterpixels tend towards the Voronoi tessellation of the markers, while c-waterpixels approach the regular grid of hexagonal cells. Both show good adherence to object boundaries, as shown in Figures 3.d, 3.e, 3.f. Of course, enforcing

  , 3.e, 3.f. Of course, enforcing regularity decreases the adherence to object boundaries (see the zoom in Figure3.f for k = 16). One advantage of waterpixels is that the user can choose the shape (and size) of resulting superpixels depending on the application requisites.

Figure 3 .

 3 Figure 3.d, for example, presents waterpixels for hexagonal (second and third columns) and square (fourth column) grids.

Fig. 3 .

 3 Fig. 3. Illustrations of waterpixels on the Berkeley segmentation database: All waterpixels images are computed with an hexagonal grid with step σ = 30 pixels and a regularization parameter k = 8, unless otherwise specified. (a): original image (middle) with corresponding m-waterpixels (left) and c-waterpixels (right). σ = 25 pixels, k = 16. (c): zooms of m-waterpixels (a) for k = 0, 4, 8, 16. (c): zooms of c-waterpixels (a) for k = 0, 4, 8, 16. (d): original image -m-wat. -c-wat. -m-wat. with square grid and σ = 40 pixels. (e): original image -m-wat. -c-wat. -zoom of c-wat.. (f): original image -m-wat. -c-wat. -zoom of m-wat. with k = 16.image dependent template, this measure is more appropriate 418

419

  To quantify the adherence to object boundaries, a classical 420 measure used in the literature is the boundary-recall (BR).421Boundary-recall is defined as the percentage of ground-truth 422 contour pixels GT which fall within strictly less than 3 pixels 423 from superpixel boundaries C:

Fig. 4 .

 4 Fig. 4. Contours missed by waterpixels: (a): original image from the Berkeley segmentation database. (b): m-waterpixels with step = 27 and k = 10. (c): c-waterpixels with step = 27 and k = 10. (d), (e), (f): zoom of (a), (b), (c) respectively. (g): zoom of the non-regularized gradient image. (h) and (i): reached (green) and missed (red) contours, respectively by m-waterpixels and c-waterpixels.

  Blue and red curves correspond to varying regularization parameters k and k respectively for waterpixels and SLIC. The values for k and k have been chosen such that they cover a reasonable portion of the regularization space between no regularization (k = 0) and a still acceptable level of regularization.

Figure 5 (

 5 Figure 5(a) shows contour density against boundary-recall for waterpixels and SLIC. The ideal case being the lowest contour density for the highest boundary-recall, we can see that the trade-off between both properties improves for decreasing regularization, as expected. On the other hand, SLIC shows another behavior: the trade-off improves, then gets worse with regularization. At any rate, it is important to note that waterpixels achieves a better "best" trade-off than SLIC (see waterpixel k = 0 and SLIC k = 15). Besides, this observation is valid for the whole family of waterpixel-methods as the zero-value regularization does not take into account d Q . In order to do a fair comparison between waterpixels and SLIC over all criteria, we choose corresponding curves in the trade-off contour density/boundary-recall, i.e. waterpixels with k = 8 and SLIC with k = 15, and compare this couple for the other criteria.

Figure 5 (

 5 Figure 5(b) shows that, for a given number of superpixels, contour density of waterpixels is more stable and most of the time lower than SLIC when varying regularization. More particularly, contour density is lower for waterpixels (k = 8)

Fig. 5 .

 5 Fig. 5. Benchmark: performance comparison between waterpixels and SLIC. (a) Contour Density against Boundary-recall. (b) Contour Density against Number of Superpixels. (c) Mismatch factor against Boundary-recall.

Figure 5 (

 5 c) shows average mismatch factor against 477 boundary-recall for waterpixels and SLIC. We can see that 478 the curves for waterpixels with k = 8 and SLIC with k = 15 479 are here again close to each other.

Fig. 6 .

 6 Fig. 6. Comparison between Waterpixels and SLIC superpixels for σ = 25 pixels on a zoom of an image from the Berkeley segmentation database. (a) SLIC k = 15. (b) Waterpixels k = 8.

Fig. 7 .

 7 Fig. 7. Computation time comparison with images of the Berkeley database.

) 590 where 1 i

 1 is the indicator function of s * i . Thus, image S cor-591 responds to the summation image of all centered superpixels. 592 Let furthermore μ A = 1/n N i=1 |s i | be the average area of 593 the considered superpixels, and let S t be the threshold of S at 594 level t: S t (x) = {x p ∈ D S(x p ) ≥ t}.

18 S

 18 ) degree in electronic engineering and the M.Eng. 720 degree in industrial automation from the Univer-721 sidad Nacional de Colombia, Manizales-Colombia, 722 in 2008 and 2011, respectively. He is currently 723 pursuing the Ph.D. degree in automatics with the 724 Universidad Nacional de Colombia. He has been a 725 Research Assistant with the Signal Processing and 726 Recongnition Group since 2008. His current research 727 interests include machine learning and signal and 728 image processing. 729 Théodore Chabardes received the degree from the 730 Engineering School, ESIEE Paris, France, in 2014, 731 as an Engineer specialized in computer science. 732 He is currently pursuing the Ph.D. degree with 733 the Centre of Mathematical Morphology, School of 734 Mines, Paris, France. His research interests include 735 mathematical morphology, image segmentation, and 736 software optimization. 737 Thomas Walter received the Diploma degree in 738 electrical engineering from Saarland University, 739 Germany, and the Ph.D. degree in mathematical 740 morphology from Mines ParisTech, France. He held 741 a post-doctoral position with the European Molecu-742 lar Biology Laboratory, Heidelberg, Germany. He is 743 currently a Team Leader in bioimage informatics 744 with the Centre for Computational Biology, Mines 745 ParisTech, and a member of the Bioinformatics Unit 746 with the Curie Institute, Paris. His most visible 747 scientific contributions have been in the field of 748 bioimage informatics, and in particular, in high content screening. He has 749 pioneered methods in the field of cellular phenotyping and phenotypic 750 clustering from live cell-imaging data. He was involved in the first genome-751 wide screen by live cell imaging in a human cell line, and co-develops the 752 open-source software cellcognition. 753 Etienne Decencière received the Engineering degree 754 and the Ph.D. degree in mathematical morphology 755 from MINES ParisTech, France, in 1994 and 1997, 756 respectively, and the Habilitation à Diriger des 757 Recherches from Jean Monnet University, in 2008. 758 He holds a research fellow position with the Centre 759 for Mathematical Morphology, MINES ParisTech, 760 where he leads several academic and industrial 761 research projects. His main research interests are in 762 mathematical morphology, image segmentation, and 763 biomedical applications. into homogeneous regions with enforced regularity and adherence 3 to object boundaries. Methods to generate these superpixels have 4 gained substantial interest in the last few years, but only a few 5 have made it into applications in practice, in particular because 6 the requirements on the processing time are essential but are not 7 met by most of them. Here, we propose waterpixels as a general 8 strategy for generating superpixels which relies on the marker 9 controlled watershed transformation. We introduce a spatially 10 regularized gradient to achieve a tunable tradeoff between the 11 superpixel regularity and the adherence to object boundaries. 12 The complexity of the resulting methods is linear with respect 13 to the number of image pixels. We quantitatively evaluate our 14 approach on the Berkeley segmentation database and compare I. INTRODUCTION UPERPIXELS (SP) are regions resulting from 19 a low-level segmentation of an image and are typically 20 used as primitives for further analysis such as detection, 21 segmentation, and classification of objects (see Figure 1

22 for27 1 )

 221 an illustration). The underlying idea is that this first 23 low-level partition alleviates the computational complexity of 24 the following processing steps and improves their robustness, 25 as not single pixel values but pixel set features can be used. 26 Superpixels should have the following properties: homogeneity: pixels of a given SP should present 28 similar colors or gray levels; 29 Manuscript received December 18, 2014; revised April 17, 2015; accepted June 15, 2015. The work of T. Walter was supported by the European Community through the Seventh Framework Programme (FP7/2007-2013) within the Systems Microscopy under Grant 258068. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Yonggang Shi.

Fig. 1 .

 1 Fig. 1. Superpixels illustration. The original image comes from the Berkeley segmentation database. (a) Original image. (b) Waterpixels.

176 1 )

 1 It gives a good adherence to object boundaries when 177 computed on the image gradient.

178 2 ) 3 )

 23 It allows to control the number and spatial arrangement 179 of the resulting regions through the choice of markers. 180 The connectivity of resulting regions is guaranteed and 181 no postprocessing is required.

182 4 )

 4 It offers linear complexity with the number of pixels in 183 the image. 184 Indeed, it has been used to produce low-level segmentations 185 in several applications, including computation intensive 186

191 1 ) 2 ) 3 ) 4 ) 5 )

 12345 that in principle the watershed transformation can be applied 192 to SP generation. 193 Here, we introduce waterpixels, a family of methods based 194 on the watershed transformation to compute superpixels. 195 III. WATERPIXELS 196 As most watershed-based segmentation methods, 197 waterpixels are based on two steps: the definition of 198 markers, from which the flooding starts, and the definition of 199 a gradient (the image to be flooded). We propose to design 200 these steps in such a way that regularity is encouraged. 201 A waterpixel-generation method is characterized by the 202 following steps: 203 Computation of the gradient of the image; 204 Definition of regular cells on the image, centered on the 205 vertices of a regular grid; 206 Selection of one marker per cell; 207 Spatial regularization of the gradient with the help of a 208 distance function; 209 Application of the watershed transformation on the 210 regularized gradient defined in step 4 from the markers 211 defined in step 2. 212 These steps are illustrated in figure 2 and developed in the 213 next paragraphs.

  214

Fig. 2 .

 2 Fig. 2. Illustration of waterpixels generation: (a): original image; (b) corresponding Lab gradient; (c): selected markers within the regular grid of hexagonal cells (step σ = 40 pixels); (d): distance function to markers; (g): distance function to cell centers; (e) and (h): spatially regularized gradient respectively with distance functions to selected markers (d) and to cell centers (g); (f) and (i): Resulting waterpixels obtained by respectively applying the watershed transformation to (e) and (h), with markers (c).

  288

Figure 3

 3 Figure 3 shows various images from the Berkeley segmentation database and their corresponding waterpixels (m-waterpixels and c-waterpixels, hexagonal and square grids, different steps). Figures 3.b and 3.c (zooms of original image presented in 3.a for m-waterpixels and c-waterpixels respectively) show the influence of the regularization parameter k (0, 4, 8, 16) for an homogeneous (blue sky) and a textured (orange rock) regions. As expected, when k → ∞, m-waterpixels tend towards the Voronoi tessellation of the markers, while c-waterpixels approach the regular grid of hexagonal cells. Both show good adherence to object boundaries, as shown in Figures 3.d, 3.e, 3.f. Of course, enforcing regularity decreases the adherence to object boundaries (see the zoom in Figure 3.f for k = 16). One advantage of waterpixels is that the user can choose the shape (and size) of resulting superpixels depending on the application requisites.

Figure 3 .

 3 Figure 3.d, for example, presents waterpixels for hexagonal (second and third columns) and square (fourth column) grids.

Fig. 3 .

 3 Fig. 3. Illustrations of waterpixels on the Berkeley segmentation database: All waterpixels images are computed with an hexagonal grid with step σ = 30 pixels and a regularization parameter k = 8, unless otherwise specified. (a): original image (middle) with corresponding m-waterpixels (left) and c-waterpixels (right). σ = 25 pixels, k = 16. (c): zooms of m-waterpixels (a) for k = 0, 4, 8, 16. (c): zooms of c-waterpixels (a) for k = 0, 4, 8, 16. (d): original image -m-wat. -c-wat. -m-wat. with square grid and σ = 40 pixels. (e): original image -m-wat. -c-wat. -zoom of c-wat.. (f): original image -m-wat. -c-wat. -zoom of m-wat. with k = 16.image dependent template, this measure is more appropriate 418

419

  To quantify the adherence to object boundaries, a classical 420 measure used in the literature is the boundary-recall (BR).421Boundary-recall is defined as the percentage of ground-truth 422 contour pixels GT which fall within strictly less than 3 pixels 423 from superpixel boundaries C:424 B R = |{ p ∈ GT, d( p, C) < 3}| |GT |(6)425where d is the L 1 (or Manhattan) distance.

Fig. 4 .

 4 Fig. 4. Contours missed by waterpixels: (a): original image from the Berkeley segmentation database. (b): m-waterpixels with step = 27 and k = 10. (c): c-waterpixels with step = 27 and k = 10. (d), (e), (f): zoom of (a), (b), (c) respectively. (g): zoom of the non-regularized gradient image. (h) and (i): reached (green) and missed (red) contours, respectively by m-waterpixels and c-waterpixels.

  Blue and red curves correspond to varying regularization parameters k and k respectively for waterpixels and SLIC. The values for k and k have been chosen such that they cover a reasonable portion of the regularization space between no regularization (k = 0) and a still acceptable level of regularization.

Figure 5 (

 5 Figure 5(a) shows contour density against boundary-recall for waterpixels and SLIC. The ideal case being the lowest contour density for the highest boundary-recall, we can see that the trade-off between both properties improves for decreasing regularization, as expected. On the other hand, SLIC shows another behavior: the trade-off improves, then gets worse with regularization. At any rate, it is important to note that waterpixels achieves a better "best" trade-off than SLIC (see waterpixel k = 0 and SLIC k = 15). Besides, this observation is valid for the whole family of waterpixel-methods as the zero-value regularization does not take into account d Q . In order to do a fair comparison between waterpixels and SLIC over all criteria, we choose corresponding curves in the trade-off contour density/boundary-recall, i.e. waterpixels with k = 8 and SLIC with k = 15, and compare this couple for the other criteria.

Figure 5 (

 5 Figure 5(b) shows that, for a given number of superpixels, contour density of waterpixels is more stable and most of the time lower than SLIC when varying regularization. More particularly, contour density is lower for waterpixels (k = 8)

Fig. 5 .

 5 Fig. 5. Benchmark: performance comparison between waterpixels and SLIC. (a) Contour Density against Boundary-recall. (b) Contour Density against Number of Superpixels. (c) Mismatch factor against Boundary-recall.

Figure 5 (

 5 c) shows average mismatch factor against 477 boundary-recall for waterpixels and SLIC. We can see that 478 the curves for waterpixels with k = 8 and SLIC with k = 15 479 are here again close to each other.

Fig. 6 .

 6 Fig. 6. Comparison between Waterpixels and SLIC superpixels for σ = 25 pixels on a zoom of an image from the Berkeley segmentation database. (a) SLIC k = 15. (b) Waterpixels k = 8.

Fig. 7 .

 7 Fig. 7. Computation time comparison with images of the Berkeley database.

552

  ing to compute different sets of waterpixels, by changing 553 design options (different cells, gradients, grid steps, etc.), 554 and then use ensemble clustering methods to obtain a final 555 segmentation[26],[27]. 556 Last buy not least, waterpixels lead to the efficient construc-557 tion of hierarchical partitions based on superpixels. Indeed, 558 the computation of the watershed can produce at the same 559 time a segmentation and a hierarchy of partitions based on 560 that segmentation, with only minor overhead computation 561 times [28]-[30]. 562 VI. CONCLUSION 563 This paper introduces waterpixels, a family of methods 564 for computing regular superpixels based on the watershed 565 transformation. Both adherence to object boundaries and regu-566 larity of resulting regions are encouraged thanks to the choice 567 of the markers and the gradient to be flooded. Different

TABLE I RECAP

 I CHART OF EXISTING METHODS TO COMPUTE REGULAR SUPERPIXELS (n IS THE NUMBER OF PIXELS IN THE IMAGE; i IS THE NUMBER OF ITERATIONS REQUIRED; N THE NUMBER OF SUPERPIXELS). "WP" CORRESPONDS TO OUR METHOD, CALLED "WATERPIXELS"

TABLE I RECAP

 I CHART OF EXISTING METHODS TO COMPUTE REGULAR SUPERPIXELS (n IS THE NUMBER OF PIXELS IN THE IMAGE; i IS THE NUMBER OF ITERATIONS REQUIRED; N THE NUMBER OF SUPERPIXELS). "WP" CORRESPONDS TO OUR METHOD, CALLED "WATERPIXELS" are typically the cells or the centers of a regular grid. the computation of the SP. Type B methods may spend more time on finding appropriate seeds, but can therefore afford not to iterate the SP generation.

96

• Type B seeds depend on the content of the image 97 (compromise between a regular cover of the image plane 98 and an adaption to the contour).

99 • Type C seeds are initially image independent, then they 100 are iteratively refined to take into account the image 101 contents. 102 If the seed does not depend on the image, an iterative 103 refinement is usually preferable, and therefore more time 104 is spent on

1057-7149 © 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

The work of T. Walter was supported by the European Community through the Seventh Framework Programme (FP7/2007-2013) within the Systems Microscopy under Grant 258068.

AUTHOR QUERIES

AQ:1 = Please check whether the edits made in the financial section are OK. AQ:2 = Please confirm the current affiliation of all authors. AQ:3 = Please confirm the postal code for "MINES ParisTech, PSL Research University, Center for Mathematical Morphology, Universidad Nacional de Colombia, Centre for Computational Biology, Institut Curie, and Inserm." AQ:4 = Table I is not cited in body text. Please indicate where it should be cited. AQ:5 = Please provide the page range for ref. [1]. AQ:6 = Please provide the publisher name and location for ref. [2]. AQ:7 = Please provide the page range and also confirm the conference title for ref. [17]. AQ:8 = Please provide the organization, location, and report no. for ref. [21].

it against the state-of-the-art.

as well as their different behaviours in terms of regularity 483 (shape, size, tortuosity). i of s i is obtained by translating s i so that its barycenter is 586 the origin of the coordinates system.

484

E. Computation Time

587

The average shape s * of the {s i } is computed as follows.

588

Let first define function S: Finally, the mean mismatch factor of superpixels {s i } 1≤i≤N is: I is not cited in body text. Please indicate where it should be cited. AQ:5 = Please provide the page range for ref. [1]. AQ:6 = Please provide the publisher name and location for ref. [2]. AQ:7 = Please provide the page range and also confirm the conference title for ref. [17]. AQ:8 = Please provide the organization, location, and report no. for ref. [21].